Landscape Ecological Risk Assessment of Saihanba under the Change in Forest Landscape Pattern

Abstract

Examining the Saihanba Mechanical Forest Farm, this study utilized Landsat remote sensing data from 1987, 1997, 2001, 2013, and 2020 to interpret land use from the Support Vector Machine (SVM) method, and to decipher evolving land use patterns over the last four decades. Grounded in landscape ecology theory, an innovative evaluation index for landscape ecological risk was introduced, leading to the delineation of 382 ecological risk evaluation units. Employing landscape pattern indices and a method of spatial autocorrelation, we analyzed the spatial and temporal distribution characteristics and spatial correlation patterns of landscape ecological risk across five distinct periods. Geostatistical approaches were used to explore the driving factors of landscape risk. The results indicate that since 1987, there have been significant changes in land use types, especially in forest landscapes, their proportion increasing from 23.19% to 74.55%. In 1987, the proportion of high-risk areas was 72.30%, but in 2020, high-risk areas had significantly decreased and clustered in specific locations. The landscape ecological risks in each period of the study area showed a positive spatial correlation and tended to gather in space. After comprehensive exploration using a geographic detector, we found that landscape type, temperature, and vegetation coverage are the main risk factors. Among them, landscape type has the greatest impact on the landscape and works together with slope, aspect, and precipitation. In forest farm management, only the adaptation and adjustment of single factors are often paid attention to, while the compound effects of multiple factors are ignored. The results of this study bring important reference value to the operation and development of forest farms.

1. Introduction

Land use/cover change (LUCC) is central to human development and utilization of the natural environment [1]. Landscape pattern change affects comprehensive geographical factors such as climate [2,3,4], soil [5,6] and water [7,8,9]. Additionally, the value of regional ecosystem services are affected [10,11] which threatens the regional ecological environment health and increases the regional risk [12,13]. Since the establishment of Saihanba Mechanical Forest Farm in 1962, the landscape pattern has undergone major changes. In order to protect the ecological environment and sustainable development, it is necessary to carry out a risk assessment of its ecological environment. Since its inception in 1962, Saihanba Mechanical Forest Farm has undergone major alterations in its landscape pattern. To safeguard the ecological environment and promote sustainable development, it is imperative to conduct a comprehensive risk assessment of its ecological milieu.

At present, landscape risk assessments locally and abroad mainly focus on rivers [14,15,16,17] and coasts [17,18,19,20]. Many scholars have conducted forest landscape risk assessments on climate change [21], strong winds [22], regulated floodplains [23], population [24], and fire [25]. Tanja et al. [26] investigated the ecological risks of the Bohai Sea landscape in China by analyzing pesticide residues in freshwater systems. Hossain et al. [27] analyzed risk and resilience in Swiss Alpine communities from the perspective of people and environment. Yanjie et al. [19] evaluated the community-level risks induced by marine ecosystems by fitting different models. Rasoul and R. [21] proposed that forest management efforts should be made to reduce the ecological risk of forest landscapes. Hua et al. [25] carried out a time series analysis on wildfire and forest landscape risk from the large-scale spread of the mountain pine beetle (MPB) epidemic, which was found to have a multitude of explanatory variables. At present, research on the ecological risk of forest landscapes is limited. We assess the landscape ecological risk drawing from the landscape pattern of Saihanba Mechanical Forest Farm, and then consider the influence of factors on the comprehensive ecological risk assessment results.

Currently, the research around the ecological risks associated with land uses two types of evaluation models. The first model uses a “source correlation–receptor evaluation–exposure and hazard evaluation–risk characterization” model. For example, scholars like Yanes [17] and Walker [15] generated an index system for ecological risk assessment from three aspects: risk source intensity, receptor exposure, and risk effect [28,29]. The second model involves direct evaluation of the landscape pattern to evaluate the landscape ecological risk. Scholars such as Ayre [16] and Dale VH [14] adopt a landscape ecology perspective, using a landscape ecological index to depict the ecological effect of LUCC change [18,30,31].

Using a geographic detector method to identify risk drivers can better ensure the accuracy of landscape risk assessment and provide scientifically sound suggestions for forest farms. The geodetector, grounded in spatial variance analysis theory, examines the correlation between factor variables and result variables. It assesses the strength of influence of each factor variable on the result variable, discerns the disparities in influence among different factor variables, and determines whether the impact of each factor variable on the result variable is independent or interactive [32]. Widely applied in diverse domains such as urbanization efficiency [33,34], population aging [35], and medicine [36], this method has recently found application in ecological research, examining insect diversity [37], soil heavy metal detection [38,39], and the intersection of air quality and social economy [40,41,42,43]. For example, Liu et al. [37] quantified the interactive effects of 15 variables on beetle distribution using the geographical detector method, successfully analyzing the associated risk drivers. Xu et al. [43] used a geographical detector to quantify the driving factors affecting air quality and establish national key functional areas. In Mingrui Li et al.’s research on the Ertix River Basin in Central Asia, they applied geographic detectors and geographically weighted regression to detect landscape risk [44]. However, there is a dearth of research on the detection of ecological risk factors based on remote sensing technology.

This study focuses on the Saihanba Mechanical Forest Farm as the research area. Since the establishment of the forest farm, the vegetation cover grade has improved, leading to further restoration of the ecological environment and consequential significant changes in the landscape pattern. To capture the diverse developmental stages of the study area, Landsat remote sensing images from 1987, 1997, 2001, 2013, and 2020 were selected. The support vector machine method was applied to generate the landscape classification map, and field survey data were utilized for validation. The NDVI was employed to obtain the vegetation index, while the landscape ecological risk index was utilized to investigate the spatial and temporal distribution characteristics of landscape ecological risk. The spatial autocorrelation method was used to assess landscape spatial correlation characteristics, and the geographical detector method was employed to explore the driving factors of landscape ecological risk. The findings derived from this research can furnish a scientific foundation for land use planning and forest management rooted in ecological security considerations.

2. Materials and Methods

2.1. Study Area

The Saihanba Mechanical Forest Farm is located in the semi-arid and semi-humid zone (Figure 1), which is positioned at the intersection of the Yinshan Mountains, the Daxing’an Mountains, and the Hunshandake Sandy Land. It falls within the transitional area both below and above the dam, marking the juncture of forest–grassland and arid–semi-arid regions. The topography is characterized by a north-to-south slope, creating a distinctive landform with higher elevations in the north and lower elevations in the south. The altitude ranges from 1010 to 1940 m, with an average annual temperature of −1.2 °C, an average annual sunshine duration of 2548.7 h, and an average annual precipitation of approximately 452.2 mm, primarily concentrated in the months of June to September.

The Saihanba Mechanical Forest Farm has strategically embraced the concept of the “four banks” of forests as a pivotal strategy for advancing high-quality development in forestry and grassland management. This approach serves as a valuable model for fostering integrated protection and systematic governance across diverse ecosystems, including mountains, rivers, forests, farmland, lakes, and grasslands. Leveraging the synergistic benefits of the forest “four reservoirs” concept, the farm aims to maximize high-level advantages, thereby promoting development and enhancing the well-being of the local population. Presently, the predominant forest resources consist of artificially cultivated pure forests, with key tree species including Larix principis-rupprechtii, Pinus sylvestris var. mongolica, Picea asperata, Betula platyphylla, among others.

In history, the Saihanba area was a hunting paddock of the Qing government, with a good natural environment and many trees. Due to a lack of finance, a large number of trees were cut down and the forest coverage rate decreased sharply. Subsequently, the Japanese invasion of China caused a serious fire, and the destruction of forest resources became more serious. After the founding of the People’s Republic of China, it began to pay attention to ecological construction.

Currently, the total value of Saihanba Forest Farm forest resources has reached CNY 20.2 billion, the total value of forest resources is about USD 2.8 billion, driving the local social income of more than USD 80 million every year. Saihanba will build a global ecological civilization and green development demonstration zone in the next, striving to create higher social and economic benefits.

The establishment of Saihanba Mechanical Forest Farm has alleviated the ecological crisis in the Beijing–Tianjin–Hebei region and even the entire north region of China. But while the external crisis has been alleviated, problems such as single tree species and fragmentation of patches have appeared inside the forest farm, and the landscape risk is unknown.

2.2. Data Source

The Landsat remote sensing images from 1987, 1997, 2001, 2013, and 2020 used in this paper were acquired from the United States Geological Survey (https://glovis.usgs.gov/). Analysis of the images was conducted using ENVI 5.1 software. The support vector machine classification method was applied to categorized the preprocessed images. To align with the research objectives, the images were segmented into five distinct classes: forest, grassland, wetland, sandy land, and construction land (Figure 2). To ensure the accuracy of the classification, a confusion matrix was established, and verification results indicated that the kappa coefficients for the interpreted land use maps during the five periods consistently exceeded 0.76, meeting the accuracy standards for medium-resolution remote sensing images.

Precipitation, soil type, DEM, vegetation cover grade, air temperature, slope aspect, slope, and subcompartment accumulation were selected from the second-class survey data of subcompartment in Saihanba Mechanical Forest Farm.

2.3. Methods

2.3.1. Analysis of Temporal and Spatial Changes in Land Use and Vegetation Cover

The SVM classification method revealed significant landscape type changes in 1987–2020. Therefore, the transfer matrix method was used to explore the amount of exchanges of each landscape type over the entire period of 1987–2020, and the NDVI was calculated to see if there was a trend of decreasing vegetation cover classes.

Normalized Vegetation Index

The normalized difference vegetation index ($\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}$) is extensively used in the research and application of vegetation cover changes, and has a linear relationship with the vegetation distribution density [45,46,47]. Its formula is:

$$\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}=(\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}\mathrm{E}\mathrm{D})/(\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}\mathrm{E}\mathrm{D})$$

The notation used in this context defines NIR as the reflection value in the near-infrared band, and RED as the reflection value in the red band. Specifically, in Landsat 2 TM images, these correspond to band 6 for NIR and band 5 for RED. In Landsat 5–8 TM images, the relevant bands are band 4 for NIR and band 3 for RED.

Transfer Matrix Method

The transition matrix model serves to elucidate the initial probability of various stages and the transitions between these stages, thereby allowing the determination of change trends before and after different periods. The formula for this model is:

$$S_{ij}=\left[\begin{array}{ccc}{S}_{11}& S_{12}& \cdots S_{1n}\\ S_{21}& S_{22}& \cdots S_{2n}\\ \begin{array}{c}⋮\\ S_{n1}\end{array}& \begin{array}{c}\cdots \\ S_{n2}\end{array}& \begin{array}{cc}\cdots & ⋮\\ \cdots & S_{nn}\end{array}\end{array}\right]$$

where $S$ reflects the coverage area for each vegetation grade; $n$ reflects the number of vegetation cover grade types; $i$ and $j$ reflect the vegetation cover grade types at the initiation and commencement of the study period, respectively.

2.3.2. Landscape Ecological Risk Analysis

Division of Landscape Ecological Risk Units

In this study, according to the national grid GIS standard “Geographical Grid” (GB12409-2009) and related research, the grid size was determined based on 2–5 times the average patch area in the study area [48]. The ArcGIS fishnet creation function was employed to partition the study area into ecological risk units of 2 km × 2 km, and these fishnets were subsequently segmented to generate a total of 382 ecological risk units. Utilizing landscape ecological risk value assignment, the semi-variation function within the Kriging interpolation method was applied for spherical fitting, yielding the landscape ecological risk level map across different periods.

Construction of Landscape Ecological Risk Index

The landscape index, through the comprehensive analysis of multiple indices, enables a quantitative reflection of the landscape pattern and its dynamic changes. In this study, landscape fragmentation, dispersion, dominance, vulnerability, and disturbance were designated as landscape ecological risk indices.

Calculating the landscape disturbance index through the landscape fragmentation index, landscape separation index, and landscape dominance index can reflect the complexity of spatial structure. The landscape vulnerability index reflects the ecological sensitivity of the study area. The higher the vulnerability, the worse the resistance to interference. The landscape loss index is calculated by combining the landscape disturbance index and the landscape vulnerability index. According to the results, the landscape ecological risk index was then calculated (

Table 1. Landscape ecological risk index calculation method.

Index	Formula	Meaning
Landscape Fragmentation Index $C_i$	$C_i=\frac{n_i}{A_i}$	$A_i$ is the total area of the landscape type $i$; A is the total area of the landscape; $n_i$ is the number of patches of the landscape type $i$; $Q_i$ = the number of squares in which the patches $\mathrm{i}$ appear/the total number of squares; $M_i$ = the number of patches $i$/the total number of patches and; $L_i$ = the area of the patches $i$/the total area of the quadrat.
Landscape Separation Index $N_i$	$N_i=\frac{A}{{2A}_i}\sqrt{\frac{n_i}{A}}$
Landscape Dominance Index $D_i$	$D_i=\frac{Q_i+M_i}{4}$+$\frac{L_i}{2}$
Landscape Disturbance Index $E_i$	$E_i=$a$C_i+$b$N_i+$c$D_i$	a, b, and c represent the weights of fragmentation, separation, and dominance indices. Referring to relevant studies [49,50], the assigned values are determined to be 0.5, 0.3, and 0.2.
Landscape Vulnerability Index $F_i$	Expert consultation method and normalized treatment	Referring to relevant research [51], the vulnerability assignments of various landscapes in the study area were determined as follows: 5 for sandy land, 4 for wetland, 3 for grassland, 2 for forest, and 1 for construction land. Forest secondary classification landscape assignment: 5 for larch, 4 for sycamore pine and birch, 3 for spruce, 2 for mixed birch forest, and 1 for mixed forest of fallen clouds.
Landscape Loss Index $R_i$	$R_i=E_i{\times F}_i$	$A_{ki}$ is the area of type i landscape in landscape ecological risk assessment unit $k$; $A_k$ is the total area of landscape ecological risk assessment unit $k$; ${ERi}_i$ is the landscape ecological risk index of landscape ecological risk assessment unit k, and its value is positively correlated with the degree of ecological risk.
Landscape Ecological Risk Index ${ERi}_i$	${ERi}_i={\displaystyle \sum _{i=1}^N}\frac{A_{ki}}{A_k}\times R_i$

).

2.3.3. Spatial Autocorrelation Analysis

Spatial autocorrelation analysis was employed to unveil the spatial correlation characteristics of a specific attribute unit and its neighboring attribute units based on eigenvalues [52]. It was used to measure the distribution characteristics and interrelationships of spatial risk data.

Global Autocorrelation

The global Moran’s I index was utilized to assess the spatial correlation of attribute values within a unit across the entire study area.

$$\mathrm{G}\mathrm{l}\mathrm{o}\mathrm{b}\mathrm{a}\mathrm{l} \mathrm{M}\mathrm{o}\mathrm{r}\mathrm{a}{\mathrm{n}}^{\mathrm{’}}\mathrm{s} I=\frac{{\sum }_{i=1}^n{\sum }_{i=1}^m{W}_{ij}\left(x_i-\overline{x}\right)\left(x_j-\overline{x}\right)}{S^{2{\sum }_{i=1}^n{\sum }_{j=1}^m{W}_{ij}}}$$

$$S^2=\frac{1}{n}\sum _{i=1}^n{\left(x_i-x\right)}^2$$

$$\overline{x}=\frac{1}{n}\sum _{i=1}^n{x}_i$$

where $x_i$ represents the observation value of the $i$-th area, $n$ is the number of grids, and $W_{ij}$ is a binary adjacency space weight matrix, which is used to represent the adjacency relationship of spatial objects. When the area $i$ and the area $j$ are adjacent, $W_{ij}$ = 1; when the area $i$ and the area $j$ are not adjacent, $W_{ij}$ = 0.

Local Autocorrelation

The local Moran’s I index was employed to reflect the correlation between the attribute value of a unit and the adjacent spatial units.

$$\mathrm{L}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l} \mathrm{M}\mathrm{o}\mathrm{r}\mathrm{a}{\mathrm{n}}^{\mathrm{’}}\mathrm{s}\mathrm{ }{I}_i=\left(\frac{x_i-\overline{x}}{m}\right){\sum }_{j=1}^n{W}_{ij}(x_i-\overline{x})$$

$$m=\left({\sum }_{j=1j\ne i}^n{x}_j^2\right)/\left(n-1\right)-x^2$$

The Moran’s I value ranges from −1 to 1. A positive Moran’s I (>0) suggests a positive correlation in the study area, indicating that the attribute values of the study units are convergent. Conversely, a negative Moran’s I (<0) signifies a negative correlation and a dispersed distribution of attribute values. When Moran’s I equals 0, it indicates the absence of spatial correlation.

2.3.4. Geographical Detector

The geographical detector model, optimized for parameters, encompasses four functions: factor detector, interaction detector, ecological detector, and risk detector, each capable of handling various types of variables [53,54]. It can be used to test the stratified heterogeneity of a single factor; it can also be used to detect a possible causal relationship between two factors by testing the coupling of their spatial distributions.

In this study, the factor detector, interaction detector, and ecological detector were employed to investigate the comprehensive influence of diverse risk factors in the forest landscape on the risk level of the Saihanba forest landscape in 2020.

Factor Detector

Factor detector was used to measure the impact strength of each risk factor on the risk of forest landscapes. The detection method is:

$$q=1-\frac{1}{N{\sigma }^2}\sum _{h=1}^L{N}_h{\sigma }_h^2$$

where $q$ is the influence of a certain factor on the spatial distribution of landscape ecological risk, the value range is [−1, 1]; L is the partition number; $h$ is composed of $N_h$ units; $N\sigma$ and $N$ are the unit numbers of layer $h$ and the whole region, respectively.

Ecological Detector

The ecological detector reflects whether the risk factors have significant differences in the spatial distribution of landscape ecological risk, and it is measured by the F test. The formula is as follows:

$$F=\frac{{N_{n=1}(N}_{n=2}-1){ \sigma }_{n=1}^2}{{N_{n=2}(N}_{n=1}-1){\sigma }_{n=2}^2}$$

where $N_{n=1}$, $N_{n=2}$ represents the sample size in the two risk factor partitions, and the null hypothesis ${\mathsf{\sigma }}_{n=1}^2$= ${\sigma }_{n=2}^2$. Rejecting the null hypothesis at the chosen significance level signifies the presence of significant differences between the two factors influencing the spatial distribution of landscape ecological risk.

Interaction Detector

In assessing whether an interaction exists among various risk factors impacting the spatial distribution of landscape ecological risk, we ascertained whether the two factors operate independently or interactively influence the spatial distribution of landscape ecological risk. This determination is achieved through a comparative analysis of the interactions between the two risk factors, encompassing the following five relationships:

If q (A∩B) < Min (q (A), q (B)), it means that the nonlinearity of the two factors is weakened;

If Min (q (A), q (B)) <q (A∩B) < Max (q (A), q (B)), it means that the single-factor nonlinearity is weakened;

If q (A∩B) > Max (q (A), q (B)), it means two-factor enhancement;

If q (A∩B) = q (A) + q (B), it means that the two factors are independent of each other;

If q (A∩B) > q (A) + q (B), it means two-factor nonlinear enhancement.

3. Research Results and Analysis

3.1. Analysis of Landscape Dynamic Change

3.1.1. Analysis of Structural Changes of Landscape Types

Figure 3 illustrates the evolution of landscape types in Saihanba from 1987 to 2020. Over this period, the forest and wetland areas had significant growth, while grassland, sandy land, and construction land areas declined. Particularly noteworthy is the significant expansion of woodland, registering a remarkable increase of 479.06 km², constituting approximately half of the total study area and elevating its proportion from 23.19% to 74.55%. The wetland area experienced a significant increase of 37.39 km², more than threefold compared to 1987, with its proportion in the total study area rising from 1.27% to 5.28%. In contrast, the areas of grassland, sandy land, and construction land decreased by 313.33 km², 192.37 km², and 10.75 km², respectively. Consequently, their proportional representation decreased from 52.13%, 22.06%, 1.36% to 18.53%, 1.43%, and 0.21%, respectively.

3.1.2. Landscape Type Structure Transfer Analysis

Changes in the landscape type transition matrix in the overall period from 1987 to 2020 are shown in

Table 2. Table of transfer area by landscape type, 1987–2020.

Matrix Transfer Area	1987
	Unit	Landscape Type	Forest	Grassland	Wetland	Sandy	Construction Land
2020	km2	Forest	169.73	397.42	5.32	119.68	5.76
	km2	Grassland	31.01	64.38	1.55	68.42	5.84
	km2	wetland	13.14	16.66	4.91	14.65	0.29
	km2	Sandy	2.12	7.43	0.07	3.05	0.41
	km2	Construction land	0.14	0.22	0.01	1.16	0.39

. In general, forests and wetlands were the transfer-in types, and grassland, sandy land, and construction land were the transfer-out types. A large amount of grassland, sandy land, and construction land were transferred to forests, while a small amount of grassland and sandy land were transferred to wetland. In addition, there were different degrees of transformation among the various landscape types.

Vegetation cover changes over the entire period 1987–2020 are shown in

Table 3. Table of transfer area by vegetation cover, 1987–2020.

Matrix Transfer Area			1987
	Unit	Vegetation Cover Grade	Very Low	Low	Middle	High	Very High
2020	km2	Very low	0.0117	0.0378	0.0549	0.2025	0.9927
	km2	Low	0.063	0.1098	0.2484	0.7776	10.5687
	km2	Middle	0.3798	0.7875	1.9584	5.3874	44.6472
	km2	High	2.4588	6.273	14.9481	32.1633	120.6054
	km2	Very high	19.5201	47.5461	85.0806	155.8125	379.8135

. The extremely high vegetation cover grade was the type of transfer-in, and the area was greatly increased. The vegetation cover grade was transferred out to different degrees. Among them, the extremely low, and low and medium coverage was greatly reduced, and the overall vegetation cover grade significantly improved.

3.2. Analysis of Temporal and Spatial Changes in Landscape Risk

To investigate the spatial and temporal distribution characteristics of landscape ecological risk in Saihanba, the ecological risk index (ERI) for each landscape risk evaluation unit served as its attribute value. The spatial and temporal distribution of landscape ecological risk in Saihanba was determined using the ordinary Kriging interpolation tool in ArcGIS. Considering the specific conditions of the study area and the ERI across five phases, the natural breakpoint method was applied to categorize the study area into five landscape ecological risk levels: high-risk area (ERI > 0.29), mid-high-risk area (0.29 < ERI ≤ 0.26), medium-risk area (0.26 < ERI ≤ 0.23), mid-low-risk area (0.23 < ERI ≤ 0.21), and low-risk area (ERI ≤ 0.21). The results are presented in Figure 4 and

Table 4. Area table of change in ecological risk level in 1987, 1997, 2001, 2013 and 2020.

Unit	Ecological Risk Levels	Low	Relatively Low	Middle	Relatively High	High
km2	1987	88.47	170.24	358.25	211.10	105.65
km2	1997	827.12	13.05	19.34	36.86	37.04
km2	2001	852.56	7.66	10.36	19.13	43.69
km2	2013	901.76	11.19	9.7	6.80	3.95
km2	2020	922.01	3.55	2.27	1.53	1.55

.

Overall (Figure 4 and Figure 5, Table 4), the landscape stability of Saihanba has been continuously improved, accompanied by a noticeable decline in landscape ecological risk. The expanse of low-risk areas has consistently expanded, while the areas of other risk grades have experienced a continual reduction over time. Examining the changes in landscape ecological risk areas across all levels in Saihanba, it is evident that in 1987, the medium-risk area constituted the largest portion, with medium-risk, mid-high-risk, and high-risk areas collectively accounting for 72.30% of the total area. The overall risk was relatively high, particularly in the central-south and northeast regions, which displayed elevated risk values. By 1997, with the exception of low-risk levels, the area of other risk levels was significantly reduced, with low-risk areas accounting for 88.55 percent of the total area, signifying a significant risk reduction. However, the risk shifted from the central-south and northeast to the west and northwest spatially. The western region was characterized by sandy land dominance, while the northwestern region was characterized by wetlands, resulting in relatively high landscape vulnerability. In 2001, the areas of both low-risk and high-risk areas experienced slight increases. The high-risk area expanded from the northwest to the north along the border line, and the regions in the south-midwest transformed from low-risk to high-risk areas. In 2013, the low- and mid-low-risk areas witnessed upward trends, increasing by 49.2 km² and 3.53 km², respectively. Spatially, the risk values in the western and northwestern regions exhibited effective improvement, while certain areas in the central and eastern regions showed enhanced risk. These areas contained farmlands, impacting the landscape’s stability, displaying a high degree of vulnerability and significant plaque fragmentation. By 2020, low-risk areas continued to grow, while other risk areas demonstrated a diminishing trend, indicating a considerable improvement in ecological risk across the study area.

3.3. Spatial Autocorrelation Analysis of Landscape Ecological Risk

3.3.1. Landscape Type Structure Transfer Analysis

Utilizing the temporal and spatial distribution data of landscape ecological risk in Saihanba from 1987 to 2020, a Moran’s I scatter plot was generated. This case is depicted in Figure 6. The landscape ecological risk values of Saihanba were 0.177, 0.120, 0.127, 0.205, and 0.072 in 1987, 1997, 2001, 2013, 2020 and 2020. Notably, the Moran’s I value peaked in 2013, signifying a relatively strong spatial agglomeration effect for landscape ecological risk. It suggests that the regional landscape ecological risk values have a significant impact on the adjacent regional landscape ecological risk values. Conversely, the Moran’s I value hit its lowest point in 1997, indicating a comparatively weak agglomeration effect. In comparison to other years, the regional landscape ecological risk value had a slightly less pronounced impact on the adjacent regional landscape ecology during this period.

3.3.2. Local Correlation Analysis

Conducting a more in-depth examination of the local spatial correlation of landscape ecological risk in Saihanba, we derived a local spatial autocorrelation clustering map (refer to Figure 7). Overall, the distribution of high-value areas in each year closely aligned with the spatial distribution observed in Kriging interpolation. In 1987, high-value areas were concentrated along the southwest boundary of the study area. The patches in this area were highly fragmented, mainly forests, with a small amount of wetland, and the degree of landscape loss was high. In contrast, low-value areas were scattered. By 1997, high-value areas shifted to the western region, with a minor distribution in the northwest. The high-value areas in the western region were dominated by sandy land, grasslands and forests. Some areas witnessed the transformation of grassland and sandy land into forests, resulting in patch fragmentation and diminished landscape stability. Sporadic low-value areas in the northeast, primarily composed of wetlands with high landscape vulnerability, were evident. In 1987, low-value areas were dispersed within high-value areas, suggesting that during this period, patch fragmentation underwent timely adjustments. Landscape patches in this area transitioned from low-fragility to high-fragility, reducing landscape loss. In 2001, high-value areas showed no significant change, but there was a noteworthy increase in low-value areas, particularly concentrated in the northeast. In 2013, the risk in high-value areas in the west and northwest became less prominent, with high-value areas exhibiting a sporadic distribution. In 2020, certain high-value areas increased in the central and western regions, predominantly composed of wetlands with relatively high landscape vulnerability. Concurrently, the distribution of low-value areas in the south transformed from being insignificant to concentrated.

3.4. Ecological Risk Analysis of Forest Landscape in 2020

3.4.1. Local Correlation Analysis

According to the forest landscape analysis of Saihanba in 2020 (Figure 8), the areas with mid-high and high-risk levels were primarily distributed in southern and northern regions, and the medium-risk regions were clustered in the southwestern, central, and western regions. From the global autocorrelation analysis, the ecological risk of forest landscape presented a clustering effect in space. Regarding local autocorrelation analysis, the spatial distribution of high-value areas and forest landscape risk levels was relatively consistent. The main distribution landscapes were large-scale larch forests and birch forests. The spatial distribution of Pinus sylvestris forests was primarily medium risk. The risk of the rest of the landscape is favorable.

3.4.2. Forest Landscape Risk Ecological Detection Results

Nine risk factors were classified using the natural breakpoint method, and the classification results are shown in the figure (Figure 9).

Risk Factor Detection and Analysis

The detection results, presented in the table, indicate that the impact strength of the risk factor on landscape ecological risk increases with a larger q value. Among the nine risk factors, the intensity of detection results of forest landscape risk factors in Saihanba is as follows: soil type > precipitation > forest landscape type > air temperature > slope > DEM > slope aspect > subcompartment volume > vegetation cover grade. From the degree of significance of differences, soil type, precipitation, forest landscape type, air temperature, slope, DEM, and slope aspect have extremely significant differences in forest landscape risk, other factors have no significant difference on forest landscape risk, as shown in

Table 5. Forest risk factor detection results.

Risk Factor	Soil Type	Precipitation	Forest Landscape Type	Air Temperature	Slope	DEM	Slope Aspect	Subcompartment Volume	Vegetation Cover Grade
q statistic	0.228	0.133	0.076	0.056	0.056	0.025	0.016	0.014	0.001
p value	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.039	0.999

.

Analysis of Risk Ecology Detection

The determination of significant differences in the spatial distribution of risk factors influencing landscape ecological risks is accomplished through ecological detection. The test, conducted at a significance level of 0.05, assesses whether two factors exhibit significant differences in influencing the distribution of landscape ecological risks. The detection results, depicted in Figure 10, reveal notable distinctions in landscape ecological risk between soil type and other factors, as well as between precipitation and other factors. Additionally, there are significant differences in other factors, such as DEM, vegetation cover grade, slope aspect, subcompartment volume and forest landscape type;slope, air temperature, and vegetation cover grade, which together play a role in landscape ecological risk. However, there is no significant difference between other factors and landscape ecological risk, indicating that each risk factor has a lesser impact on forest landscape ecology in the study area, maintaining overall stability.

Risk Interaction Detection Analysis

Using the method of interactive detection, we can detect whether different risk factors have interactive effects on the spatial distribution of landscape ecological risks (

Table 6. Forest risk interaction detector results.

	Forest Landscape Type	DEM	Vegetation Cover Grade	Precipitation	Slope	Slope Aspect	Air Temperature	Soil Type	Subcompartment Volume
Forest landscape type	0.076
DEM	0.136 #	0.025
Vegetation cover grade	0.081 #	0.035 #	0.001
Precipitation	0.204 *	0.203 #	0.143 #	0.133
Slope	0.139 #	0.178 #	0.067 #	0.241 #	0.056
Slope aspect	0.113 #	0.077 #	0.030 #	0.182 #	0.109 #	0.016
Air temperature	0.219 #	0.160 #	0.077 #	0.461 #	0.209 #	0.124 #	0.056
Soil type	0.265 *	0.313 #	0.237 #	0.345 *	0.265 *	0.253 #	0.451 #	0.228
Subcompartment volume	0.096 #	0.059 #	0.021 #	0.161 #	0.089 #	0.051 #	0.100 #	0.246 #	0.014

Note: * denotes double factor enhancement, # denotes nonlinear enhancement.

). The interaction between the two factors was significantly higher than the single factor alone. Among them, the nonlinear enhancement effect was more significant than the double-factor enhancement. It also showed that there was no independent effect between the factors, indicating the impact of forest landscape ecological changes was the result of multiple factors. Forest landscape type and precipitation, forest landscape type and soil type, precipitation and soil type, slope and soil type have double-factor enhancement effects, and the other factors are nonlinear enhancement effects, which alone affect the spatial distribution of ecological risks.

4. Discussion

4.1. Landscape Pattern and Spatial Scale Change in Landscape Ecological Risk

In this research, we used NDVI and the transition matrix method, and found that the landscape pattern of Saihanba has undergone substantial changes from 1987 to 2020, and the vegetation cover grade area has increased significantly. Among diverse landscape types, there was an upward trajectory in the area of forests and wetlands, while the areas of grassland, sandy land, and construction land exhibited a declining trend.

In 1987, the overall risk was relatively high, and the risk was significantly reduced in 1997. This was related to the decrease in grassland and sand and the substantial increase in forest area from 1987 to 1997. Before 1987, China began extensive afforestation: as an important strategic location in China, a large number of forests were planted in Saihanba. Until 1997, tree growth resulted in increased crown density, so that through remote sensing imagery we found that the landscape changed markedly during this period, and this change was mainly due to the trees that were planted on grasslands and sandy land. Based on the spatial correlation analysis, the changes to landscape risk within the five periods all showed a positive correlation trend. This was consistent with the results of many scholars that changes in land use have a significant effects on landscape ecological risks [1,55,56,57].

4.2. Evaluation Method of Landscape Ecological Risk Index

The assessment of regional landscape ecological risk primarily relies on landscape pattern index evaluation [58,59] and ecosystem service value [60,61]. Given the significant changes in the landscape pattern over more than 40 years in the study area, the variation in the landscape ecological risk index is primarily attributed to alterations in the landscape pattern. The landscape ecological risk index exhibited a continuous decrease from 1987 to 2020. Notably, from 2013 to 2020, the high-risk area of landscape ecology remained relatively stable with no noticeable transfers, indicating an inward contraction trend in the high-risk area. However, diverse quantitative methods yielded different outcomes for the landscape ecological risk assessment, and the subjective nature of landscape vulnerability classification impacted the spatial distribution of landscape ecological risk. Consequently, the assessment of landscape ecological risk in this study remains relatively one-sided. Moving forward, it is imperative to further account for the influence of additional factors on landscape ecological risk.

Building on the evaluation of the landscape ecological risk index from 1987 to 2020, a more focused investigation into the risk of forest landscapes in 2020 was conducted. Through a comprehensive analysis of landscape ecological risk in the study area, the vulnerability of the forest landscape was determined based on the risk level associated with each forest landscape, combined with the expert consultation method. The findings revealed that mid-high and high-risk areas were dispersed in the south, northeast, and central parts of the study area, primarily characterized by extensive larch forests and birch forests. This study suggests the need to adjust stand structures within a reasonable range in these areas.

4.3. Analysis of Forest Landscape Risk Drivers and Suggestions for Future Development

Addressing the ecological risk status of the forest landscape in 2020, the geographical detector method was employed to delve into its governing factors. Recognizing that the internal structure of the forest landscape and various geographical factors exert a substantial impact on landscape risk, nine factors with a significant influence on the Saihanba forest landscape were chosen based on the management status of the Saihanba Mechanical Forest Farm. The study revealed that the forest landscape type poses the greatest risk to landscape ecology, aligning with the findings of numerous experts. Simultaneously, this observation is connected to the classification method of landscape vulnerability. Currently, research on landscape ecological risk predominantly focuses on river basins and coastlines, leaving a research gap in the context of forest landscapes [62]. Therefore, building upon the landscape risk in the study area, this study compared the risk levels of each forest landscape and classified landscape vulnerability based on expert advice.

Soil type, precipitation, forest landscape type, temperature, slope, DEM, and slope direction have a significant impact on forest landscape ecological risk. Among these factors, there are significant differences in soil type and precipitation compared to other factors. Hence, this area should be designated as a pivotal zone for risk prevention to avert an upward trend in risk. Analyzing the factors that have a pronounced impact on the Saihanba Mechanical Forest Farm, including landscape type, slope, slope direction, precipitation, vegetation cover grade, and subcompartment volume, reveals that, apart from fixed geographical factors (slope, slope direction, precipitation), the vegetation cover grade has reached an optimal level, with low-value areas sparsely distributed, mainly found in water areas and construction land. Subcompartment volume is related to the forest landscape type and forest management. Therefore, formulating a reasonable forest management plan and gradually adjusting the forest landscape type, specifically the stand structure, is imperative. Currently, the stand structure is monotonous, necessitating the introduction of diverse tree species. Simultaneously, adhering to the principle of suitable tree species and updating the forest patch configuration, as suggested by this study, not only contributes to mitigating the risk of the forest landscape but also holds significant ecological importance for forest carbon storage and enhancing forest site quality.

Utilizing the landscape ecological risk index and spatial autocorrelation analysis, this study assessed the ecological risk pattern of the Saihanba Mechanical Forest Farm from the viewpoint of landscape spatial structure. Additionally, a geographic detector was employed to analyze risk-driving factors. From the results, the effect of multi-factor combination is not simply one-plus-one, but far greater than the single factor, so it should be fully considered in the management and development of forest farms. Since the establishment of Saihanba, extensive tree planting efforts have significantly improved the ecological environment, resulting in substantial changes to the landscape pattern. Given that landscape pattern alterations inherently impact ecosystem functions, utilizing landscape pattern indices to investigate the ecological risk pattern is both feasible and effective. Consequently, it is recommended that the Saihanba Mechanical Forest Farm maintains and adjusts its existing landscape structure and addresses high-risk areas in alignment with its specific conditions. At the same time, to pay attention to the diversity of forest structure and its spatiotemporal variation characteristics is necessary [62]. The enhancement of the landscape integration pattern should prioritize natural factors. During resource development, adhering to the principle of ecological priority is crucial. While safeguarding river waters and promoting forestry development, it is imperative to regulate the uncontrolled development of local tourism. Zoning management, aligned with the current state, should be implemented, coupled with the establishment of an effective ecological compensation system.

5. Conclusions

In this study, the landscape ecological risk of the Saihanba Mechanical Forest Farm was comprehensively evaluated from four aspects: land use change and transfer, spatial-temporal distribution, spatial correlation of landscape ecological risk, and risk factor detection. The results are shown as follows:

(1)

Between 1987 and 2020, the Saihanba forest landscape area increased, along with a corresponding rise in wetland coverage. Conversely, grassland, sandy land, and construction land areas experienced a decrease. The forest landscape emerged as the predominant type in Saihanba, exerting a significant influence on the overall landscape pattern change.

(2)

Between 1987 and 2020, the Saihanba area experienced a continuous reduction in the extent of medium-, mid-high-, and high-risk levels, coupled with a significant increase in the low-risk area. Spatially, the overall risk level was relatively high in 1987. In the period from 1997 to 2001, the risk level demonstrated a significant decrease, but there was a spatial shift in the risk areas. During this time, the high-risk area was limited to the border between the west and northwest, indicating an overall decreasing trend of landscape ecological risk. From 2013 to 2020, the landscape predominantly comprised low-risk areas. The medium-, mid-high-, and high-risk areas were confined to a small, aggregated distribution, without outward spread. The overall landscape ecological risk tended to stabilize during this period.

(3)

Throughout the study period, the global autocorrelation Moran’s I values of the watershed were 0.177, 0.120, 0.127, 0.205 and 0.072 respectively. The landscape ecological risks were positively correlated and tended to be clustered in spatial distribution. The landscape ecological risk level showed greater alignment with the high and low values of the local autocorrelation pattern.

(4)

Analyzing the landscape ecological risk distribution map in 2020 enabled a detailed examination of the risk associated with various forest landscapes. The vulnerability was classified based on the risk value of each forest landscape. It was observed that the southern and northern regions exhibited higher risk levels, predominantly characterized by pure forests. Mitigating the risk of the forest landscape in these areas requires measures to enhance the stand structure.

(5)

In the analysis of regional risk factors for each landscape risk level in 2020, it was determined that soil type, precipitation, forest landscape type, temperature, slope, DEM, and slope aspect exerted significant influence on forest landscape ecology risk. Among these factors, soil type had the most substantial impact, and its interaction with precipitation, as well as with other factors, played a crucial role. Notably, the interactions between forest landscape type and precipitation, forest landscape type and soil type, precipitation and soil type, and slope and soil type were considerably stronger than the influence of each factor in isolation.