Multi-Factor Collaborative Analysis of Conservation Effectiveness of Nature Reserves Based on Remote Sensing Data and Google Earth Engine

Abstract

Protected areas (PAs) play a crucial role in safeguarding biological resources and preserving ecosystems. However, the lack of standardized and highly operational criteria for evaluating their conservation effectiveness, particularly across different ecological types, remains a significant gap in the literature. This study aims to address this gap by constructing a conservation effectiveness evaluation model for two distinct types of PAs in Heilongjiang Province, China: the Zhalong National Nature Reserve (ZlNNR), a wetland ecological reserve; and the Mudanfeng National Nature Reserve (MdfNNR), a forest ecological reserve. We employed various methods, including land use dynamic index, visual analysis of landscape patterns, remote sensing inversion, and a multi-factor comprehensive assessment model, to assess changes in conservation effectiveness from 2000 to 2020. Our findings reveal a contrast between the two PAs. In the ZlNNR, croplands and water bodies increased significantly by 4069.4 ha (K = 1.5820%) and 2541.58 ha (K = 3.2692%). In the MdfNNR, impervious lands increased greatly by 65.35 ha (K = 7.4021%), whereas forest lands decreased by 125 ha (K = −0.067%). The core area of the two PAs displayed increased landscape regularity, whereas the experimental area showed heightened landscape diversity. In ZlNNR, the MPS_L value increased by 134.91%, whereas the PD_L value decreased by 57.43%, indicating a more regular landscape pattern. In MdfNNR, the SHDI_L value decreased by 110.7%, whereas the PD_L value increased by 52.55%, indicating a more fragmented landscape pattern. The area with improved vegetation trends in ZlNNR was 8.59% larger than in MdfNNR, whereas the area with degraded vegetation trends was 4.86% smaller than in MdfNNR. In all years, the high effectiveness area was larger in ZlNNR than in MdfNNR, whereas the medium and low effectiveness areas were smaller in ZlNNR compared to MdfNNR. This study provides a scientifically rigorous assessment method for evaluating the conservation effectiveness of different types of PAs, laying a solid theoretical foundation and practical guidance for future conservation strategies.

1. Introduction

Protected areas (PAs), as delineated in scholarly literature, are unique geographical zones that are not only legally recognized and managed, but are also committed to the enduring preservation of nature, inclusive of the associated ecosystem services and cultural values [1]. At present, the establishment of PAs is a pivotal global strategy for biodiversity conservation, climate change mitigation, and the reduction of human impact [2,3,4,5,6,7]. It also serves as a significant policy instrument for the protection of natural resources and the promotion of sustainable development [1,8,9]. China, being one of the countries with the most diverse ecosystems globally, has instituted more than ten categories of PAs, which encompass 18% of its terrestrial area and 4.1% of its marine area [10,11,12]. However, numerous PAs are still under considerable pressure, confronting challenges such as the exploitation of agricultural and forestry products, illegal hunting, and infrastructure development [13]. In essence, merely augmenting the quantity and size of PAs is no longer adequate to cater to the present and future requirements of biodiversity conservation [12,14]. The most pressing concern at this juncture is to amplify the effectiveness of existing PAs to ensure their preservation [15].

The effectiveness of PAs is contingent on a multitude of political, economic, and local environmental factors. There are authoritative domestic and international assessment standards for PAs, such as the Standard for Conservation Effectiveness Assessment of Ecology and Environment in Nature Reserves, promulgated by the Ministry of Ecology and Environment of China, and Evaluating Effectiveness: A Framework for Assessing Management Effectiveness of Protected Areas, published by the International Union for Conservation of Nature (IUCN), Gland, Switzerland. In practice, the assessment criteria are relatively sophisticated, necessitating specialized technical expertise, rendering implementation challenging. Furthermore, numerous studies have utilized multi-temporal historical land use data to appraise the land use/land cover changes and landscape pattern dynamics within PAs [16,17,18]. Indexes such as the landscape development intensity index and landscape pressure index have been employed by academics to assess the effectiveness of safeguarding marine PAs [19]. In addition, alterations in ecosystems, as well as the distribution and abundance of key species, also serve as indicators of conservation effectiveness [20,21,22]. Nevertheless, these existing evaluation techniques are often labor-intensive and time-consuming [23,24], or necessitate plenty of field investigations [21].

Remote Sensing (RS) is an advanced data acquisition technology that can swiftly, comprehensively, and accurately generate thematic data. It facilitates repeated observations of large territories and inaccessible areas, proving to be more cost-effective than field observations [25,26]. In recent years, research using satellite imagery to predict sustainable development outcomes has achieved rapid growth, especially studies employing machine learning methods. Satellite data with greatly improved spatial, temporal, and spectral resolutions have provided richer information on the Earth’s surface. RS meets the needs for monitoring sustainable development outcomes, with the advantages of mature technology, rapid data updates, and extensive coverage, which can satisfy the requirements of large-scale, high-precision effectiveness monitoring in PAs [27]. The Geographic Information System (GIS), with its robust information extraction technology and spatial analysis capabilities, can accurately and efficiently address geographic issues related to resources and the environment [28]. The significance of this technology continues to escalate, with ongoing enhancements to its role and functionality. Consequently, the integration of RS and GIS has emerged as the preferred methodology for evaluating the conservation effectiveness of PAs.

From our perspective of PAs as habitats for various species, the quality trends, whether improving or deteriorating, are an important reflection of the changing effectiveness of conservation efforts. We utilize an integrated approach of RS and GIS techniques to conduct monitoring of PA quality changes. This allows for more cost-effective and extensive spatio-temporal analyses compared to conventional methods based on localized ground surveys. Furthermore, our multi-factorial integrated assessment model employs a weighted summation method, which can effectively mitigate the limitations of single-factor evaluation lacking accuracy. This represents a methodological innovation in our research. These represent the innovations in our research ideology and analytical methodology. Thus far, many appraisals of conservation effectiveness have focused solely on either a single category of PA [21,29] or involved a plethora of disparate PA categories [19,30,31,32,33,34]. Therefore, our study aims to analyze differences in the effects of various PA types through a comparative analysis of two major types (forest and wetland PAs). By concentrating the analysis on these two predominant habitat types, we can gain deeper insights into the relative performance of PAs designated for forests versus wetlands in terms of conservation outcomes and effectiveness.

Land use change is a pivotal factor driving global environmental change [35,36]. It provides a direct perspective on the impact of human activities on surface ecosystems [37], effectively mirroring the ecosystem status of PAs [38]. Landscape pattern serves as a vital means to depict landscape structure [39]. Landscape indexes encapsulate essential information about landscape patterns, and can serve as macro-scale surrogate indicators to perceive and comprehend changes in species and their habitats [21]. Landscape indexes have been extensively employed to assess the effectiveness of PAs [40,41,42]. The spatial and temporal variation of the Normalized Difference Vegetation Index (NDVI) is one of the most cost-effective methods for assessing ecological effectiveness [34,43,44,45]. The Google Earth Engine (GEE) Cloud Platform can swiftly and accurately process remote sensing images using its robust data calculation function. It can also obtain the best-quality images through minimum cloud amount screening for research [46,47,48], thereby yielding the results of time series changes of NDVI. Consequently, employing the aforementioned factors to evaluate the effectiveness of PAs in a more systematic manner holds significant guiding implications.

Frequently, forests are often referred to as the “lungs of the earth”, playing a pivotal role in the global carbon cycle, regulating the climate system, and protecting biodiversity [49,50,51,52]. Wetlands, on the other hand, are commonly referred to as the “kidney of the earth”, providing a habitat for a vast array of species [53], maintaining the ecological balance [54,55,56], and mitigating the impact of extreme weather and other outside destructive forces on surrounding areas [57]. Both ecosystems are critically important, yet are under severe threat [51,53]. The establishment of forest ecosystem PAs and wetland ecosystem PAs is a key approach to protecting forests and wetlands, and a universally powerful tool to reverse the rapid imbalance and degradation of these two ecosystems [58]. Given disparities in the developmental status and geographic contexts of these two PA types, they may encounter differing degrees of threat and exhibit variations in biodiversity composition and structure [30]. Consequently, their responses and outcomes regarding conservation management could diverge as well. Therefore, to what extent do differences exist in the conservation effectiveness between these two types of PAs? What are the driving factors contributing to these differences? Through comparison, we can ascertain relative superiority in habitat quality conservation outcomes, elucidate the adaptive propriety of conservation management across contrasting ecosystems, and evaluate the necessity of establishing specialized PAs grounded in particular ecological systems.

In light of the above, the objectives of this study are delineated into the following components: (1) the development of a conservation effectiveness evaluation model that incorporates parameters such as land use, landscape pattern, and NDVI; (2) the exploration of the differences in conservation effectiveness between Zhalong National Nature Reserve (ZlNNR), Qiqihaer, China and Mudanfeng National Nature Reserve (MdfNNR), Mudanjiang, China from the period of 2000 to 2020; (3) the identification of the primary driving factors influencing the conservation effectiveness of each PA, and the provision of scientifically sound and reasonable suggestions for the effective management of PAs.

2. Materials and Methods

2.1. Study Areas

The ZlNNR study area is situated in the Songnen Plain in the west of Heilongjiang Province, southeast of Qiqihar City (Figure 1), within the reed marsh belt in the lower reaches of the Wuyuer River. It belongs to the wetland ecosystem type, and focuses on protecting rare bird species such as red-crowned cranes and wetland ecosystems [59,60]. It is the most complete, primitive, and open wetland ecosystem in the same latitude area of northern China [61]. The reserve, irregularly shaped like an olive and stretching from northeast to southwest, spans a total area of 2100 km². Its geographical coordinates are 123°47′~124°37′E and 46°52′~47°32′N. The region experiences a temperate continental monsoon climate, with an annual average temperature of 2~4 °C, annual precipitation of approximately 426 mm, and annual evaporation of 1307.7 mm [61,62].

The MdfNNR study area is located in the low-altitude Laoyeling mountains, southeast of Mudanjiang City (Figure 1), and is 15 km away from the city center [63]. It is a forest ecosystem, covering a total area of 210 km², with geographical coordinates of 129°40′~129°53′E and 44°20′~44°30′N. The region is characterized by a temperate continental monsoon climate, featuring dry and windy springs, followed by warm and wet summers. The average annual temperature is a moderate 3.6 °C, whereas the annual precipitation averages at 550 mm. The majority of this precipitation occurs during the growing season, which typically spans from June to August [64].

Both the NNR and buffer (ZlNNR: 10-km-wide, MdfNNR: 3-km-wide) around them are the concerned areas in this paper.

2.2. Data Sources

In this study, land cover/land use data were derived from the GlobeLand30 data platform [65,66,67]. These maps have recently been considered useful references for land use studies and have been used to classify land use [68]. The remote sensing image data use the GEE platform programming (JavaScript API) [69] to call the Landsat-5 TM, Landsat-7 ETM+, and Landsat-8 OLI surface reflection data set to screen the cloud amount of the whole year’s image. A total of 2012 available Landsat images were obtained across the two PAs and their surroundings. The images are then mosaicked to find the minimum cloud cover image of the research area from 2000 to 2020. Notably, since the Landsat-5 data set only contains the images from 2003 to 2011, and the Landsat-7 satellite failed in 2003, images after that year cannot be used and the images from 2012 are discarded during the analysis.

2.3. Method

Figure 2 shows the specific processes of assessing the conservation effectiveness and the included indexes.

2.3.1. Land Use Change

The analysis of land use change is conducted from two perspectives: temporal series variation and differences inside and outside PAs. The ArcGIS 10.8 software (ESRI (Environmental Systems Research Institute). 2018. ArcGIS Desktop: Release 10.8. Redlands, CA, USA: Environmental Systems Research Institute. https://www.esri.com/en-us/arcgis/products/arcgis-desktop/overview, accessed on 9 April 2023) is used to make statistics on the land use area and transfer matrix characteristics to analyze the process and degree of change among different land use types in the study areas. The dynamic degree model of land use is used to analyze the change degree of various types of land use [21,70,71]. The dynamic degree index (K) is calculated as follows:

$$K=\frac{U_a-U_b}{U_b}\times \frac{1}{T}\times 100\%$$

where, K refers to the dynamic degree of a certain land use type during the study period; U_b and U_a are the areas of a certain land use type before and after the study period, respectively (unit: ha); and T is the number of years between the study period.

2.3.2. Landscape Pattern

The landscape pattern index value can effectively encapsulate information pertaining to landscape development. It can accurately reflect the composition and dynamic changes of the landscape structure, and closely aligns spatial changes of the landscape with the temporal process [72]. In this study, eight indexes were selected from the patch level and landscape level to analyze the changes in landscape pattern, referencing previous studies. Additionally, these indexes at the patch level were categorized into two groups: patch characteristics and patch shape. At the landscape level, they were divided into three categories: landscape characteristics, landscape shape, and landscape diversity [73]. The calculation formula for each index is presented in

Table 1. The index of landscape pattern.

Landscape Index	Patch-Level	Landscape-Level	Formula	Ecological Significance
Patch density (PD)	patch characteristics	landscape characteristics	$PD=\frac{{{\displaystyle \sum }}_{i=1}^M{N}_i}{A} \left(PD>0\right)$	PD [78] is used to describe the fragmentation and uniformity of land use in the study area.
Largest patch index (LPI)			$LPI=\frac{Max\left(a_1{a}_2\dots \dots a_n\right)}{A}\times 100$	LPI [78] is used to reflect the impact of the patch with the largest area on the landscape pattern in the study area.
Main patch size (MPS)			$MPS=\frac{A_i}{N_i}\times {10}^6$	MPS [79] is used to reflect the degree of landscape aggregation and fragmentation of various landscape types in the study area.
Edge density (ED)	patch shape	landscape shape	$ED=\frac{L_i}{A_i}\times {10}^6$	ED [73] is used to reflect the complexity of plaque boundary in the study area.
Landscape shape index (LSI)			$LSI=\frac{0.25L_i}{\sqrt{A_i}}$	LSI [78] reflects the complexity of patch shape and the possible evolution trend of shape characteristics of landscape spatial structure. The larger the value, the more complex and irregular the patch shape.
Contagion (CONTAG)			$CONTAG=1+\frac{{{\displaystyle \sum }}_{i=1}^n{{\displaystyle \sum }}_{j=1}^n{P}_{ij}{log}_2{P}_{ij}}{2lnn}$	CONTAG [80] is used to describe the connectivity between patches in the study area.
Shannon’s diversity index (SHDI)		landscape diversity	$SHDI=-{\displaystyle \sum _{i=1}^m}\left(p_i\times ln{p}_i\right)$	The higher the SHDI [81], the higher the fragmentation degree, indicating that the distribution of each patch type is more balanced.
Shannon’s evenness index (SHEI)			$SHEI=\frac{-{{\displaystyle \sum }}_{i=1}^m\left(p_i\times ln{p}_i\right)}{lnm}$	SHEI [73] is used to describe the dominance of one or several landscape types in the landscape structure.

Note: N_i refers to the number of patches of the i landscape type. A is the total area of the study area. M represents the number of landscape types in the study area. A_i is the total area of patch type i, and L_i represents the total length of the patch boundary in the study area. n is the total number of landscape types. P_ij represents the probability that the i and j land use types are adjacent. m is the total number of type i patches. N is the total number of patches in the landscape. A_ij represents the patch area of the i and j landscape types in the study area and N_ij represents the patch perimeter of the landscape type. p_i is the proportion of the area occupied by landscape type i.

. The grid analysis method was employed to visually analyze the landscape level of the landscape pattern indexes of both NNRs. Primarily, referring to the relevant literature [74,75,76,77], considering the area and scale of each study area, ZlNNR and MdfNNR were divided into 2 km × 2 km and 0.5 km × 0.5 km grid cells, respectively. Furthermore, the landscape indexes of each study area were interpolated and analyzed using the Kriging interpolation method.

2.3.3. Trend Analysis of NDVI

NDVI is suitable for characterizing the growth and spatial distribution of vegetation within our study area [82,83]. We extracted the annual average NDVI value for each pixel within the study areas through the GEE platform [84,85]. However, when conducting a univariate analysis of NDVI, there are several considerations that must be taken into account. NDVI may exhibit anomalous values in non-vegetated regions. Furthermore, water bodies result in negative NDVI values. Therefore, when performing an isolated NDVI analysis, croplands, water bodies, barren areas, and impervious areas were excluded from the study areas. This was done to avoid potentially skewed results arising from these land cover types which do not accurately represent vegetation dynamics. The time series analysis of NDVI was conducted pixel by pixel using the Theil–Sen median trend analysis [86,87] and Mann–Kendall trend test [88,89]. The Theil–Sen Median method is a robust non-parametric statistical trend calculation method. This method boasts high computational efficiency, is insensitive to measurement errors and clustered data, and is apt for the trend analysis of long-term series data. The Mann–Kendall test is a non-parametric time series trend test method, which does not require measurement values to conform to a normal distribution and is unaffected by missing values and outliers. It is suitable for the trend significance test of long-time series data. Both methods are vital for determining the trend of long-time series data and can be effectively combined. This combination has gradually been used to analyze and reflect the trend change of each pixel in the long-time series [84,90,91,92,93].

The calculation formula of NDVI is:

$$NDVI=\frac{NIR-R}{NIR+R}$$

where NIR refers to near-infrared band value. R represents red band value.

The calculation formula of the Theil–Sen median method is:

$$\beta =Median\left(\frac{NDV{I}_j-NDV{I}_i}{j-i}\right),\forall j>i$$

where Median() represents the median value, and NDVI_i and NDVI_j represent the NDVI values in years i and j. If β > 0, the NDVI is on an upward trend; otherwise, it is a downward trend.

The calculation formula of the Mann–Kendall test is as follows. The calculation formula of test statistic S is:

$$S=\sum _{i=1}^{n-1}\sum _{j=i+1}^{n}sgn(NDV{I}_j-NDV{I}_i)$$

where sgn() is a symbolic function, and the calculation formula is:

$$sgn\left(NDV{I}_j-NDV{I}_i\right)=\{\begin{array}{c}+1 NDV{I}_j-NDV{I}_i>0\\ 0 NDV{I}_j-NDV{I}_i=0\\ -1 NDV{I}_j-NDV{I}_i<0\end{array}$$

The calculation formula of test statistic Z is:

$$Z=\{\begin{array}{c}\frac{S}{\sqrt{Var\left(S\right)}}(S>0)\\ 0 \left(S=0\right)\\ \frac{S+1}{\sqrt{Var\left(S\right)}}(S<0)\end{array}$$

where the formula for calculating Var is:

$$Var\left(S\right)=\frac{n\left(n-1\right)\left(2n+15\right)}{18}$$

where n is the length of data in the sequence.

2.3.4. Evaluation and Analysis of Conservation Effectiveness

The evaluation indicator system is presented in

Table 2. Evaluation index system.

Level 1 Evaluation Indexes	Level 2 Evaluation Index
Land use types	Cropland
	Forest land
	Grassland
	Wetland
	Water
	Impervious
	Barren
Landscape pattern indexes	PD
	LPI
	MPS
	ED
	LSI
	CONTAG
	SHDI
	SHEI
NDVI	NDVI

. We selected evaluation indicator data for the years 2000, 2010, and 2020 to conduct a comprehensive evaluation of conservation effectiveness. The weight of the secondary evaluation index was calculated first, followed by the weight of the primary evaluation index using the coefficient of variation results of the secondary evaluation index. The coefficient of variation method [93,94] was employed to determine the weight of each indicator. The formula for calculating the coefficient of variation method is:

$$v_i={\sigma }_i/\overline{x_i}$$

$$W_i=v_i/{\displaystyle \sum }_{i=1}^n{v}_i$$

where W_i refers to the weight of index i; v_i represents the coefficient of variation of index i; σ_i represents the standard deviation of index i; $\overline{x_i}$ is the average of index i.

We then standardized the evaluation factors between [0–10] and used the relative evaluation method to evaluate the conservation effectiveness of the NNRs [95]. The indicator data of the NNRs were collated into a database for use. According to the weight of each evaluation indicator, all evaluation indexes were spatially superimposed by the method of weighted summation to calculate the conservation effectiveness indicator of the NNRs. Subsequently, we employed an equal spacing segmentation method and divided the NNRs into three level intervals, high effectiveness (HE), medium effectiveness (ME), and low effectiveness (LE), to realize the comprehensive evaluation of the conservation effectiveness of the NNRs [96]. The model calculation formula is as follows [97]:

$$E={\displaystyle \sum }_{i=1}^n{W}_i\times X_i$$

where E is the conservation effectiveness index of the study area; X_i refers to the quantitative expression value of the evaluation index i element after standardization; and n is the number of thematic index elements participating in the evaluation.

3. Results

3.1. Land Use

Wetland constitutes the primary land use type in ZlNNR, accounting for approximately 50% of the total area (Figure 3a and

Table 3. The change of land use in ZlNNR from 2000 to 2020.

Inside and Outside	Type	Area/ha			Change Rate (K)/%
		2000	2010	2020	2000–2010	2010–2020	2000–2020
Inside	Cropland	25,723.60	25,720.10	29,793.00	−0.0014	1.5835	1.5820
	Forest land	19,330.00	19,324.30	17,881.80	−0.0029	−0.7465	−0.7492
	Grassland	26,524.90	26,539.20	27,133.90	0.0054	0.2241	0.2296
	Wetland	127,938.00	127,936.00	122,823.00	−0.0002	−0.3997	−0.3998
	Water	7774.32	7774.40	10,315.90	0.0001	3.2691	3.2692
	Impervious	1732.97	1731.23	2200.61	−0.0100	2.7113	2.6985
	Barren	4507.66	4506.79	3384.49	−0.0019	−2.4902	−2.4917
Outside	Cropland	124,200.00	124,191.00	124,245.00	−0.0007	0.0043	0.0036
	Forest land	20,653.00	20,649.60	16,738.40	−0.0016	−1.8941	−1.8954
	Grassland	61,684.80	61,687.50	63,247.20	0.0004	0.2528	0.2533
	Wetland	17,125.80	17,127.40	16,480.00	0.0009	−0.3780	−0.3771
	Water	15,110.80	15,111.30	16,120.50	0.0003	0.6678	0.6682
	Impervious	6769.50	6768.12	10,619.10	−0.0020	5.6899	5.6867
	Barren	9638.28	9646.68	7726.32	0.0087	−1.9907	−1.9837

). Over the span of twenty-one years, the wetlands experienced a reduction of 5115 ha (K = −0.3998%). However, the areas of croplands and water bodies significantly increased, by 4069.4 ha (K = 1.5820%) and 2541.58 ha (K = 3.2692%), respectively. The change in grasslands was the least significant (K = 0.2296%). Outside the NNR (Figure 3b), impervious lands (K = 5.6867%) underwent the largest change, whereas the croplands, being the primary land use type outside of ZlNNR, exhibited the smallest change range (K = 0.0036%).

A total of 25,481.16 ha and 43,366.68 ha of land underwent transformation inside and outside ZlNNR, accounting for approximately 11.94% and 17.01% of their respective total areas (Figure 4a,b). Among these, inside the ZlNNR, the area of grasslands experienced the largest transformation, with 8115.78 ha transferred out, accounting for 30.60% of the total area of grasslands, primarily converted into forest lands (3410.08 ha) and croplands (2894.82 ha). The area transferred from impervious lands was the smallest, with 285.76 ha transferred out, accounting for 16.51% of the total area of impervious lands, mainly converted into croplands (177.65 ha). Outside the ZlNNR, the area of grasslands transferred out was the largest (14,932.17 ha), accounting for 24.23% of the total grasslands area, and primarily converted into croplands (6615.25 ha). The area transferred from impervious lands was the smallest (877.59 ha), accounting for 12.98% of the total area of impervious lands, and mainly converted to croplands (577.69 ha). The area of wetlands both inside and outside of ZlNNR has declined.

Forest lands constitute the primary land use type both inside and outside of MdfNNR, covering 86.98% of the total area inside the reserve (Figure 3c,d and

Table 4. The change of land use in MdfNNR from 2000 to 2020.

Inside and Outside	Type	Area/ha			Change Rate (K)/%
		2000	2010	2020	2000–2010	2010–2020	2000–2020
Inside	Cropland	2534.86	2594.20	2591.64	0.2341	−0.0099	0.2240
	Forest land	18,644.40	18,514.40	18,519.40	−0.0697	0.0027	−0.0670
	Grassland	6.73	74.55	8.13	100.8333	−8.9094	2.0871
	Water	16.39	16.46	18.29	0.0428	1.1098	1.1573
	Impervious	88.28	91.08	153.62	0.3174	6.8667	7.4021
Outside	Cropland	5462.07	5695.66	5600.35	0.4277	−0.1673	0.2532
	Forest land	17,175.10	16,998.50	16,943.00	−0.1028	−0.0326	−0.1351
	Grassland	10.16	45.54	26.28	34.8274	−4.2290	15.8699
	Water	28.87	31.67	32.73	0.9708	0.3351	1.3384
	Impervious	475.86	380.72	551.56	−1.9994	4.4873	1.5907

). Inside the MdfNNR, the area of impervious lands exhibited the largest increase of 65.35 ha (K = 7.4021%). Conversely, the area of grasslands showed the smallest increase, expanding by a mere 1.4 ha (K = 2.0871%). However, forest lands were the only land type that decreased in size, with a total reduction of 125 ha (K = −0.067%). Outside the MdfNNR, only forest lands decreased in area, with a total reduction of 232.1 ha (K = −0.1351%). The area of croplands underwent the most significant change, 138.28 ha (K = 0.2532%). The most substantial change was observed in grasslands (K = 15.8699%).

In total, 761.04 ha and 1582.36 ha of land inside and outside of MdfNNR underwent changes, accounting for approximately 3.58% and 6.85% of the total area inside and outside of the NNR, respectively (Figure 4c,d). Among these, the largest transfer inside the MdfNNR was observed in forest lands, with an area of 392.26 ha and a ratio of 2.11% of the total forest lands area, primarily converted into croplands (382.24 ha). The area of grasslands transferred was the least (5.96 ha), with a ratio of 88.54% of the total grasslands area, all of which were converted into forest lands. Outside the MdfNNR, the area of forest lands transfer was the largest (740.58 ha), with a ratio of 4.33%, which was mainly converted into croplands (690.47 ha). The grasslands contributed the smallest area change of 9.25 ha, but the largest conversion ratio of 91.03% of its total area, mainly transformed into forest lands (8.13 ha).

3.2. Landscape Pattern

The patch level of the landscape pattern indexes of both NNRs are presented in

Table 5. Landscape index of patch metrics ZlNNR in 2000, 2010, and 2020.

Inside and Outside	Type	Year	PDP	LPIP	MPSP	EDP	LSIP
Inside	Cropland	2000	0.1222	2.1914	98.5048	6.6391	23.5042
		2010	0.0574	1.0334	98.8698	3.3280	23.4995
		2020	0.0519	1.2677	126.7836	3.0631	20.0990
	Forest land	2000	1.2982	1.1861	6.9728	17.6603	68.3538
		2010	0.6141	0.5593	6.9472	8.3920	68.3333
		2020	0.5275	0.6785	7.4879	7.5987	64.3016
	Grassland	2000	2.2569	1.7338	5.5059	26.6686	88.6510
		2010	1.0663	0.8176	5.4983	12.7586	88.6053
		2020	0.9523	0.8460	6.2854	11.3034	77.7368
	Wetland	2000	0.0150	59.5367	3998.1375	6.4252	10.0239
		2010	0.0073	28.0698	3876.9109	3.1660	10.0201
		2020	0.0245	26.8231	1106.5735	3.2972	10.6496
	Water	2000	0.1339	0.4646	27.1570	2.5533	15.7177
		2010	0.0636	0.2190	26.9675	1.2253	15.7296
		2020	0.0590	0.4581	38.6558	1.3597	15.1372
	Impervious	2000	0.0309	0.0974	26.2800	0.8348	11.0468
		2010	0.0146	0.0459	26.2541	0.4065	11.0360
		2020	0.0358	0.0760	13.5289	0.5788	13.9585
	Barren	2000	0.7919	0.1617	2.6690	5.6902	45.9799
		2010	0.3751	0.0762	2.6552	2.7299	45.9978
		2020	0.3277	0.0766	2.2895	2.1325	41.3830
Outside	Cropland	2000	0.1368	15.8733	355.8554	17.2783	33.4166
		2010	0.0458	5.2816	353.7959	6.1416	33.4174
		2020	0.0473	5.3034	342.3769	5.5302	30.0774
	Forest land	2000	1.6153	0.3581	5.0064	18.5621	83.2265
		2010	0.5380	0.1191	5.0007	6.2442	83.3100
		2020	0.4728	0.1075	4.6118	5.1588	76.4710
	Grassland	2000	2.5668	7.8802	9.4235	38.2741	99.7248
		2010	0.8531	2.6219	9.4354	12.9340	99.7737
		2020	0.6570	2.7245	12.5530	10.9510	83.4448
	Wetland	2000	0.0310	1.7465	216.8806	2.7310	14.7411
		2010	0.0102	0.5812	219.6958	1.0069	14.7423
		2020	0.0123	0.5469	175.4301	0.8924	13.3092
	Water	2000	0.1611	2.0709	36.7673	3.0695	16.7073
		2010	0.0537	0.6892	36.6798	1.0706	16.6878
		2020	0.0464	0.8132	45.2682	0.9491	14.3188
	Impervious	2000	0.0874	0.0930	30.3417	2.4834	19.8907
		2010	0.0291	0.0309	30.3381	0.8544	19.8925
		2020	0.0575	0.0941	24.0976	1.4573	27.0654
	Barren	2000	1.4402	0.2621	2.6158	10.6079	69.6651
		2010	0.4787	0.0869	2.6213	3.5640	69.6529
		2020	0.3985	0.0892	2.5273	2.8682	62.5410

and

Table 6. Landscape index of patch metrics MdfNNR in 2000, 2010, and 2020.

Inside and Outside	Type	Year	PDP	LPIP	MPSP	EDP	LSIP
Inside	Cropland	2000	0.1268	6.4395	94.1133	9.3318	10.8754
		2010	0.0530	3.2146	112.9304	5.0644	10.7794
		2020	0.0646	3.2608	92.4718	5.3503	11.3735
	Forest land	2000	0.0986	83.5368	887.4257	8.8498	4.7420
		2010	0.0461	40.6976	925.3890	6.8028	5.4300
		2020	0.0438	40.4178	974.7805	5.9945	4.7731
	Grassland	2000	0.1080	0.0051	0.3013	0.2706	5.3333
		2010	0.2764	0.0369	0.6330	1.0557	12.9492
		2020	0.0415	0.0037	0.4650	0.1384	5.0000
	Water	2000	0.0470	0.0275	1.6380	0.2818	3.7407
		2010	0.0230	0.0137	1.6470	0.1396	3.6071
		2020	0.0161	0.0129	2.5586	0.1315	3.2759
	Impervious	2000	0.0423	0.1420	9.8000	0.6341	4.3651
		2010	0.0276	0.0696	7.5750	0.4076	4.6094
		2020	0.0277	0.1036	12.7875	0.5895	5.1325
Outside	Cropland	2000	0.2160	7.0758	109.2690	16.2769	14.3367
		2010	0.0640	2.5252	123.8302	5.7580	13.6865
		2020	0.0695	2.3748	112.1436	5.9110	14.1780
	Forest land	2000	0.2246	67.9482	330.2256	14.276	9.1602
		2010	0.0807	21.0858	293.0416	6.8114	9.4241
		2020	0.0834	20.9960	282.2940	6.6256	9.1809
	Grassland	2000	0.1210	0.0066	0.3696	0.2890	5.4545
		2010	0.1363	0.0091	0.4711	0.4078	10.6522
		2020	0.0931	0.0039	0.3775	0.2326	8.2059
	Water	2000	0.1037	0.0350	1.1663	0.5119	5.5278
		2010	0.0348	0.0114	1.2168	0.1878	6.0811
		2020	0.0361	0.0113	1.2669	0.1943	5.9744
	Impervious	2000	0.0734	0.4949	27.9900	2.4104	6.9932
		2010	0.0376	0.0849	14.0800	0.7437	6.8538
		2020	0.0334	0.2850	22.7550	0.9180	7.0577

. There was a negligible difference in the average values of PD_P inside and outside ZlNNR. Notably, the LPI_P and MPS_P of inside wetland patches were much larger than those of outside wetland patches. Over the past 21 years, the change extents of the three indexes were larger outside wetland patches than those of the inside wetland. Conversely, the average values of LPI_P and MPS_P outside were generally more dwarf than those inside, especially for forest land types, and the average value of PD_P outside MdfNNR was also generally larger than that inside. The change extent of PD_P and LPI_P values was generally greater outside than inside in MdfNNR, whereas the change extent of MPS_P value shows the opposite trend, with the extent of change in outside indexes being smaller than that of the inside. Whether in ZlNNR or MdfNNR, in terms of patch shape, the average values of ED_P and LSI_P are generally higher outside than inside, and the variation ranges of the two indexes were the same.

The landscape level of the landscape pattern indexes for each study area are presented in

Table 7. Landscape index of landscape metrics ZlNNR in 2000, 2010, and 2020.

Inside and Outside	Year	PDL	LPIL	MPSL	EDL	LSIL	CONTAGL	SHDIL	SHEIL
Inside	2000	4.6491	59.5367	21.5097	33.2357	40.0213	61.5834	1.2795	0.6575
	2010	2.1988	29.8795	45.4803	16.3355	28.4977	66.2232	1.2948	0.6226
	2020	1.9791	36.0420	50.5292	14.9965	26.2452	65.8515	1.3148	0.6323
Outside	2000	6.0385	15.8733	16.5603	46.5032	61.9886	54.4824	1.4657	0.7532
	2010	2.0089	27.8439	49.7789	16.3367	36.7712	70.2768	1.1238	0.5404
	2020	1.6924	27.8517	59.0880	14.3315	32.3780	70.4998	1.1236	0.5403

and

Table 8. Landscape index of landscape metrics MdfNNR in 2000, 2010, and 2020.

Inside and Outside	Year	PDL	LPIL	MPSL	EDL	LSIL	CONTAGL	SHDIL	SHEIL
Inside	2000	0.4228	83.5368	236.5390	9.6841	5.1603	85.3014	0.4011	0.2492
	2010	0.4284	50.9712	233.4431	7.8288	5.0759	73.0289	0.9025	0.5037
	2020	0.2006	40.4178	498.3910	7.1997	4.7459	73.2023	0.9002	0.5024
Outside	2000	0.7387	67.9482	135.3763	16.8821	9.7542	75.9286	0.6538	0.4062
	2010	0.3589	29.6175	278.6023	8.3511	6.5954	74.6234	0.8411	0.4694
	2020	0.3210	36.3360	311.5052	8.3445	6.5939	74.4851	0.8460	0.4722

. Regarding the landscape characteristics, the change extent of the PD_L and MPS_L values was larger outside ZlNNR than inside, with the LPI_L value being greater inside ZlNNR than outside. The change extent of the PD_L and LPI_L values outside MdfNNR was larger than inside, whereas the three-year mean of LPI_L and MPS_L values was smaller outside compared to inside, with PD_L showing the opposite trend. As for changes in landscape shape, the four indexes demonstrate greater fluctuation outside the ZlNNR than inside, whereas the three-year mean of ED_L and LSI_L values was higher outside than inside, and CONTAG_L was similar between inside and outside. Additionally, the change extent of the ED_L and LSI_L values was smaller inside MdfNNR than outside, and the three-year mean of ED_L and LSI_L values was higher outside than inside. Regarding landscape diversity, the change extent of the SHDI_L and SHEI_L values was smaller inside ZlNNR than outside, and the average values of both indexes were larger inside ZlNNR than outside, but vice versa for MdfNNR.

Comparing the patch indexes between the two PAs, we found that from 2000–2020, the PD_P, LPI_P, and ED_P values increased in MdfNNR, with water bodies increasing the most (65.74%, 53.09%, and 53.34%, respectively). The MPS_P value only increased for cropland, but the increase was small, at just 1.74%, whereas other patch types decreased. Water bodies decreased the most (56.20%). The LSI_P value change was negligible. In ZlNNR, both LPI_P and ED_P value increased, with the LPI_P value of cropland increasing the most (54.95%) and the ED_P value of barren land increasing the most (62.52%). The PD_P value of wetlands decreased the most (63.33%), whereas MPS_P increased the most for wetlands (72.32%). The LSI_P value change was also negligible. Overall, wetland and impervious surface patches became more regular in shape in ZlNNR, whereas patch types became increasingly fragmented. In MdfNNR, all patch types became increasingly fragmented, especially forest, grassland, and water body patches. For landscape index, the MPS_L, CONTAG_L, and SHDI_L values of ZlNNR increased by 134.91%, 6.93%, and 2.76% respectively, whereas the rest of the index values decreased. Notably, the PD_L value decreased by 57.43%. For MdfNNR, the MPS_L, SHDI_L, and SHEI_L values decreased by 110.7%, 124.43%, and 101.61%, respectively, whereas the other index values increased. Notably, the PD_L value increased by 52.55%. This indicates that the MdfNNR landscape became increasingly fragmented with weaker connectivity, reduced diversity, and irregularity. In contrast, the ZlNNR landscape showed smaller distribution evenness changes, enhanced connectivity, increased aggregation, and more regular shapes.

The landscape-level landscape pattern indexes for both NNRs are visualized in Figure 5 and Figure 6. Figure 5 reveals that the CONTAG and MPS indexes in ZlNNR exhibit a noticeable increasing trend, whereas the other indexes show no significant changes. The gap between the experimental zone and the core zone is prominent, with the extreme values of each landscape index observed in these areas. For instance, the ED and LSI of the core zone are lower, whereas the LPI is larger. Additionally, the PDs of the core zone and experimental zone are both lower, suggesting that the closer to the edge of ZlNNR, the higher the fragmentation of the landscape. However, SHDI and SHEI follow the pattern of experimental zone > buffer zone > core zone, indicating greater landscape diversity in the experimental zone. Figure 6 shows that the indexes inside MdfNNR have changed significantly in the southeast, displaying an increasing trend from 2000 to 2010 and a decreasing trend from 2010 to 2020. Among them, CONTAG, ED, LSI, SHDI, and SHEI all exhibit a pattern of experimental zone > buffer zone > core zone, whereas LPI and MPS follow the pattern of core zone > buffer zone > experimental zone. This suggests greater landscape diversity in the experimental zone and a less fragmented and more regular landscape in the core zone. Comparing the landscape indexes of the two NNRs reveals that MdfNNR experiences greater landscape fragmentation, whereas ZlNNR exhibits a more regular landscape shape.

3.3. NDVI

3.3.1. The Temporal Variation Characteristics

The interannual variation of the NDVI, as depicted in Figure 7, demonstrates that the annual average NDVI for each study area oscillated within a specific range and exhibited consistent temporal trends. The average annual NDVI inside and outside MdfNNR and ZlNNR ranged from 0.436–0.712, 0.439–0.690, 0.221–0.452, and 0.183–0.370, respectively. The annual average NDVI has been higher inside ZlNNR than outside over a period of nineteen years. The difference in NDVI between the inside and outside of MdfNNR is relatively small, but, overall, the inside has slightly higher NDVI values. The annual average NDVI of the MdfNNR was higher than that of the ZlNNR, primarily due to the predominance of forest land inside the MdfNNR. More specifically, the annual average NDVI of the MdfNNR experienced a steady decline from 2000 to 2003, although the rate of decline diminished each year. In 2010, the MdfNNR’s annual average NDVI peaked. From 2018 to 2020, the annual average NDVI of the MdfNNR consistently increased, with a more rapid increase observed inside than outside MdfNNR. The annual average NDVI of the ZlNNR consistently increased from 2003 to 2007, with a more rapid increase observed inside than outside ZlNNR. In 2007, the ZlNNR’s annual average NDVI reached its peak. Overall, both NNRs exhibited an upward trend in their annual average NDVI.

3.3.2. The Spatial Variation Characteristics

The distribution of improved, stable, and degraded vegetation inside ZlNNR accounted for 96.93%, 1.11%, and 1.96%, respectively (

Table 9. Trend of NDVI of ZlNNR, and MdfNNR.

β	Z	NDVI Trend	Area Percentage
			ZlNNR		MdfNNR
			Inside	Outside	Inside	Outside
≥0.0005	≥1.96	Significantly improved	64.17%	62.40%	39.71%	38.88%
≥0.0005	−1.96–1.96	Slightly improved	32.76%	30.91%	48.64%	44.08%
−0.0005–0.0005	−1.96–1.96	Stable	1.11%	2.10%	4.84%	5.55%
<−0.0005	−1.96–1.96	Slightly degraded	1.67%	4.04%	6.74%	11.00%
<−0.0005	<−1.96	Severely degraded	0.28%	0.55%	0.08%	0.50%

). In contrast, outside ZlNNR, they are 93.31%, 2.10%, and 4.59%. The vegetation improvement area was 3.62% greater inside than outside, and the vegetation degradation area was 2.63% less inside than outside. As shown in Figure 8a, the area of vegetation improvement in the ZlNNR was far greater than the area of vegetation degradation inside and outside. The areas with improved vegetation were predominantly located in the central region of the NNR and the southeast outside the NNR. Areas with slight vegetation improvement were predominantly distributed in the northeast and southwest regions of the NNR and various regions outside the NNR. Areas with stable vegetation were mainly found in the southwest region. In contrast, areas with slight and significant vegetation degradation were primarily located near water bodies and impervious lands.

The areas with improved vegetation, stable vegetation, and degraded vegetation in MdfNNR account for 88.35%, 4.84%, and 6.82%, respectively (Table 9). Conversely, outside MdfNNR, they are 82.96%, 5.55%, and 11.49%. The improved vegetation area was 5.39% greater inside than outside, and the degraded vegetation area was 4.68% less inside than outside. As shown in Figure 8b, the area of vegetation improvement was far greater than the area of vegetation degradation, whether inside MdfNNR or outside. Areas of notable vegetation improvement were predominantly situated in the southwest, whereas areas of slight improvement and stable vegetation were more dispersed. Regions with minor vegetation degradation were primarily found near roads, and areas with severe vegetation degradation were chiefly located near impervious land.

The time series analysis comparing the two PAs shows that ZlNNR has a higher percentage of areas with improving NDVI trends than MdfNNR. Specifically, the percentage is 8.59% higher in ZlNNR. For significantly improving trends, ZlNNR has 24.46% more area than MdfNNR. For slightly improving trends, ZlNNR has 15.88% less area than MdfNNR. On the other hand, the percentage of areas with degrading NDVI trends is 4.86% lower in ZlNNR than in MdfNNR. For severely degrading trends, ZlNNR has 0.06% more area than MdfNNR. For slightly degrading trends, ZlNNR has 6.96% less area than MdfNNR. However, the percentage of areas with stable NDVI is 3.45% lower in ZlNNR than in MdfNNR, indicating relatively stable NDVI in MdfNNR. In summary, the overall NDVI condition is better in ZlNNR compared to MdfNNR.

3.4. Comprehensive Analysis of Conservation Effectiveness

The weights assigned to each indicator in the study area were determined using the coefficient of variation method, as detailed in

Table 10. Weights of each evaluation index in ZlNNR and MdfNNR.

Primary Evaluation Indicators	Weights of Primary Evaluation Indicators		Secondary Evaluation Indicators	Weights of Secondary Evaluation Indicators
	ZlNNR	MdfNNR		ZlNNR	MdfNNR
Land use type	0.1691	0.4972	Cropland	0.0682	0.0152
			Forest land	0.1336	0.0050
			Grassland	0.0225	0.7307
			Wetland	0.0370
			Water	0.1654	0.0532
			Impervious	0.3462	0.1959
			Barren	0.2270
Landscape pattern index	0.6266	0.3549	PD	0.2406	0.1955
			LPI	0.1366	0.0759
			MPS	0.1808	0.1401
			ED	0.2257	0.1836
			LSI	0.1147	0.0958
			CONTAG	0.0321	0.0370
			SHDI	0.0287	0.0418
			SHEI	0.0388	0.0527
Remote sensing index	0.2043	0.1479	NDVI	1.0000	1.0000

. Figure 9 presents the calculated results of conservation effectiveness for each study area. The proportions of the area corresponding to the three levels of conservation effectiveness in each study area were also computed, as depicted in Figure 10.

For inside ZlNNR, the area characterized by HE gradually increased, whereas the area of ME progressively decreased. The area of LE initially shrank and then increased. The relative proportions of the three conservation effectiveness levels in each period followed the pattern ME > HE > LE. Conversely, outside ZlNNR, the area of LE initially increased and then decreased, whereas the areas of both ME and HE first decreased and then increased. The relative proportions of the three conservation effectiveness levels in each period followed the pattern ME > LE > HE. Moreover, the inside had a higher area proportion of HE and ME than the outside each period, whereas LE was the opposite.

For inside MdfNNR, the area of LE continues to decrease, whereas the area of ME continues to increase, and the area of HE initially increased and then decreased. The relationship between the areas occupied by the three conservation effectiveness levels in each period was ME > LE > HE. Conversely, for outside MdfNNR, the area occupied by LE initially decreased and then increased, whereas the areas occupied by both ME and HE increased and then decreased. The relationship between the areas occupied by the three conservation effectiveness levels in each period was ME > HE > LE. LE areas inside MdfNNR were less than those outside each period, and the ME areas were higher than those outside each period. For HE areas, the difference between the inside and outside of the PAs was minimal.

Both ZlNNR and MdfNNR have served their purpose of protecting vegetation and landscape to a certain extent. As expected, the core zones offer the highest conservation effectiveness, followed by buffer and experimental zones. The areas of HE in ZlNNR were more significant than those in MdfNNR in each period, and the areas of ME and LE in ZlNNR were smaller than those in MdfNNR in each period. Therefore, the conservation effectiveness of ZlNNR was superior to that of MdfNNR.

4. Discussions

Our research contributes to the assessment of conservation effectiveness, to quantify the differences in conservation effectiveness across various types of PAs. We found that the conservation effectiveness of ZlNNR surpassed that of MdfNNR. This discrepancy could be attributed to a multitude of factors, including the size and location of the PAs, the level of economic development in the surrounding areas, and local population density. Future studies could delve deeper into these factors, providing valuable insights that could inform the creation of more effective conservation strategies and policies.

4.1. Analysis of Changes in Conservation Effectiveness and Their Driving Forces

Our research findings reveal that the degree of land use change in ZlNNR was relatively less than that observed in MdfNNR. This disparity can be attributed to the unique geographical locations and climatic characteristics of the two NNRs. ZlNNR, situated in a water-abundant environment, has limited human activity, and the implementation of green policies has further mitigated damage to wetlands. Consequently, changes in land use are primarily driven by climatic factors. Conversely, MdfNNR is an urbanized area with extensive transport networks and a high rate of economic growth. Here, societal and economic factors primarily shape land use transformations. Numerous studies have suggested that intense human intervention often leads to significant land degradation [21,71,84].

Our research indicates that MdfNNR experienced a higher degree of landscape fragmentation, largely due to its extensive transport infrastructure and dense network of roads [39]. This fragmentation of the internal landscape hinders species movement, thereby obstructing the effective preservation of biodiversity. Conversely, ZlNNR also saw a surge in road construction during the 1990s, primarily due to the development of highways and railways, which influenced landscape fragmentation [98,99]. However, the linear characteristics of highways and railways facilitate ecotourism growth and make it easier to mitigate their adverse impacts on the ecosystem. Moreover, ZlNNR’s diverse land use types and robust ecological features enhance its resilience to landscape fragmentation. Since the implementation of various protective measures in 2005 [61], the degree of landscape fragmentation in ZlNNR has significantly decreased, and the landscape diversity and connectivity have improved. This improvement can be attributed to the enhanced interconnectivity between the core area, buffer zone, and external environment, achieved through the construction and restoration of ecological corridors, the building of ecological bridges, and the implementation of slope greening measures. These initiatives have effectively mitigated the negative impacts of previous road construction, thereby establishing a solid foundation for the sustainable conservation of regional biodiversity.

Our NDVI study results indicate that the percentage of areas exhibiting an improving trend is higher for ZlNNR than for MdfNNR, whereas the percentage of areas exhibiting a degrading trend is lower for ZlNNR compared to MdfNNR. Overall, the NDVI conservation effectiveness of MdfNNR was inferior to that of ZlNNR. This slightly lower conservation performance of MdfNNR could be attributed to the greater impact of human activities and a prolonged drought from 1999 to 2002 [100,101]. Although ZlNNR also experienced this concurrent drought, a long-term water replenishment mechanism implemented since 2005 contributed to the earlier recovery of the NDVI [60,61].

The core area of NNRs typically serves a pivotal role in safeguarding biodiversity, ecosystems, and other related aspects [12]. This protective function facilitates the congregation and survival of various species, thereby enhancing the effectiveness of biodiversity and ecosystem conservation within the area. In our study, the core area, experiencing minimal external interference, presented a relatively stable and intact ecological environment, leading to a concentration of HE regions. The overall conservation effectiveness of ZlNNR surpassed that of MdfNNR, a finding that is consistent with our observations on land use, landscape patterns, and NDVI in the two NNRs. This superior performance can be attributed to ZlNNR’s robust natural attributes, more comprehensive ecological environment, and higher landscape connectivity. In contrast, MdfNNR faced significant ecological and environmental disturbances due to human activities, which inevitably curtailed its conservation effectiveness, leading to differential impacts between the two regions. To enhance the conservation effectiveness of MdfNNR, it is imperative to minimize the damage inflicted on its ecological environment by human activities [102]. The implementation of practical measures such as controlling unreasonable development and usage, restoring ecological corridors, and improving landscape connectivity could significantly bolster the ecological security and biodiversity conservation effectiveness of MdfNNR. Climate change, as a natural factor, can also exert substantial impacts on biodiversity and ecosystems in different regions [103]. For the effective and long-term protection of biodiversity in PAs, it is necessary to mitigate both the negative impacts of human activities and natural factors concurrently.

Holistic, multi-scale, and interconnected contemplation, strategizing, and management are important, yet still relatively lacking, aspects of current PA management [32]. The results of this study emphasize the necessity of adaptive management tailored to different ecosystem types. The studied wetland PA and forest PA require different management measures. Wetlands have higher sensitivity and require strict protection measures and proposing ecological restoration projects [33], while also needing more proactive restoration measures to improve degraded forests. Comparative results show management plans must consider the unique ecological characteristics, species compositions, and ecological functions of each ecosystem. The inherent complexities and incompatibilities between systems imply specialized PAs for particular ecosystems may allow more customized and effective protection. Establishing ecosystem PAs, safeguarding the unique attributes of habitats such as wetlands, forests, and grasslands, maintains their integrity. It enables tailored monitoring, restoration projects, and policy fine-tuning to local ecology. However, we still know little about future habitat suitability in PAs for colonization by alien species. Under climate change, predicting and assessing the receptivity of different ecosystems to alien organisms provides a basis for vulnerability assessments and adaptive strategy development. Additionally, concrete assessments of native species in PAs of varying ecosystems are needed to identify highly vulnerable species within different PA types, enabling timely targeted translocation [33]. Furthermore, our research found that both core and buffer zones of both PAs also experienced varying degrees of anthropogenic disturbance, although impacts were smaller in core zones. This illustrated that the core and buffer zones of different types of PAs had not been completely free from disturbance. In the future, the administration of PAs and higher-level departments will need to invest more manpower and material resources, increase patrol efforts, ensure the safety of protected objects, and achieve the purpose of on-site protection. Therefore, while establishing multifunctional PAs, developing specialized PAs based on specific ecosystems remains an important strategy for comprehensive and adaptive protection management.

4.2. Evaluation on Research Method Selection and Application

Land use serves as a significant indicator of the interaction between human activities and ecosystems [104], with extensive cropland areas serving as a stark testament to human influence. The transfer matrix of land use reveals the trends and principles of land use changes by examining the exchanges among various land use categories at different time points. The land use dynamic index is a comprehensive metric that encapsulates the direction, pace, and intensity of land use changes, offering superior practicality and advantages [70]. These two methods have been widely utilized in land use studies within PAs. Yan [71] and Mu [21] employed the land use dynamic index to monitor and evaluate the spatiotemporal changes of wetlands, revealing a lack of effectiveness in wetland conservation, a finding that contrasts with our results. Using a transition matrix from a different land use classification system, which subdivides forests into three categories, Bonilla observed increases in forest cover of 5% and 3% inside and outside PAs, respectively [105]. In contrast to previous research [73,105,106], we utilized the aforementioned two methods to conduct a comprehensive and accurate analysis of land use variation. In our study, forests both inside and outside the MdfNNR exhibited a decreasing trend. However, the declining trend inside the reserve was less pronounced than outside, which, to some extent, underscores the effectiveness of conservation efforts.

We used the GlobeLand30 data, which has emerged as a widely accepted resource for analyzing land use patterns [8,67], to investigate land use changes. In contrast, the previous research [61] conducted relied on the China Multi-temporal Land Use Remote Sensing Monitoring Dataset of the Resource and Environmental Science and Data Center, whereas other scholars opted for Landsat and Sentinel images for land use classification [107,108]. Though each of these data sources have their unique advantages, they vary in their estimations of the extent of different land use types, primarily due to differences in data quality and classification methods. Our investigation indicates that the degree of land use changes in ZlNNR was minimal between 2000 and 2010, with more noticeable variations emerging between 2010 and 2020. This sudden shift between these two periods differs from the findings of previous research [61,107,108]. This discrepancy is primarily due to inconsistencies in small-scale land use classification across different datasets. Therefore, we propose that future research should more effectively consider the scale and location of the study area when selecting data sources. This approach will enable a more accurate analysis of land use patterns.

In contrast to previous studies, the landscape index chosen for this research is distinguished by its comprehensive scope and systematic approach. Earlier studies, such as Mu [21], evaluated the conservation effectiveness of wetland conservation areas based solely on landscape fragmentation, whereas Lu [38] relied exclusively on the Landscape Development Intensity Index for their assessments. In this study, we utilized ten landscape indexes spanning five different aspects, including patch characteristics, patch shape, landscape characteristics, landscape shape, and landscape diversity. This approach allowed us to create a comprehensive quantitative characterization of spatial and shape characteristics across patches and landscapes. This methodology enabled us to accurately analyze the landscape patterns of PAs within a scientific and comprehensive framework. Moreover, we also adopted a visual analysis of landscape indexes, which offers superior analytical efficiency compared to previous studies [109,110,111,112], benefiting both researchers and readers. We deviated from the traditional moving window method [76,113] and instead embraced a combination of grid analysis and Kriging interpolation to provide accurate landscape boundary analyses. This approach aligns more closely with the latest trends in research.

The combination of the Theil–Sen median trend analysis and the Mann–Kendall test has proven to be a reliable approach for time series analysis. This method offers several advantages over univariate linear regression analysis, including the lack of requirement for data to follow a normal distribution, reduced sensitivity to outliers and extreme values, the effective detection of sequence trends and change points, and lower susceptibility to errors. Furthermore, this method is relatively straightforward to calculate, enhancing its practical applicability [84,114]. In their research on vegetation coverage in the Yellow River Basin, Jiang [84] also employed the Hurst index. We attempted to use this method as well, but found that the Hurst index is more suitable for studies involving large-scale research areas.

To provide a comprehensive and accurate assessment of ecological environment changes within the PA, this study utilized the NDVI as the primary remote sensing evaluation indicator. NDVI is recognized for its ability to quantitatively characterize the growth status and distribution changes of vegetation. Its widespread application in the scientific community also enhances the comparability and credibility of research results [34,45,84]. Though the Fraction of Vegetation Cover (FVC) and Remote Sensing Ecological Index (RSEI) were initially tested for their applicability in detecting ecological environment changes, NDVI proved to be more effective. This was due to its superior ability to filter out external environmental interference and accurately extract vegetation information. Conversely, FVC and RSEI were found to be more susceptible to significant errors caused by other environmental factors, making it difficult to extract useful data related to the research objectives. Consequently, NDVI was chosen as the sole remote sensing indicator to detect spatiotemporal changes in the ecological environment of the PA and evaluate the effectiveness of its protection. Our research findings indicated that a higher percentage of vegetation exhibited an improvement trend inside the PA than outside, suggesting that the PA has a certain degree of effectiveness. This conclusion aligns with previous studies, which demonstrated that nature reserves contribute to the improvement of sustainable and stable geographical spaces [34].

In addition, this study leveraged the GEE platform, which has gained prominence in the realm of remote sensing image processing, thanks to the rapid advancement of computer technology. The datasets within the platform require no preprocessing, thereby enhancing research efficiency [47,115]. Our study highlights the efficacy of remote sensing as a tool for evaluating the conservation effectiveness of PAs. Interestingly, the results derived from remote sensing were found to be almost identical to those obtained through traditional ground-based methods. This congruence underscores the validity and reliability of remote sensing as a method for conservation assessment. However, the advantage of remote sensing over ground-based methods lies in its timeliness and convenience. Remote sensing allows for real-time data collection and analysis, which can significantly enhance the effectiveness of conservation efforts. Furthermore, it provides a convenient and less intrusive means of monitoring PAs, reducing the need for physical intrusion that could potentially disrupt the ecosystems under protection.

In this study, we developed a generalizable framework for assessing the conservation effectiveness of PAs by evaluating three key factors: land use, landscape pattern, and NDVI. These three factors apply to various types of PAs and can be visualized to improve the clarity and intuitiveness of the analysis. The framework can effectively compare the conservation effectiveness of different types of PAs, such as forest and wetland types. We objectively determined the weight of each factor using the coefficient of variation method and constructed a conservation effectiveness assessment model using weighted summation, making the assessment results more scientific and convincing. This highly operable model provides a standardized analytical framework and assessment method for evaluating the conservation effectiveness of different types of PAs. The assessment model established in this study possesses expansive applicability for appraising conservation effectiveness in extant PAs of all typologies, and proffers references for prospective PA planning, hence possessing salient practical value.

4.3. Limitations and Future Trends

NDVI can become saturated in areas of high-density green biomass, making it challenging to detect changes in regions with 100% vegetation coverage. Additionally, NDVI may not accurately reflect critical vegetation characteristics during seasonal peaks [48,116,117]. Furthermore, within the ZlNNR, when wetlands transitioned to grasslands, the NDVI may exhibit an increasing trend, which is a false indication of effectiveness. This is a limitation of NDVI as a single assessing index. Nevertheless, we used a multi-factor collaborative framework to diminish or compensate for such deficiencies in this study. Other factors capturing changes in land use type and habitat fragmentation could reveal when apparent “improvements” in vegetation greenness may counter-intuitively signal deteriorating conditions in PAs. Our study underscores the importance of examining conservation outcomes through diverse lenses to obtain a comprehensive perspective. Therefore, in our future work, it will be necessary to corroborate the results obtained using NDVI with other data sources, such as high-resolution imagery and field survey data. To achieve more precise NDVI results and research conclusions, it is crucial to select appropriate observation periods to minimize the impact of seasonal changes in spectral signals on a temporal scale. High-resolution data can enhance the discernibility of various vegetation and ecological features within pixels at the spatial scale, partially compensating for the limitations of NDVI data in areas of high vegetation coverage. However, achieving high temporal and spatial resolutions simultaneously presents a challenge. Moreover, existing spectral indicators have limitations in accurately inverting vegetation feature parameters [117]. Therefore, we strongly recommend that future research employs advanced remote sensing methods, such as microwave remote sensing, to monitor and evaluate the vegetation dynamics in PAs more accurately and comprehensively.

Apart from the types of PAs analyzed in this study, there are other important ecosystems such as wildlife and grasslands that merit attention. However, our current conservation effectiveness evaluation system does not encompass all types of PAs, presenting a challenge for future research. We aim to direct future research efforts towards addressing this limitation. To develop a comprehensive conservation effectiveness evaluation system, it is necessary to conduct extensive and in-depth studies on methodologies, indicators, and management techniques. This endeavor requires not only building upon existing research findings, but also collaborating with experts from various fields to establish a more comprehensive and scientifically rigorous evaluation framework and research paradigm. We are committed to achieving these objectives to provide increasingly accurate, comprehensive, and advanced scientific support for PA management. Furthermore, it is crucial to pay close attention to minimizing the spillover effects of PAs [118,119,120], an issue that will be a key focus of our future research endeavors.

5. Conclusions

The quantitative evaluation model established in this study utilizes technologies including land use dynamic index, transition matrix, landscape pattern analysis, and NDVI time series to accurately measure changes in ecological quality and conservation effectiveness across PAs with varying ecosystems. By comparing reserves protecting different ecosystems, this research validates the methodological approach of using multiple quantitative techniques to comprehensively evaluate indicators of conservation effectiveness. The novel ecological and technological methodologies introduced establish a scientific system to assess conservation effectiveness. We found that significant disparities were observed in land use changes between the two PAs. ZlNNR experienced a substantial impact from wetland reduction, whereas MdfNNR was predominantly affected by forest reduction. The landscape patterns inside and outside PAs exhibited notable differences, characterized by higher internal landscape quality and reduced fragmentation. Both reserves exhibited smaller internal fluctuations in landscape pattern indexes compared to external fluctuations, indicating relative ecological stability. However, ZlNNR demonstrated higher landscape quality and stability compared to MdfNNR. Additionally, the percentage area of NDVI in ZlNNR exhibited a more pronounced increasing trend versus MdfNNR. Our evaluation method holds substantial theoretical and practical value, enriching the foundation and techniques for evaluating nature reserve efficacy and offering innovation for ecology and environmental management fields.