1. Introduction
The establishment of protected areas is a means of maintaining biodiversity while assuring the sustainable provision of ecosystem services [
1,
2,
3]. In the tenth meeting of the Conference of the Parties regarding the Convention on Biological Diversity, one of the most critical Aichi Targets was set: at least 17% of terrestrial land and 10% of waters should be effectively protected by 2020. In the fifth edition of the Global Biodiversity Outlook, it was concluded that this target had been partially achieved [
4]. An updated target, increasing these protected areas to 30% by 2030, has been proposed in the first drafts of the fifteenth meeting of the Conference of the Parties [
5] and could be recommended as a new deal for humans and nature. However, the debate over the necessary proportion of protected areas seems likely to persist because of the complexity of the recommendations provided by studies with different spatiotemporal patterns and selection methods and focus on different ecological processes. A recent review indicated that the protected area coverage called for by previous studies has a large range (30–70%) at a global level [
6].
Management methods and the corresponding intensity of land-use regulation is another important point in setting up protected areas, in addition to their spatial extent. Although hierarchical and zoning-based protected area management has reached international consensus, owing to the complexity of social development and ecosystems, there are almost no universal solutions applicable to all countries or regions. In some developed countries, such as the United States and New Zealand, national parks have achieved remarkable results. However, this approach may not be appropriate in some developing countries, owing to their large rural populations which rely on natural resources and have a low tolerance for the large predators that play essential roles in many ecosystems [
7]. Poor management has led to poaching and logging in some protected areas, turning them into “paper parks” e.g., [
8,
9], a severe issue that requires awareness and action [
10].
Despite these difficulties, some countries and regions are still attempting to develop new methods of establishing and governing protected areas. In recent years, China has been implementing a national ecological conservation policy: the Ecological Conservation Redlines (Eco-redline) policy, which refers to the use of red lines to demarcate areas that have important ecological functions within the scope of ecological space and which must therefore be strictly protected [
11]. By assessing the ecosystem services and ecological sensitivity of the target areas, areas with important ecological functions and ecological sensitivity are identified and then superimposed and calibrated with existing protected areas (e.g., nature reserves and national parks) to form Eco-redline areas (ERAs).
As one of the most important ecological policies in China, the Eco-redline policy will have a lasting influence on many aspects of social and economic development. Most government departments fall under the scope of the policy, whether they make development plans or environmental protection plans. Land use and land cover (LULC) result from human–nature interactions, and changes may be beneficial or detrimental to human beings, profoundly affecting human well-being and welfare [
12]. Therefore, evaluating the impact of the Eco-redline policy on LULC is crucial to the future implementation and optimization of this policy.
In this regard, scenario analysis is considered an effective method of assessing the medium- and long-term impacts of a policy [
13]. Scenario analyses in protected area studies often consider alternative policy interventions for protected areas and the drivers of change (e.g., climate change, economic development, and population growth) that influence the operation of the protected areas. For example, using data on the area and percentage of protected areas in countries worldwide from 1950–2005, McDonald and Boucher [
14] modeled future projections of protected areas in 2030 to examine the effectiveness of two contrasting strategies: strict conservation and multiple-use (in which resource extraction is partly permitted). Based on LULC change between 1990 and 2001, Martinuzzi, et al. [
15] quantified areas of urban landscapes, croplands, and natural vegetation around protected areas in the United States under business as usual, forest incentive, high crop demand, and urban containment scenarios. Velazco, et al. [
16] modeled plant species losses resulting from the continuous emission of greenhouse gases in Bolivia, Brazil, and Paraguay in 2050 and 2080, and found that the current protected area network is not sufficient to safeguard the most valuable Corrado plant species, even in the most optimistic scenario.
There are also several studies on the demarcation and management of Eco-redlines, some of which used scenario analysis. For instance, taking Shanghai as an example, Bai, et al. [
17] analyzed multiple land-use scenarios with different policy interventions to explore their impacts on LULC and ecosystem services and assess their implications for the effective implementation of the Eco-redline policy. Using the CLUE-S model, Jia, et al. [
18] explored the impacts of the Eco-redline policy on spatiotemporal land-use changes in Beijing. They concluded that the Eco-redline policy could improve the spatial integrity and connectivity of ecological functions. Ju, et al. [
19] projected the urban expansion of the Beijing–Tianjin–Hebei megaregion in China by 2030 with and without the Eco-redline policy to show the effects of the policy on runoff.
However, there is not much discussion on how to further optimize the Eco-redline policy. Many studies compared scenarios in which the Eco-redline policy either exists or not, to illustrate the positive effects or deficiencies of the policy but did not provide specific adjustment suggestions. In addition, restricted by practical conditions such as financial support especially in developing countries or regions, there are usually trade-offs between different approaches. Therefore, which approach should be prioritized needs to be considered in policymaking. About the study area, existing studies most often use the northern plains or the developed eastern regions of China as examples. These studies have paid little attention to China’s central and western regions, which are dominated by mountains and hills and remain relatively underdeveloped. Thus, we aimed to examine the effectiveness of alternative policy interventions and their implications on future LULC patterns, focusing on the capital of Chongqing, the only municipality under direct control of the central government in western China. We employed scenario analysis with different ERA spatial extents and management intensities, using LULC simulation and landscape indices to analyze how these interventions influence ERA effectiveness. According to the implementation outline of Eco-redline policy from the central government [
11], each province or city in China needs to regularly (e.g., every 5 or 10 years) evaluate the effectiveness of its ERAs and make appropriate adjustments if necessary. By comparing the impact of different measures on the landscape, we try to evaluate which approach is more critical to promoting the Eco-redline policy; thus, our results will provide theoretical support in subsequent policymaking.
2. Materials and Methods
2.1. Study Area
Chongqing is in southwest China, in the upper reaches of the Yangtze River (
Figure 1). The capital of Chongqing is its administrative and economic center. It covers a total area of 5470 km
2 and has nine districts under its jurisdiction, with a total population of nearly 8 million. The landforms of Chongqing are dominated by mountains and hills, of which mountains account for 76%, giving it the name “mountain city” [
20]. Since the most influential development strategy in China, named Reform and Opening Up in 1978, Chongqing has undergone great changes in social and economic aspects over four decades. From 2000 to 2018, almost every year, the annual GDP growth rate of Chongqing exceeded 10%. In 2019, the GDP of Chongqing was USD 342 billion, with a medium to high growth rate (6.3%), and GDP per capita just exceeded USD 10,000 [
21]. Economic development has brought rapid urbanization, with the urban population increasing from 35.6% in 2000 to 65.5% in 2018. According to the latest Eco-redline delimitation plan released by Chongqing Municipal People’s Government [
22], the ERAs in Chongqing capital cover 912 km
2, accounting for 17% of the total area. As there was no GIS data of the ERAs available, we produced vector data from a digital Eco-redlines map, showing the spatial distribution of the ERAs. The digitized ERAs covered a total area of 906 km
2, which is very close to the official datum (912 km
2).
2.2. Framework of Analysis
Scenario analysis is the combined application of models and scenarios [
23]. With the development of GIS (geographic information system) technology, it becomes more and more convenient to obtain and analyze LULC data or use it for modeling. In the current study, we used the land change modeler (LCM) on TerrSet (version: 18.31) for land change simulation, which allows users to consider different policy interventions such as zoning regulations, land development plans, and road developments. We used 2050 as the time horizon for future projection and LULC maps from 2000 and 2010 for between different LULC classes, based on Chongqing’s socio-economic development. In addition, an LULC map from 2015 was used for the verification of the land change simulation.
Figure 2 demonstrates the framework of the analysis. Generally, it can be divided into the following steps: the construction of the land change model, the setting and conversions of the scenarios, the generation of future LULC images, and the calculation and statistical test of the landscape indices of the LULC images. First, cloud-free historical satellite images of 2000, 2010, and 2015 were acquired. The satellite images were then used to conduct LULC classification. The classified 2000 and 2010 LULC images were then used for land change modeling. A predicted 2015 LULC was created using the model and compared with the real 2015 LULC image. The model was then readjusted multiple times until satisfactory accuracy was obtained. Once the model was ready, the second step was the setting and converting of scenarios. Six scenarios based on different assumptions which incorporated the change of area and management intensity were proposed. The interference layers corresponding to different scenarios were generated and loaded into the model to simulate the future LULC images in 2050. Next, to compare the differences between these LULCs under different scenarios, several key forest patch landscape indices were calculated for each scenario. Finally, the landscape indices were used for statistical tests to find the effectiveness of the change of area and management intensity of ERAs. The following sections will provide more detailed explanations.
2.3. Data Acquisition of Historical Satellite Images
Landsat surface reflectance data images provided by the United States Geological Survey were sourced using the Google Earth Engine (GEE) platform. The Landsat surface reflectance data has a spatial resolution of 30 m, a temporal resolution of 16 days, and has been atmospherically corrected. GEE provides a function to filter for cloud-free pixels (
https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LT05_C01_T1_SR accessed on 20 September 2021) by reading the description of the cloud cover in the metadata. In addition to this function, the IMAGE_QUALITY and CLOUD_COVER image properties were also used during pixel filtering.
To minimize phenological influences, the 2000, 2010, and 2015 images were preliminarily screened by month. After many attempts, we accepted that almost no image of a single month could cover the entire study area. In addition, cloudless images in summer and autumn generally covered no more than half the area. Therefore, we decided to use multiple images from adjacent periods for the subsequent land classification and then overlap them to obtain an LULC map of the whole study area. The following two criteria were applied in selecting images: first, the date range needed to be within three years of the target year; second, images from two months were used for each target year. Compared with the strategy of fusing images from different periods first and then classifying them, this strategy of first selecting images with higher homogeneity for separate classification and then overlapping (simply classify first then overlap) can effectively reduce the negative influences of phenology and sensor differences on later classification. The final images used for the target year 2000 were from July 2000 and July 2001; for 2010, they were from August 2010 and June 2008; and for 2015, they were from August 2015 and July 2016. A topographic illumination correction method proposed by Poortinga, et al. [
24] was applied to all the images to reduce visual interpretation errors in the following LULC classification.
2.4. LULC Classification
We encountered two challenges in the LULC classification of the study area. One was that the LULC was fragmented, and the other was the influence of cloud, fog, and haze. After many attempts to classify the satellite images using various common classification methods (e.g., supervised classification, unsupervised classification, and segmentation classification), decision tree classification was employed to minimize the errors and inconsistencies in images from different years and classify the images. The indices and calculation methods are shown in
Table 1. Considering convenience during later modeling and the reality of LULC in Chongqing capital, the satellite images were classified into six LULC categories: urban, cropland, forest, shrubland, grass, and water, which refer to the categories defined by the Chinese Academy of Sciences and were slightly adjusted [
25].
Figure 3 illustrates the classification process, taking August 2010 as an example. The digital numbers of each band were linearly stretched from 0–255 and then used to calculate the indices. The thresholds of the indices of all six images are presented in
Table 2.
The normalized difference vegetation index (NDVI) was first used for the separation of vegetation and non-vegetation. The modified normalized difference water index (MNDWI) is very effective in extracting water bodies [
26]. It was thus used to distinguish between urban landscapes and water bodies in the non-vegetation areas. Originally, the normalized difference building index (NDBI) was used to extract buildings e.g., [
27,
28], as it can also be used in conjunction with the NDVI to distinguish vegetation, built-up areas, and bare soil. Through visual interpretation of the satellite images and comparison with high-resolution historical images provided by Google Earth, we found that a small amount of grassland with low vegetation cover and soil as a background existed in the study area. The NDBI values of these grasslands were significantly different to those in areas with high vegetation cover, like forests. Urban areas were extracted using the NDVI in the previous step. Therefore, the NDBI was used to distinguish grasslands from areas with higher vegetation cover. The subdivision of the areas with high vegetation cover was challenging because their spectral characteristics are similar. As overly steep land is generally unsuitable for agricultural planting, slope degree was often used to distinguish woodland from cropland, e.g., [
29,
30]. As mountains and hills dominate Chongqing and several projects on the conversion and consolidation of sloping land have been underway for many years, e.g., the “Grain for Green” program [
31], we were careful to use 30° as a threshold to initially extract some forest and shrubland areas. The relatively high threshold of 30° was used to minimize the risk of classifying cropland as woodland. After analyzing the band values, we found that the band values from infrared to red differed between forest, shrubland, and cropland. Therefore, the simple calculation of the band value “near infrared (NIR) + shortwave infrared (SWIR) − 2 × RED” was used to separate forest, shrubland, and cropland. Forests usually had the lowest “NIR + SWIR − 2 × RED” values, while shrublands had the highest, and the cropland values usually fell in between. Forests and shrublands with slopes less than 30° were distinguished using a fixed value (100), which was determined by trial and error. For areas with a slope greater than 30°, two thresholds were used to divide them into forests, croplands, and shrublands in turn (
Figure 3,
Table 2). In this study, the GREEN, RED, NIR, and SWIR bands correspond to b2, b3, b4, and b5 in Landsat 5 images, and b3, b4, b5, and b6 (SWIR 1) in Landsat 8 images, respectively.
The thresholds for the decision tree were based on the statistical analysis of the regions of interest (ROIs) in each category. ROIs were randomly selected from the segmented images by referring to high-resolution historical images on Google Earth. Six ROIs, corresponding to the six classification categories, were carefully selected to ensure balance in the quantity of each category and uniformity of distribution. The means and standard deviations (SD) of the indices used in the decision tree for each ROI were calculated. The distance between the threshold and the mean value of the target extraction category should preferably be at least twice the SD and should not be less than one SD in special cases, to ensure that the different categories can be easily distinguished. The meticulous extraction of the ROIs was used to determine the threshold values and for accuracy validation.
Finally, the classified image was generalized to remove isolated pixels using a 3 × 3 kernel-mode filter then validated using the ROIs as references. The overall accuracy for the 2000, 2010, and 2015 images was 80.2%, 86.3%, and 81.1%, respectively, while the corresponding kappa index of agreement was 0.74, 0.81, and 0.76, respectively.
Figure 4 shows the classified LULC images for 2000, 2010, and 2015.
2.5. Land Change Modeling
The most important steps in land change modeling using LCM in TerrSet are the identification of major LULC transitions and the screening of corresponding explanatory variables. The rest settings such as new development plan by the government, infrastructure changes, road growth, and adjustment of change difficulty are selected depending on the actual needs. Our modeling process is comprised of the following major parts.
2.5.1. Identification of Major LULC Transitions
It is difficult and often unnecessary to analyze and simulate all LULC transitions. The LCM allows users to focus only on major transitions or those which may have a significant influence. Therefore, change analysis was first conducted to identify the major LULC transitions. We used land-use accounts and a change matrix to identify the main change trends and major transitions. Land-use accounts are an effective tool to explain the flow of each land category by explaining their gain, loss, persistence, net change, and turnover [
32]. In the land-use accounts table, “gain” means a new formation from other categories, “loss” means consumption by other categories and persistence means no change, “net change” is the gain minus the loss, and “turnover” is the sum of gain and loss, reflecting all areas that underwent changes. Thus, a land-use accounts table was first used to show the general change trends of each LULC category (
Table 3). The specific changes between LULC categories were displayed using a change matrix (
Table 4).
The selection of transitions for modeling ultimately depends on the research purpose. Usually, a threshold transition area is set to distinguish major transitions. In our study, after careful analysis of the land-use accounts and change matrix, we decided to use 20 km
2 as the threshold value, and transitions with an area less than this were ignored during change modeling. Seven major transitions were eventually included during change modeling (
Table 5). A further advantage of the LCM is that it provides another tool, the multi-layer perceptron neural network tool, which can model several or even all changes at once that have similar driving forces (with the exception of logistic regression) [
33]. Therefore, in this study, we divided the seven major transitions into three sub-models: urbanization, reclamation, and conservation, to represent three change trends. However, after many trials, the accuracy of the conservation model was always lower than 0.60. Thus, the transitions from cropland to forest and from cropland to shrubland were ultimately not merged into a sub-model but were calculated using separate models named “Conservation 1” and “Conservation 2,” respectively (
Table 5).
2.5.2. Determination of Explanatory Variables
Transition potential modeling is used to find suitable explanatory variables and generate transition potential maps based on the calculation of these variables [
34]. The preliminarily selected variables are shown in
Table 6, including evidence likelihood–normalized past changes, topographic factors (elevation and slope degree), the distance to each category, the distance from roads, and the distribution density of each category. The elevation and slope degree were calculated using the NASA Shuttle Radar Topography Mission Digital Elevation (30 m) dataset [
35], which was downloaded through the GEE. The road maps were from the National Earth System Science Data Center, the National Science & Technology Infrastructure of China (
http://www.geodata.cn, accessed on 7 July 2021). There were five categories of roads included on the maps: railways, expressways, national roads, provincial roads, and county roads. For the LCM, the roads needed to be divided into three major levels; thus, railways and expressways were grouped into primary roads, national and provincial roads were grouped into secondary roads, and county roads were classified as tertiary roads. The distance parameters were Euclidean distance and the map density values were created using a 7 × 7 kernel on TerrSet.
Not all explanatory variables must be employed in the final modeling. The selection of suitable explanatory variables can be divided into two steps. First, LCM provides a quick tool, Cramer’s
V, to measure the explanatory power of variables. Cramer’s
V ranges from 0 to 1; the larger the value, the stronger the potential association. In our study, variables with a Cramer’s
V less than 0.1 were eliminated (
Table 6). Second, the multi-layer perceptron neural network analysis provides a more comprehensive report, including the detailed explanatory power of each variable. According to their influence order and a stepwise backward-elimination analysis, variables with less influence were dropped to achieve a more parsimonious model (
Table 7). It should be noted here that “accuracy” refers to the model simulation accuracy when the corresponding variable is forced as a constant; thus, in general, the lower the accuracy, the higher the influence ranking of the variable. Once the accuracy of the selected explanatory variables was acceptable (larger than 70% in our study), the transition potential map could be created automatically by the modeler.
2.5.3. Incorporation of Government-Led Land Development
The Chongqing Municipal Government announced a land development plan for the capital urban area from 2007–2020 [
36], which provided a map showing the planned construction zones in the Chongqing capital. Therefore, we digitized the construction zones and converted them into an incentive layer for urban expansion (
Figure 5). Because the deadline for this development plan was 2020 and the target projection year was 2050, we added buffer zones of 2 km around the 2020 construction plans, indicating long-term development areas, and assigned them a value of 2 to indicate the higher urbanization probability within these areas. It should be noted that, before model accuracy assessment, the development plan layer was also imported into the model for land change simulation, but we did not incorporate the 2-km buffer zones around the construction areas because the planning period of the development plan, i.e., “2007–2020,” covered the later year “2010” and the verification year “2015”.
2.5.4. Validation of the Model
The projected 2015 LULC image was generated using the above settings. Meanwhile, a soft-prediction image was also created for validation. Soft prediction yielded a map indicating the change potential for LULC transitions. Thus, it can be used to quantify the predictive power of the model by comparing it to the image showing real change. The area under the curve (AUC, ranging from 0–1) was calculated by comparing the soft prediction image and the Boolean image showing the difference between the projected and real LULC images for 2015 [
37]. Usually, transition potential modeling requires many attempts (selecting different explanatory variables, etc.). Considering the objective of our study, when the AUC exceeded 0.7, the model was used for future projection.
2.5.5. Road Growth Settings
Over relatively short periods, road growth simulation is optional, for example, in the 2015 validation projection. However, in this study, the target year 2050 is 40 years after the later year of the model, 2010. The proximity to roads may be a strong factor in LULC change, especially considering that China is still a developing country and has a considerable demand for roads. The LCM provides a tool for dynamic road development which can help determine how road networks may grow. The logic behind the road development tool is that primary roads can be extended and grow secondary roads, secondary roads can also be extended and grow tertiary roads, and tertiary roads can just be extended. The primary roads are generated automatically according to the internal calculations of the tool. Users need to decide the “road spacing” and “road length” of the secondary and tertiary roads, reflecting the frequency with which the lower-grade roads grow along the higher-grade roads and their degree of extension. In our study, the “road spacing” and “road length” values for secondary roads were 12 and 5 km, respectively, while those for tertiary roads were 12 and 3 km.
2.5.6. LULC Projection
The LCM produces an LULC change map through a multi-objective land allocation procedure which determines the change potential of the involved transitions by maximizing the suitability of the land for all the objectives [
38]. The LULC transition quantities are calculated by an internal Markov model, then these transitions are allocated by a multi-objective land allocation algorithm. Each scenario is converted to a set of suitability maps for the decision process, and these can be loaded into the LCM to project the corresponding LULC map. The settings and conversions of the scenarios are described below in detail.
2.6. Scenario Settings and Conversions
The primary objective of this study was to explore the impact of the Eco-redline policy and assess various policy implications for the further improvement of policy implementation. In our analysis, two factors were integrated with the settings and analysis of future scenarios: 1) the spatial extent of the ERAs and 2) their management intensity. The “Normal Eco-redline” scenario was first set to represent a scenario in which the policy is implemented according to the current situation (
Table 8). The “Normal Eco-redline” scenario served as a baseline scenario. The “No Eco-redline” scenario was set to explore land-use consequences in the absence of the Eco-redline policy. The “Less ERAs” and “More ERAs” scenarios were derived by adjusting the “Normal Eco-redline” scenario; “Less ERAs” incorporated smaller ERAs than those in the “Normal Eco-redline” scenario, while “More ERAs” incorporated larger ERAs. The “Loose management” and “Strict management” scenarios refer to the adjustment of management intensity, i.e., strengthening and loosening management, respectively.
The above-mentioned assumptions for each scenario were translated into the model. We achieved the manipulation of the spatial extent of the ERAs by employing a 500-m buffer zone, based on the current ERAs (
Figure 6). The “Less ERAs” scenario subtracted a negative buffer zone, while the “More ERAs” scenario added a positive buffer zone. The LCM can introduce constraints or incentive coefficients for each LULC transition to simulate the difficulty of change. A value of 1 represents no impact while a value of 0 represents an absolute constraint. Values between 0 to 1 were treated as constraints and values larger than 1 acted as incentives. Thus, with regard to management intensity, the major LULC transitions were assigned different constraint or incentive coefficients to represent corresponding management intensities (
Table 9). For example, within the ERAs, cropland to urban transition is absolutely prohibited at normal management intensity, so it was assigned a 0 value, but it is slightly allowed under loose management, so it was assigned a value of 0.2. In the Strict management scenario, the transition of cropland to forest is highly encouraged in ERAs, so it was assigned a value of 2.
2.7. Future Projection for 2050
Before projection, recalculation stages must be specified. The period from 2010–2050 is 40 years, four times the length of the model period (2000–2010), so four recalculation stages were set in this study. The dynamic roads settings and recalculation stages were the same for all scenarios. The LULC images for all the scenarios were then generated, according to their different ERAs and management intensity coefficients.
2.8. Landscape Index
The landscape index is one of the most popular indicators to address the spatial and temporal characteristics, change trends, and driving factors of LULC. According to the explanation of the Eco-redline policy by the Chongqing Municipal People’s Government [
22], the primary protection targets of Chongqing’s Eco-redlines are forests, wetlands, and grassland ecosystems. Our land change modeling indicated that the water and grass areas in Chongqing capital are relatively small and stable, and the changes related to them were not included in the major transitions (
Table 5). Therefore, the forest (including both the forest and shrubland categories in this study) landscape indices are the best indicators to identify the effectiveness of the Eco-redline policy in Chongqing capital.
Many landscape indices have similar meanings and can be strongly correlated. To reduce this potential correlation, three aspects of forest patches representing distinctly different landscape features were considered in the current study. To be more precise, shape complexity, contrast with other patch categories, and aggregation of the same LULC class (i.e., forest) were calculated to compare the effectiveness of ERAs under the different scenarios. Shape complexity, which reflects the shape and size of a patch and their potential interaction, is highly related to ecological processes both within the patch and along its edges, e.g., [
39]. The degree of contrast reflects the quality of the microhabitats and climates at the edge of the patch, which is critical for the survival and migration of some species [
40]. In particular, an increase of urban patches around a forest, caused by urban expansion, may seriously affect the overall functioning of the forest. Therefore, it is often suggested that sufficient buffer space be left around protected areas to ensure their conservation value, e.g., [
41]. As a comprehensive indicator of patch distribution, forest aggregation is highly correlated with connectivity; thus, it can be used as an indicator of the quality of ecological corridors [
42]. Although there is still debate as to whether to maintain a few large forest patches (i.e., land sparing) or a large number of small forest patches (i.e., land sharing), the aggregation degree of an excessively fragmented landscape is generally smaller and is considered detrimental to ecological functions [
43]. In practice, it is usually difficult to quantitatively determine the optimal values or intervals of landscape indices. Considering that the topography of Chongqing capital is fragmented, land-use change occurs frequently and human–land conflict is severe. For forest patches, which dominate most ecological processes, relatively low shape complexity, low contrast, and high aggregation are thought to be preferable, and these served as the criteria when comparing the effectiveness of the various ERA and management-intensity scenarios.
From FRAGSTATS version 4.2 [
44], three class-level forest indices, the perimeter-area fractal dimension (PAFRAC), total edge contrast index (TECI), and aggregation index (AI), were selected to represent shape complexity, contrast, and aggregation, respectively (
Table 10). As the proportion of shrubland to forest was small (see
Table 3 and
Table 4), the shrubland was integrated into the forest for the calculation of the landscape indices to simplify the calculations and comparisons. In addition, a table indicating the magnitude of the contrast is needed in the calculation of the contrast index, with values from 0–1 reflecting increasing contrast (
Table 11). An exhaustive sampling strategy was selected, using 10 × 10-km squares to divide the whole area of Chongqing capital into 69 tiles. The landscape indices of the forest patches in these 69 tiles were calculated to represent the forest landscape characteristics of the whole of Chongqing capital.
2.9. Statistical Analysis
Owing to the large differences in LULC in the 69 tiles, a normality test was first conducted on the landscape indices. According to the results of the Shapiro–Wilk test, none of the indices were normally distributed (p < 0.05). The nonparametric Friedman rank-sum test was therefore used to examine the significance of the variance in all pairwise combinations of the six groups of landscapes indices. All analyses were carried out on SPSS v22.0 (IBM, USA). Box plots were used to display the median, first quartile, third quartile, minimum, and maximum values of the landscape indices for each scenario. The letters on the boxes indicate the significance of statistical analysis, while indices sharing the same letter indicate that the differences between them were not statistically significant (p < 0.05).