Spatial-Temporal Evolution and Pattern Characteristics of Green Infrastructure in the Pisha Sandstone Area of the Middle Reaches of the Yellow River Basin, China

1. Introduction

As a distinctive lithological type, pisha sandstone is characterized by high porosity, loose structure, and weak cementation [1]. This inherent property causes it to disintegrate into mud when wetted and disperse as sand when exposed to wind [2]. resulting in a typical landscape of fragmented terrain and longitudinal gullies in its distribution area, with extremely low stability of the ecosystem substrate [3]. From a global perspective, pisha sandstone resources are scarce and concentrated, mainly distributed continuously in the middle reaches of the Yellow River, China, covering an area of approximately 16,700 km² across the junction of Shaanxi, Shanxi, and Inner Mongolia provinces [4]. It is the world’s largest and most intensely eroded concentrated pisha sandstone area, with an erosion modulus as high as 30,000–40,000 t/(km²·a) [5]. Although accounting for only 2% of the Yellow River Basin area, it contributes 1/4 of the coarse sediment (particle size >0.05 mm) deposited in the lower reaches of the Yellow River [6]. Consequently, it has become the core source area of coarse sediment in the Yellow River and the most challenging area for soil erosion control in the basin [7]. Previous studies and management practices in the pisha sandstone area have primarily focused on micro-scale processes such as erosion mechanisms [8,9], action principles [10], and the effects of vegetation and water systems on soil and water conservation [11,12,13,14]. While these efforts have deepened the understanding of local ecological processes and achieved certain management results, they lack systematic attention to regional patterns, spatial configurations, and macro-quantitative and qualitative issues. The traditional management model mainly adopts fragmented approaches such as point, line, and surface-based governance, making it difficult to achieve the overall optimization of the ecosystem. As a spatial ecological system integrating “patch-corridor-network”, Regional Green Infrastructure (RGI) can identify priority conservation areas, systematically connect fragmented ecological units, and construct a complete ecological protection network [15,16],providing a novel solution for soil erosion control in pisha sandstone areas. Relevant studies have also demonstrated that the optimization of RGI can weaken soil erosion intensity, reduce terrain fragmentation, and mitigate soil erosion [17,18].

GI was first proposed by the GI Task Force jointly organized by the U.S. Conservation Foundation and the U.S. Department of Agriculture Forest Service [19,20]. It constitutes an interconnected green space network composed of various functional landscape types, including natural and artificial ecological areas such as greenways, wetlands, forests, and native vegetation [21,22]. It emphasizes the connectivity and systematic integration of ecological elements[23], enabling the construction of comprehensive ecological networks at the regional scale[24] to achieve the holistic enhancement of ecological functions. GI can specifically connect fragmented grassland and woodland ecological patches formed by gully cutting and human development in the pisha sandstone area, restore soil and water conservation function corridors interrupted by mining activities and urban expansion, and enhance landscape connectivity. Through the protection of core patches, it can improve soil erosion resistance and water conservation capacity, strengthen the soil and water conservation service function of the ecosystem, and provide systematic support for regional ecological security. Therefore, in-depth research on the spatiotemporal evolution of GI in the pisha sandstone area and analysis of its pattern characteristics can not only make up for the deficiencies of existing studies in terms of systematic ecological patterns, spatial configurations, and macro-quantitative and qualitative issues in the pisha sandstone area but also provide a theoretical basis for territorial spatial planning, soil erosion control, and human habitat improvement in the pisha sandstone area. Meanwhile, it facilitates targeted and phased ecological space management, ensuring regional ecological security and overall ecosystem optimization, and providing scientific support for exploring sustainable development paths.

The development of patch-corridor-network frameworks is becoming increasingly comprehensive, with research primarily focusing on the spatiotemporal evolution and prediction of GI [25,26,27], the construction and optimization of ecological networks [28,29,30,31,32], and the evaluation of ecosystem services and ecological security patterns [33,34,35,36,37,38,39].For instance, scholars such as Ran Yue[40], Ana Medeiros[41], Wang Jiening[42], Duan Pengting[43] have examined land use changes under various scenarios to forecast future trends. Their analyses of GI’s spatiotemporal evolution reveal its characteristics over different periods, the development of core areas, and its ecological impact. Researchers like Gu Kangkang[44], Dimas Danar Dewa[45], Guan Huiyuan[46], Zhang Haoxiang[47], Lu Yujia[48] have identified ecological sources, resistance surfaces, and corridors to construct and optimize ecological networks. Meanwhile, Duy X. Tran[49], Matthew G.E.[50], Gao Yaling[51], Zerong Li[52] explored the spatiotemporal evolution of ecosystem services in different geographical environments and their interactions with ecological risks and landscape structures, emphasizing the importance of optimizing land use and landscape patterns to enhance ecosystem services.Based on the above studies, the research framework of “ecological source identification-resistance surface construction-ecological corridor extraction” has become a paradigm for GI research, and related research methods have shown a diversified trend. Its technical path can be summarized as follows: taking land use change as the research basis, identifying the laws of land use type conversion by analyzing its dynamic change characteristics; using Guidos software with an 8-neighborhood setting to identify GI types through MSPA and analyze their spatiotemporal evolution. Among them, the identification of ecological source areas—key patches supporting regional ecological processes and functions—has been realized through methods such as Morphological Spatial Pattern Analysis (MSPA)[53], Landscape Connectivity Analysis (LCA) [54], Ecosystem Service Value (ESV) assessment [55], Ecological Sensitivity Assessment (ESA) [56], and the Integrated Valuation of Ecosystem and Land use Tradeoffs (InVEST) model[57]. Studies have shown that combining landscape connectivity as a source screening index with the MSPA method [58,59,60] can improve the accuracy and rationality of ecological source selection, making it suitable for areas with severe landscape fragmentation. On this basis, landscape pattern indices such as the number of patches, patch density, largest patch index, landscape shape index, and aggregation degree [61] ,can be further introduced for quantitative analysis to comprehensively evaluate the distribution characteristics and spatial patterns of ecological sources.The construction of ecological resistance surfaces aims to form an evaluation carrier that comprehensively reflects the impact of spatial landscape heterogeneity on ecological flows [62,63,64]. It is common to use a combination of natural and anthropogenic factors to determine the resistance levels and corresponding resistance values of different influencing factors. The construction of ecological corridors mainly adopts methods such as the Least Cost Path (LCP) model [65], Minimum Cumulative Resistance (MCR) model [66], Gravity model [67], and circuit theory [68]. Among them, the MCR model has significant advantages in calculating the coupling of natural and human disturbances and can be paired with the Gravity model to identify important corridors [69,70], enabling in-depth analysis of network structural characteristics.

Therefore, this study selected the pisha sandstone area in the middle reaches of the Yellow River as the research area, based on land use data from 2003, 2013, and 2023 to provide objective data support for the identification and analysis of GI. Combined with MSPA and landscape connectivity analysis, ecological sources in the three periods were identified; the landscape index method was used to analyze changes in the landscape pattern of GI; the MCR model and Gravity model were employed to screen ecological corridors in the three periods, thereby constructing the GI network in the study area and evaluating its connectivity. Ultimately, the spatiotemporal evolution laws and pattern characteristics of GI in the study area over the past two decades were revealed (Figure 1). At the theoretical level, based on the special geomorphological and ecological vulnerability characteristics of the pisha sandstone area, this study systematically identifies and analyzes the regional GI, contributing to research on macro-ecological patterns, spatial configurations, and quantity-quality synergy in the pisha sandstone area. At the practical level, closely aligned with the strategy of ecological protection and high-quality development of the Yellow River Basin, this study clarifies the spatiotemporal evolution of GI, identifies key ecological sources and important ecological corridors, and determines spatial conservation priorities. It not only promotes the transformation from traditional fragmented governance to network-coordinated ecological restoration but also effectively addresses prominent regional issues such as soil erosion control and landscape fragmentation restoration, realizing the precise allocation of ecological protection resources.

2. Materials and Methods

2.1. Study Area

The pisha sandstone area in the middle reaches of the Yellow River is located at the junction of Shaanxi, Shanxi, and Inner Mongolia, between 38°10’-40°10’N latitude and 108°45’-111°31’E longitude, covering a total area of 16,700 km² (Figure 2).The terrain slopes from northwest to southeast, with an altitude ranging from 1000 to 1800 m. The topography and landscape are fragmented with longitudinal gullies. The area belongs to a typical semi-arid temperate continental climate, with an average annual temperature of 6-9°C and annual precipitation of approximately 280-400 mm, concentrated in July-September, accounting for 71.2% of the total annual precipitation, mostly in the form of short-duration heavy rains. This precipitation characteristic makes summer a high-risk period for soil erosion in the area. The rugged terrain and exposed bedrock pose significant challenges to soil and water conservation and management.

2.2. Data Sources

Land use data for the study area in 2003, 2013, and 2023 were obtained from the Chinese Land Cover Dataset (CLCD) developed by the research team of Professors Yang Jie and Huang Xin from Wuhan University using Google Earth Engine (https://zenodo.org/records/5816591), with a spatial resolution of 30 m. Landsat series remote sensing images with a 30 m resolution from the U.S. Geological Survey (USGS; http://espa.cr.usgs.gov) in the same periods were used as the primary data source to verify the accuracy of land use classification in the CLCD data through visual interpretation sampling. After verification, the data met the accuracy requirements for land use data in this study, ensuring the reliability of subsequent analysis results. Using ArcGIS software, combined with the shapefile of the pisha sandstone area, mask extraction was performed to obtain land use data within the pisha sandstone region (Figure 3). Monthly Normalized Difference Vegetation Index (NDVI) data were sourced from the MOD13A3 dataset regularly released by NASA (https://www.earthdata.nasa.gov/). Annual average values were calculated to generate the required annual NDVI raster data, with a spatial resolution of 1 km and coordinates in WGS1984. Digital Elevation Model (DEM) data were obtained from GEBCO’s 2024 global dataset covering land elevation and ocean depth (https://www.gebco.net/data_and_products/gridded_bathymetry_data/). Slope data were derived from elevation data using the Surface Analysis tool in the Spatial Analyst module of ArcGIS. Due to the lack of 2003 road network data, 2000 road network data from the National Basic Geographic Information Center of Tianditu were used as a substitute. Road data for 2013 and 2023 were sourced from OpenStreetMap.

2.3. Research Methods

2.3.1. Land Use Change Analysis Methods

The single land use dynamic degree refers to the magnitude of change in a specific land use type per unit time, which can reflect the rate of quantitative change in a particular land-use type at the micro-scale and compare regional differences in land-use change [71]. The formulas for calculation are as follows:

$$\mathrm{K}={\displaystyle \frac{{\mathrm{U}}_2-{\mathrm{U}}_1}{{\mathrm{U}}_1}}\times {\displaystyle \frac{1}{{\mathrm{T}}_2-{\mathrm{T}}_1}}\times 100$$

(1)

In the formula, $\mathrm{K}$ is the single land use dynamic degree ; ${\mathrm{U}}_1$ and ${\mathrm{U}}_2$ denote the land area of a particular land category at the beginning and end of the study period, respectively; ${\mathrm{T}}_1$ and ${\mathrm{T}}_2$ represent the first and last years of the study period, respectively.

To comprehensively understand the conversion between land use types, the land use transfer matrix was adopted. The land use transfer matrix is mainly used to analyze the mutual conversion of different land use types in different periods, revealing the direction and quantity of conversion between different land use types [72]. The formulas for calculation are as follows:

$${\mathrm{S}}_{\mathrm{i}\mathrm{j}}=\left[\begin{array}{ccc}{\mathrm{S}}_{11}& \dots & {\mathrm{S}}_{1\mathrm{n}}\\ ⋮& \ddots & ⋮\\ {\mathrm{S}}_{\mathrm{n}1}& \dots & {\mathrm{S}}_{\mathrm{n}\mathrm{n}}\end{array}\right]$$

(2)

In the formula: ${\mathrm{S}}_{\mathrm{i}\mathrm{j}}$ indicates the land area converted from land use type $\mathrm{i}$ to land use type $\mathrm{j}$ during the study period; $\mathrm{i}$=$\mathrm{j}$ represents the area where land use type $\mathrm{i}$ remained unchanged throughout the study; $\mathrm{n}$ refers to the total number of land use types considered.

2.3.2. Identification of Ecological Sources Based on MSPA

Using the reclassification tool in ArcMap, grassland, forest land, and water bodies in the study area were assigned a value of 2 as the foreground, and other land use types were assigned a value of 1 as the background. After a series of preprocessing steps, raster data of the study area were obtained [73]. Using Guidos software with an 8-neighborhood setting, MSPA analysis was performed to classify the foreground into 7 landscape types (Table 1). The edge effect width needs to consider factors such as the study area’s conditions and target species. Existing studies have shown that corridors with a width of 30-60 m can basically meet the migration and diffusion needs of plants and animals [74,75,76,77], and corridors within this width range can maintain river ecosystems and prevent soil erosion [78,79]. By comparing the current scale of the study area with edge widths of 10 m, 30 m, 60 m, 90 m, and 120 m, it was found that the 30 m precision achieved relatively comprehensive landscape type identification, and the core area patch size met the research requirements. Therefore, 30 m was determined as the edge effect width.

Based on the ecological significance of MSPA elements, core areas in the study region were identified as the basic ecological source areas. To further ensure the connectivity function of these source areas, landscape connectivity index analysis was introduced [80]. The Integral Index of Connectivity (IIC) and Probability of Connectivity (PC) were selected as core metrics. Conefor Sensinode software was used to calculate the IIC and PC values of each core area patch,use equal weights to calculate the two types of indexes. Finally, ecological patches with a synthesized index value greater than 1 were selected as core ecological sources, identifying those with significant impacts on landscape connectivity [81]. The formulas are as follows:

$$\mathrm{P}\mathrm{C}={\displaystyle \frac{\sum _{\mathrm{i}=1}^{\mathrm{n}}\sum _{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{a}}_{\mathrm{i}}\times {\mathrm{a}}_{\mathrm{j}}\times {\mathrm{P}}_{\mathrm{i}\mathrm{j}}^{\ast }}{{\mathrm{A}}^2}}$$

(3)

$$\mathrm{I}\mathrm{I}\mathrm{C}={\displaystyle \frac{\sum _{\mathrm{i}-1}^{\mathrm{n}}\sum _{\mathrm{j}-1}^{\mathrm{n}}\left[\left({\mathrm{a}}_{\mathrm{i}}{\mathrm{a}}_{\mathrm{j}}\right)/\left(1+\mathrm{n}{\mathrm{l}}_{\mathrm{i}\mathrm{j}}\right)\right]}{{\mathrm{A}}_{\mathrm{L}}^2}}$$

(4)

In this context, $\mathrm{n}$ signifies the number of patches, while ${\mathrm{a}}_{\mathrm{i}}$ and ${\mathrm{a}}_{\mathrm{j}}$ represent the areas of Patch $\mathrm{i}$ and Patch $\mathrm{j}$, respectively; $\mathrm{A}$ is the total landscape area; ${\mathrm{P}}_{\mathrm{i}\mathrm{j}}^{\ast }$ indicates the maximum potential for diffusion between Patch $\mathrm{i}$ and Patch $\mathrm{j}$; and the values of $\mathrm{I}\mathrm{I}\mathrm{C}$ and $\mathrm{P}\mathrm{C}$ are proportional to landscape connectivity.

2.3.3. Spatiotemporal Analysis of GI Landscape Patterns Using Landscape Pattern Index Methods

Fragstats software was used to calculate landscape pattern indices. Combining the ecological significance of landscape pattern indices, six indices were selected: Number of Patches (NP), Largest Patch Index (LPI), Landscape Shape Index (LSI), Contagion Index (CONTAG), Shannon’s Diversity Index (SHDI), and Shannon’s Evenness Index (SHEI) to quantitatively analyze the spatiotemporal changes in the GI landscape pattern in the study area from 2003 to 2023.

2.3.4. Potential Corridor Simulation Based on the MCR Model

The MCR model combines quantitative analysis with spatial visualization capabilities. It identifies optimal pathways for species migration and energy exchange by calculating resistance values between ecological source areas and target areas. Currently, it is the most widely used, convenient, and rapid landscape connectivity assessment model [82,83]. As illustrated in Equation：

$$\mathrm{M}\mathrm{C}\mathrm{R}={\mathrm{f}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\sum _{\mathrm{j}=\mathrm{n}}^{\mathrm{i}=\mathrm{m}}{\mathrm{D}}_{\mathrm{i}\mathrm{j}}\times {\mathrm{R}}_{\mathrm{i}}$$

(5)

In the equation, $\mathrm{M}\mathrm{C}\mathrm{R}$ is the minimum cumulative resistance value; ${\mathrm{f}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ is a function indicating that the MCR is positively correlated with ecological processes; ${\mathrm{D}}_{\mathrm{i}\mathrm{j}}$ is the distance from ecological source $\mathrm{j}$ to landscape unit $\mathrm{i}$; $\mathrm{R}$ is the resistance coefficient of landscape unit $\mathrm{i}$ to the expansion of ecological sources.

The identification of source areas and the construction of ecological resistance surfaces are key to simulating minimum resistance corridors using the MCR model, and the selection of appropriate resistance factors is crucial for constructing resistance surfaces. Combining the general situation of the pisha sandstone area in the middle reaches of the Yellow River and existing research foundations, five resistance factors were selected: land use type, elevation, slope, distance to roads, and vegetation coverage [84] to construct a comprehensive resistance surface (Table 2). Discrete data such as land use types were standardized (assigned values 1-5, followed by positive standardization)[85].Continuous data such as slope, elevation, vegetation cover, roads, and rivers underwent positive and inverse range standardization. Both approaches employed dimensionless data transformation. For instance, steeper slopes indicate greater migration difficulty, thus receiving positive standardization values. Conversely, higher vegetation cover facilitates migration, resulting in negative standardization values [86].The formulas are as follows:

For positive indicators:

$$\overline{\mathrm{X}}={\displaystyle \frac{\mathrm{X}-{\mathrm{X}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}{{\mathrm{X}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathrm{X}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}}$$

(6)

For negative indicators:

$$\overline{\mathrm{X}}={\displaystyle \frac{X_{max}-X}{X_{max}-X_{min}}}$$

(7)

Among these: $\overline{\mathrm{X}}$ is the normalised indicator value,$\mathrm{X}$ is the original indicator data, ${\mathrm{X}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ and ${\mathrm{X}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ are the maximum and minimum values of the indicator, respectively.

Finally, the weights of resistance factors were determined using the Analytic Hierarchy Process (AHP), and the comprehensive resistance surface was obtained using the raster calculator in ArcGIS 10.8 software to provide cost data for the MCR model.

2.3.5. Identification and Evaluation of Key Corridors Based on Gravity Models

The Gravity model quantitatively evaluates the interaction intensity between source areas by calculating the mutual gravitational force between patches. Higher interaction intensity indicates lower resistance to the exchange of materials, energy, and information between source areas, and the corridors connecting these source areas are more important [87], based on which important ecological corridors can be extracted. The formula is as follows:

$${\mathrm{G}}_{\mathrm{i}\mathrm{j}}={\displaystyle \frac{{\mathrm{N}}_{\mathrm{i}}{\mathrm{N}}_{\mathrm{j}}}{{\mathrm{D}}_{\mathrm{i}\mathrm{j}}^2}}={\displaystyle \frac{\left[\frac{1}{{\mathrm{P}}_{\mathrm{i}}}\times \mathrm{l}\mathrm{n}\left({\mathrm{S}}_{\mathrm{i}}\right)\right]\left[\frac{1}{{\mathrm{P}}_{\mathrm{j}}}\times \mathrm{l}\mathrm{n}\left({\mathrm{S}}_{\mathrm{j}}\right)\right]}{{\left(\frac{{\mathrm{L}}_{\mathrm{i}\mathrm{j}}}{{\mathrm{L}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}\right)}^2}}={\displaystyle \frac{{\mathrm{L}}_{\mathrm{m}\mathrm{a}\mathrm{x}}^2\mathrm{l}\mathrm{n}\left({\mathrm{S}}_{\mathrm{i}}\right)\mathrm{l}\mathrm{n}\left({\mathrm{S}}_{\mathrm{j}}\right)}{{\mathrm{L}}_{\mathrm{i}\mathrm{j}}^2{\mathrm{P}}_{\mathrm{i}}{\mathrm{P}}_{\mathrm{j}}}}$$

(8)

In the equation, ${\mathrm{G}}_{\mathrm{i}\mathrm{j}}$ denotes the interaction force between patches $\mathrm{i}$ and $\mathrm{j}$, ${\mathrm{N}}_{\mathrm{i}}$ and ${\mathrm{N}}_{\mathrm{j}}$ are the weight values of the two patches, ${\mathrm{D}}_{\mathrm{i}\mathrm{j}}$ is the standardized value of the potential corridor resistance between patches $\mathrm{i}$ and $\mathrm{j}$, ${\mathrm{P}}_{\mathrm{i}}$ and ${\mathrm{P}}_{\mathrm{j}}$ are the resistance values of patches $\mathrm{i}$ and $\mathrm{j}$, ${\mathrm{S}}_{\mathrm{i}}$ and ${\mathrm{S}}_{\mathrm{j}}$ are the areas of patches $\mathrm{i}$ and j, and ${\mathrm{L}}_{\mathrm{i}\mathrm{j}}$ is the cumulative resistance value of the corridor between patches $\mathrm{i}$ and $\mathrm{j}$. ${\mathrm{L}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is the maximum value among all corridor resistances.

Based on the important ecological corridors identified by the Gravity model, three indices were used to quantitatively characterize the ecological network characteristics of the study area: network closure index (α), line-point ratio (β), and network connectivity index (γ) [88]. The formulas are as follows:

$$\begin{array}{rr}& \mathsf{\alpha }={\displaystyle \frac{\mathrm{L}-\mathrm{n}+1}{2\mathrm{n}-5}}\\ & \mathsf{\beta }={\displaystyle \frac{\mathrm{L}}{\mathrm{n}}}\\ & \mathsf{\gamma }={\displaystyle \frac{\mathrm{L}}{3\left(\mathrm{n}-2\right)}}\\ & \end{array}$$

(9)

In the formula, $\mathrm{L}$ represents the number of corridors and nodes. $\mathsf{\alpha }$ higher index value signifies more closed loops within the network, suggesting smoother internal material and energy flow. $\mathsf{\beta }$ the larger the index value, the more corridors connected to each node within the network, indicating a more complex corridor network structure; $\mathsf{\gamma }$ the larger the index value, the closer the actual number of corridors connecting nodes is to the possible number of corridors, indicating higher network connectivity [89].

3. Results

3.1. Analysis of Land Use Structure and Transfer Changes

In the pisha sandstone area in the middle reaches of the Yellow River, the grassland area was 13,872.73 km², 14,449.10 km², and 14,185.06 km² in 2003, 2013, and 2023, accounting for 82.70%, 86.14%, and 84.57% of the total land area, respectively. The proportion of impervious Surface increased steadily from 1.41% to 3.02%. The total area of farmland decreased by 433.75 km², with its proportion declining from 13.80% to 11.22%. The proportion of unused land decreased from 1.86% to 0.82%. The forest land area increased by 22 km². The water body area first increased and then decreased, with an overall increase of 2.65 km² and its proportion rising from 0.22% to 0.24% (Table 3 and Table 4). Summarizing the above data, the conversion between various land use types was frequent and significantly different from 2003 to 2023. Grassland dominated the land use in the pisha sandstone area, while the areas of farmland and unused land decreased, and the areas of forest land and impervious surface continued to increase.

The chronological map of land use conversion in the pisha sandstone region (Figure 4),reveals significant changes from 2003 to 2013. During this period, farmland was predominantly converted out of the region, with a total converted-out area of 2,315.26 km² compared to a converted-in area of 427.02 km², leading to a net loss of 1,888.24 km². Specifically, 934.35 km² of farmland transitioned to grassland, 39.95 km² to impervious surface, 1.73 km² to water bodies, and 0.81 km² to forest.In the subsequent decade, from 2013 to 2023, grassland experienced the most significant outflow, totaling 868.68 km². Of this, 661.91 km² was converted to farmland, 99.46 km² to impervious surface, and 13.17 km² to forest land.The main land types converted to grassland, from highest to lowest, are farmland 526.93km2, unused land 75.90km2, impervious surface 6.09km2, water body 1.67km2, impervious surface 0.11km2, forest 0.02km2. The overall amount of grassland conversion outflow is greater than the amount of grassland conversion inflow, and the amount of grassland decreases.

In general, before 2013, farmland had the largest outflow area and showed a decreasing trend, which was directly related to the national “Grain for Green” policy launched in 2002. As a typical erosion-sensitive area, the pisha sandstone area in the middle reaches of the Yellow River became a key implementation area for this project, resulting in the large-scale conversion of farmland to forest land and grassland. After 2013, grassland had the largest outflow area, and impervious surface expanded significantly, with other land use types also undergoing corresponding adjustments. This was mainly driven by the advancement of urbanization and regional economic development.

3.2. Spatio-Temporal Analysis of GI Based on MSPA

From a temporal perspective, the GI area accounted for 84.66% to 87.70% of the study area between 2003 and 2023 (Figure 5). While the total GI area initially increased before declining, it maintained an overall upward trend. The proportion rose from 84.66% in 2003 to 89.21% in 2013, representing a 4.55% increase. From 2013 to 2023, it decreased slightly by approximately 1.5%, reaching 87.70% in 2023.From the perspective of GI types, the core area was the largest, accounting for approximately 40%-54% of the total GI area. Based on the ecological significance of the core area, it can be inferred that the study area has contiguous large-scale green spaces. Secondly, the areas of bridges, edges, perforations, and loops accounted for relatively large proportions of the total study area, at 20%, 7%, 10%, and 7%, respectively. The existence of bridges and other elements enhances the connectivity of the GI ecological network and ensures the integrity of the green spaces in the study area. In addition, the proportions of islets and branches were relatively small, at 0.40% and 1.50%, respectively.

From a spatial perspective (Figure 6), the core areas were mainly distributed in the northern, northwestern, and southern boundaries of the study area from 2003 to 2023. Among these, the core areas in Dongsheng, Dalate Banner, Hangjin Banner, Kangbashi, and Yijinhuoluo Banner were widely distributed. Types such as bridges, branches, and islets were mainly distributed in the southern, eastern, and southeastern parts of the study area, including Shenmu City, Fugu County, and the eastern part of Jungar Banner. The topography of the study area slopes from northwest to southeast. The southeastern part consists of high hills with steep slopes, which are less suitable for plant growth. The eastern part of the Pisha Sandstone Formation along the west bank of the Yellow River features low mountain terrain with severe soil erosion and a relatively small core area. It is prone to forming types such as connecting bridges, loops, and isolated islands.(Figure 6: C, D).The northwest region consists of low hills with higher elevations, featuring larger areas of green patches and a widespread core distribution.A distinct boundary line appears between the core and other types. As shown in Figure 2, this boundary line aligns closely with the contour lines. From the perspective of background areas, the background areas of urban centers such as Dongsheng, Shenmu City, Baode County, and Jungar Banner showed an expanding trend (Figure 6: A, B). All the above factors result in the wide distribution of core areas in the northwest.

Overall, the changes in GI types and their distribution are mainly attributed to two factors. On the one hand, the fragmented topography and gully development in the pisha sandstone area are not conducive to the integrity and connectivity of the ecosystem. On the other hand, frequent human activities and energy development, especially coal mining, have hindered vegetation restoration and intensified landscape fragmentation. Despite the decrease in the area of some GI types, the gradual expansion of the core area reflects its important role in the GI ecological network. This indicates that with the promotion of ecological protection and governance measures, the ecological space in the study area is developing in a positive direction despite the dual pressures of human activities and natural conditions, but there is still room for improvement in the systematicity of the entire ecological space.

3.3. GI Landscape Pattern Analysis Based on Landscape Indices

Using Fragstats software to calculate landscape pattern indices, we screened indicators based on their ecological significance. Six key metrics were selected: NP, LPI, LSI, CONTAG, SHDI, and SHEI to analyze landscape patterns in the study area from 2003 to 2023.

From 2003 to 2013, the NP increased from 10,242 to 11,259, showing a relatively stable growth trend. The LPI decreased from 48.5126 to 47.2493, indicating a downward trend. The LSI increased from 100.5386 to 112.8338, indicating an upward trend. The CONTAG slightly increased from 98.5602 to 98.8147. SHDI decreased from 0.0329 to 0.0282, and SHEI decreased from 0.0183 to 0.0158. The special geomorphology of the study area has led to unstable surface conditions. Small-scale land development and early mineral resource utilization by humans have intensified soil erosion to a certain extent, resulting in a certain degree of landscape fragmentation. The core patch area has been reduced under the combined influence of human activities and soil erosion, and the complexity of the landscape shape has increased. However, the overall continuity of the landscape remains high, and the diversity and evenness have slightly decreased. This reflects that the interference between human activities, natural geomorphology, and soil erosion factors is in the initial stage and has not yet caused drastic impacts on the landscape pattern.

From 2013 to 2023, the NP increased sharply from 11,259 to 21,811. The LPI continued to decline to 45.4556. The LSI further increased to 119.7655. The CONTAG significantly decreased to 93.5171. SHDI significantly increased to 0.2062, and SHEI also increased to 0.1059. During this period, the intensity of human activities increased significantly. Large-scale mineral resource development led to drastic changes in land use, further damaging the protective layer of surface vegetation and aggravating soil erosion. Under the dual effects of human activities and soil erosion, the contiguous landscape was severely fragmented, the degree of fragmentation increased sharply, the core patch area was further eroded, and the landscape shape became increasingly complex.

Figure 7. Area and proportion of landscape types in 2003, 2013, and 2023.

Figure 7. Area and proportion of landscape types in 2003, 2013, and 2023.

3.4. GI Network Construction

3.4.1. Ecological Resistance Surface Analysis

A comprehensive resistance surface system for the study area was constructed by integrating resistance factors such as land use type, elevation, slope, distance to roads, and NDVI, resulting in resistance surfaces for 2003, 2013, and 2023 (Figure 8).Over the past 20 years, the high-resistance area has shown an expanding trend, with an overall spatial distribution characterized by higher resistance in the northwest and lower resistance in the southeast. Although the northwestern part of the pisha sandstone area is dominated by low hills with high elevations and relatively gentle terrain, its low vegetation coverage directly increases the resistance to ecological flows, leading to overall high resistance values in the northwest.

Meanwhile, urban centers and energy bases such as Dongsheng, Shenmu City, Baode County, and Jungar Banner characterized by high urbanization levels, dominated by impervious surface and farmland, intensive human activities, dense road networks, and large-scale mineral extraction and urban expansion that encroach on ecological space,have restricted biological migration and ecological process diffusion, becoming concentrated areas of high resistance values. Low-resistance areas are distributed across all counties and cities, with widespread distribution in the southeastern region. Despite the lower elevation in the southeast, the ecological background of steep slopes and discontinuous vegetation still restricts local resistance values, making it difficult to form low-resistance connectivity zones. The spatiotemporal evolution and pattern characteristics of GI analyzed in this study can effectively reduce soil erosion, avoid environmental damage caused by human activities, protect forest resources, and provide a theoretical basis for the coordination of economic development and environmental protection.

3.4.2. Simulation of Potential GI Corridors

Based on the resistance surface of the study area, the minimum cost distance method was used to calculate the minimum cumulative cost paths between GI source areas and construct potential species migration corridors. As shown in Table 5 and Figure 9, there were 190 GI corridors in the study area in 2003, densely distributed in the central and northwestern regions where the terrain is relatively flat and patch connectivity is high. In 2013, the number of GI corridors decreased by 36.8% to 120 compared with 2003, but the phenomenon of corridor intersection slightly increased. This is mainly due to the merging of some ecological source areas into patches and the increase in ecological source areas in the northeast and south, leading to more corridor intersections. This phenomenon indicates that despite the decrease in the number of corridors for species exchange, local ecological connectivity has become more complex. In 2023, the number of GI corridors further decreased by 24.2% to 91 compared with 2013. The area of ecological green spaces in the study area increased, and the distribution of corridors tended to be simplified.

Over the past 20 years, the number of GI corridors in the study area has shown a decreasing trend, which is mainly attributed to the fragile characteristics of the regional ecosystem and human interference: 1) The pisha sandstone area is characterized by severe soil erosion, leading to the fragmentation and damage of ecological corridors and a continuous decrease in their area. The loose soil texture is susceptible to wind and water erosion, making it difficult for vegetation to grow stably, further destroying the connectivity of ecological corridors and exacerbating soil erosion in the pisha sandstone area. 2) With the socio-economic development of cities, the continuous expansion of human activities and impervious surface has continuously encroached on GI source areas in the study area, reducing the connection starting points of corridors and leading to a decrease in the number of GI corridors. 3) The GI source areas in the southeast have fragmented into small patches, leading to a temporary increase in corridor intersections.However, this connection does not indicate improved connectivity but rather a sign of weakened corridor functions. Eventually, the continuous reduction of source areas leads to a decrease in corridor intersections and an overall decline in connectivity.

3.4.3. Construction of GI Networks Based on Gravity Model Analysis

The Gravity model was used to calculate the interaction force between GI source areas, construct a gravity matrix of GI source areas, and distinguish the importance of corridors. Based on the results of the gravity matrix analysis, GI corridors with an interaction force ≥20,000 were classified as important corridors, and those with an interaction force <20,000 were classified as general corridors. In 2003, there were 190 GI corridors in the study area, including 38 important corridors and 152 general corridors. In 2013, the total number of GI corridors in the study area was 120, including 21 important corridors and 99 general corridors. In 2023, the number of GI corridors in the study area continued to decrease to 91, including 16 important corridors and 75 general corridors(Figure 10).

Important corridors are usually located in areas with stable ecological environments, gentle terrain, and good connectivity, with weak impacts of water and wind erosion, relatively slight soil erosion, less external interference, favorable connectivity backgrounds of ecological sources, and low resistance values. They are mainly distributed in the northwest of the pisha sandstone area. In contrast, general GI corridors are numerous and widely distributed, limited by natural conditions and disturbed by human activities, making them unsuitable for large-scale species exchange activities. Based on the identified GI source areas and important corridors, a potential GI network was constructed.

3.4.4. Evaluation of GI Network Structure

From 2003 to 2023, the number of ecological source sites decreased from 20 to 14, and the number of important ecological corridors decreased from 38 to 16 (Table 6). The α value decreased from 0.54 to 0.13, indicating a reduction in closed loops within the network. The β value remained above 1, suggesting a mesh-like structure of the regional network. However, the connectivity level of this ecological network is generally moderate, with limited resilience to external pressures. The γ value decreased from 0.70 in 2003 to 0.44 in 2023, indicating that there is room for improvement in network connectivity. Based on this series of network structure indices and comparative analysis with relevant literature [90,91,92], it is evident that the overall spatial structure of the regional network is relatively underdeveloped. Its layout is biased toward the northwest, presenting a fragmented pattern characterized by local concentration and overall sparseness.

The evaluation of network connectivity shows that the overall connection tightness and connectivity of the GI network have decreased. The evolution of the GI network in the study area from 2003 to 2023 indicates that the main driver of changes in the GI network pattern is the change in GI source areas. Therefore, changes in the number and area of GI source areas have important impacts on network nodes and corridors, ultimately altering the original network structure.

4. Discussion

4.1. Necessity of GI Identification and Analysis in the Middle Reaches of the Yellow River

As a key ecologically fragile area in the Yellow River Basin, the pisha sandstone area in the middle reaches of the Yellow River requires the identification and systematic analysis of GI, which is an essential requirement for addressing core regional ecological challenges and overcoming the limitations of traditional governance models. Its necessity is reflected in three aspects: 1) From the perspective of regional ecological characteristics, this area is the largest concentrated pisha sandstone area in China with the highest erosion intensity. Its high porosity and weak cementation cause it to disintegrate into mud when wetted and disperse as sand when exposed to wind, with an erosion modulus of 30,000–40,000 t/(km²·a). Although accounting for only 2% of the Yellow River Basin area, it contributes 1/4 of the coarse sediment deposition in the lower reaches. This area is a key challenge for soil erosion control in the basin, with extremely poor inherent ecosystem stability, urgently requiring GI as a systematic carrier to establish targeted protection systems. 2) Analysis of the limitations of traditional governance models shows that past research has mainly focused on micro-scale erosion processes and the soil and water conservation effects of individual vegetation types. Governance has primarily adopted fragmented approaches such as point, line, and surface-based measures, making it difficult to achieve overall ecosystem optimization. Such fragmented approaches fail to connect ecological patches fragmented by gully erosion and mining activities or restore soil and water conservation corridors interrupted by urban expansion, thus being ineffective in curbing landscape fragmentation and ecological degradation. The integrated “patch-corridor-network” structure of GI precisely addresses this deficiency. 3) From a practical perspective, this region faces dual pressures of ecological protection and economic development. Urbanization and energy extraction have led to frequent land use conversions, increasing ecological uncertainty. Through GI identification and analysis, priority ecological spaces and management focuses can be clarified, providing a scientific basis for territorial spatial planning and preventing the escalation of ecological and economic conflicts. In summary, conducting GI research in this area is not only an inevitable choice for addressing unique ecological challenges but also a key support for the ecological protection and high-quality development strategy of the Yellow River Basin.

4.2. Driving Mechanism and Ecological Significance of GI Spatiotemporal Evolution in the Pisha Sandstone Area

The spatiotemporal evolution of GI in the pisha sandstone area in the middle reaches of the Yellow River is the result of the combined effects of natural constraints and human activities, and its change process is directly related to the stability of regional ecological functions. From the perspective of driving mechanisms, the proportion of the total GI area increased from 84.66% to 89.21% between 2003 and 2013, with a simultaneous expansion of core areas. This trend is directly attributable to the “Grain for Green” policy implemented in 2002, which promoted the large-scale conversion of farmland to grassland and forest land, directly improving vegetation coverage and laying the foundation for the initial establishment of the GI network. From 2013 to 2023, the GI area decreased to 87.70%, with grassland becoming the primary outflow land use type. This trend is closely linked to urban expansion and energy development, as human activities encroach on ecological space, leading to the fragmentation of GI patches. Furthermore, natural terrain imposes rigid constraints: the gently sloping low hills in the northwest with slight erosion form the core concentration area of GI, while the steep high hills in the southeast with severe erosion have fragmented GI patches, resulting in a pattern of “concentration in the northwest and sparseness in the southeast”.Ecologically, these structural changes in GI directly affect soil and water conservation and ecological connectivity. From 2003 to 2023, the expansion of core areas enhanced resistance to soil erosion, potentially reducing the amount of coarse sediment entering the Yellow River. However, landscape fragmentation intensified, with ecological source sites decreasing from 20 to 14 and important corridors shrinking from 38 to 16. The connectivity of the GI network declined, weakening the capabilities of species migration and material cycling and increasing the risk of soil erosion during extreme precipitation events. Meanwhile, the spatial differentiation of GI indicates that the southeastern region with severe erosion coincides with GI vulnerable zones, which require prioritized ecological patch supplementation and corridor restoration, providing a scientific basis for determining restoration priorities. The findings of this study are largely consistent with existing research conclusions. For example, Liu et al. [93] found that the ecological environment quality in the pisha sandstone area has generally improved over the past two decades but is highly sensitive to human activities. Yao et al. [94] observed that although the ecological environment quality in the pisha sandstone area was generally poor between 2000 and 2010, it showed significant improvement by 2010. These findings all indicate that the ecological environment in the pisha sandstone area faces severe challenges.

4.3. Limitations and Future Research Directions

Compared with traditional GI identification and analysis methods, the combination of MSPA and the MCR model enables more efficient and accurate core area selection, source site extraction, and corridor screening. Incorporating a resistance model to evaluate corridor importance enhances the scientific rigor and rationality of ecological network construction. The edge width setting in MSPA analysis affects the total area of core zones. Multiple tests were conducted by adjusting the edge width threshold, with the raster pixel size set to 30 m × 30 m (edge width value of 3), ensuring sufficient resolution while adequately representing landscape types.The identification of ecological source areas based on MSPA and landscape connectivity methods shows that the distance threshold setting in Conefor 2.6 software affects the IIC and PC values. In this study, the threshold was set to 20,000 with a connectivity probability of 0.5. By comparing the effects of different IIC and PC values on the total area and number of ecological source patches, use equal weights method, and values greater than 1 were ultimately selected as the screening threshold for ecological sources, enabling landscape connectivity analysis and ecological function evaluation.Although this study determined distance thresholds, connectivity probabilities, and ecological source screening thresholds through comprehensive comparative analysis, the pisha sandstone area exhibits diverse topography and severe soil erosion. Therefore, future research should integrate the spatial variation characteristics of topographic differentiation and soil erosion in the pisha sandstone area to develop more precise threshold determination methods. Regarding the construction of ecological resistance surfaces, there is no unified standard in domestic and international ecological network research. Most current studies assign baseline resistance values to landscapes based on expert scoring methods. This research selected resistance factors including land use, elevation, slope, distance to roads, and vegetation coverage, with the selection of resistance factors and weight settings relying on previous studies, which may introduce a certain degree of subjectivity. Future research should adopt more scientific and standardized quantitative methods.For the construction of ecological corridors, this study used the MCR model to extract ecological corridors and the Gravity model to identify important ecological corridors. Previous research on the Gravity model has shown that threshold selection tends to be somewhat subjective. While the theoretically derived data are scientific and objective, practical construction processes must consider various factors, including human-induced elements. Further research is needed to determine how to establish distance thresholds that effectively integrate real-world conditions. Future research may propose optimized GI schemes tailored to specific sub-regions based on the topographic variations and soil erosion patterns in the pisha sandstone area of the middle reaches of the Yellow River.

5. Conclusions

Based on land use data from 2003, 2013, and 2023, this study integrated MSPA, landscape index method, and MCR model to identify and analyze GI in the pisha sandstone area in the middle reaches of the Yellow River in three periods. The main conclusions are as follows:

(1)

From 2003 to 2023, grassland remained the dominant land use type in the study area, accounting for 82.70% to 84.57% of the total land area. Impervious surface increased by 270.48 km², while farmland and bare land decreased by 433.35 km² and 174.11 km², respectively. Forest land and water bodies expanded by 22 km² and 2.65 km², respectively. Before 2013, farmland was the primary outflow land use type, with bare land decreasing and forest land, water bodies, impervious surface, and grassland increasing,driven by the “Grain for Green” policy. After 2013, grassland became the primary outflow land use type, with water bodies decreasing and farmland, forest land, impervious surface, and bare land increasing,characterized by significant impervious surface expansion influenced by urbanization and economic development.

(2)

From 2003 to 2023, the total area of GI in the study area showed a trend of initial increase followed by decrease, maintaining a proportion between 84.66% and 87.70%. Spatially, it exhibited a pattern of concentration in the northwest and sparseness in the southeast. As the primary component of GI, core areas showed an annual increase in area, accounting for 47.07% to 54.42% of the total GI area, with a spatial distribution concentrated in the northwest of the study area. Other GI types, such as bridges, edges, and loops, accounted for smaller proportions and were scattered between core areas. Despite certain achievements in ecological restoration, the GI landscape has become increasingly fragmented and heterogeneous due to inherent low erosion resistance, development activities, and energy extraction.

(3)

From 2003 to 2023, the study area exhibited a spatial distribution pattern of higher resistance values in the northwest and lower values in the southeast, with high-resistance areas continuously expanding. The number of ecological source sites decreased from 20 to 14, the number of general ecological corridors reduced from 152 to 75, and the number of important ecological corridors declined from 38 to 16. This reduction in the number of corridors is directly related to soil erosion fragmentation and source area fragmentation.

(4)

From 2003 to 2023, the network closure index (α index) decreased from 0.54 to 0.13, the line-point ratio (β index) dropped from 1.90 to 1.14, and the network connectivity index (γ index) fell from 0.70 to 0.44. This indicates that the overall spatial structure of the regional GI network is relatively underdeveloped, presenting a fragmented pattern characterized by local concentration and overall sparseness. Changes in the number and area of GI source areas are the primary drivers of network pattern evolution.