1. Introduction
Nature reserves are intended to protect natural resources, such as rare wildlife, and to produce value in the form of ecological services, scientific research, and cultural reserves. Therefore, they are effective for protecting ecological resources and biodiversity [
1]. Nature reserves are typically divided into functional zones in China, namely core, buffer, and experimental zones [
2]. The control requirements of the three zones are similar; thus, the core and buffer zones are typically merged into core protected areas. The experimental zone is retained as a general control area [
3].
Designing nature reserves is complex; the goal of these designs is to scientifically divide functional zones and to improve conservation efficiency [
4,
5,
6,
7]. Ecological and operational algorithms can be combined to further research this area [
8]. Most studies first predict the geographical distribution of species and then use planning algorithms to select grids because predicting species’ distribution is key for humans to make rational use of natural and economic resources in constructing protected areas [
9,
10]. Species distribution modeling methods use presence–absence data, climate data, and data on biotic interactions [
4,
11,
12,
13,
14,
15]. The Maxent model is the most widely used model for a highly accurate prediction using a small amount of data [
7].
Network selection algorithms are applied for nature reserve planning; for example, the prediction results for species distribution and other multisource remote sensing data can be used to achieve a multiarea division. Some studies have used intelligent heuristic algorithms, such as genetic algorithms, simulated annealing algorithms, and the tabu search algorithm, to determine and divide the scope of protected areas [
16,
17,
18]. However, a heuristic algorithm cannot ensure that the obtained solution is optimal [
19]. By contrast, optimization algorithms can achieve an optimal objective function solution under various constraints; integer programming (IP) algorithms are commonly used. An IP model can transform the conservation objectives into objective functions and conditions; for example, nonoverlapping regions can be transformed into modeling constraints. Therefore, some studies have used IP in network selection to express planning nature reserves as an optimization problem [
20].
In early operations research, researchers proposed using the set covering problem (SCP) model and maximal covering problem (MCP) model [
21,
22] with species numbers as input. Moreover, studies have used the 0–1 programming model and graph theory as the core ideas associated with spatial information, such as spatial continuity, site minimum distance, compactness, or both continuity and compactness, to develop linear IP nature reserve selection models [
23,
24,
25,
26,
27]. These models represent the relationships between grids by including linear inequalities to achieve spatial continuity and to increase the constraints on the total perimeter and grid number in the objective function to achieve a more compact distribution of the selected grids. Other researchers have considered the effects of environmental changes on the benefits of protection and have improved the transshipment model and other dynamic IP models [
28,
29,
30,
31]. They indicate that dynamically updating the cost function and terminal set or changing the selection objects in each period according to changes in the social economy and environment can produce more realistic results in the planning process than a general IP can. In addition, a previous study proposed a programming algorithm for certain ecosystem attributes, such as those of forest ecosystems, that combines biodiversity and spatial characteristics. Lin et al. [
32] associated ecological information with spatial properties to generate a space–ecology SCP (SeSCP) for designing a reserve network in Daiyun Mountain, China. Their SeSCP model outperformed three popular site selection models: the species SCP (SSCP) problem, tail length problem (TLP), and network flow problem (NFP). Moreover, they further proposed a dual-flow mechanism to overcome the problem of sparseness in site selection and to outperform three popular site selection models: SSCP, NFP, and simulated annealing (SA) [
33]. These methods can all be used to systematically and successfully design nature reserves, but most algorithms attempt to achieve spatial continuity and thus choose numerous grids with little protection benefit, resulting in inefficient land use.
Moreover, the few studies that have employed a multiperiod reserve planning method have focused on the characteristics of specific ecosystems and on dynamic changes in the biological distribution. Jafari et al. [
30] first investigated the reserve network design problem in the context of a multiperiod decision problem. Their proposed 5-year planning scheme attempts to achieve a certain total value, and the value for each year could differ. Finally, selections are compared in these different periods to choose a final optimal strategy. The central concept of dynamic programming is obtaining the best decision scheme by considering all possible changes in a certain period. However, applying dynamic programming to the multiperiod zoning of protected areas has three main requirements that are difficult to satisfy: (1) appropriately selecting grids in each period that consider their current status and potential value to determine the optimal global solution, (2) maintaining the continuity of the protected areas in different periods—that is, the selected girds in each period should remain spatially connected with those in the previous period—and (3) determining the zoning of the protected area with dynamic multiphase programming. To the best of our knowledge, no relevant model has achieved all three of these requirements. Thus, a discussion of dynamic programming of the multiperiod nature reserve by fully considering the dynamic changes in ecological suitability value and the overall benefits of the nature reserve is necessary.
In this paper, a multiperiod dynamic programming algorithm (MDP) for nature reserves is proposed. The MDP algorithm includes a comprehensive objective function with three aspects: the benefits of the nature reserve, the perimeter of the nature reserve, and the number of sites. The MDP algorithm achieves spatial continuity, compaction, and zoning using the primary optimization constraints and multiperiod constraints based on virtual points; it also determines an optimization scheme for the nature reserve function area. Moreover, the MDP algorithm still includes three-zone planning as the zoning basis to achieve more accurate and adequate protection. The research results significantly improve reserve planning algorithms and scientific methods for optimizing nature reserves.
2. Study Area and Data
2.1. Study Area
The study site, Quanzhou Bay Estuarine Wetland Nature Reserve (QBEWNR), is located in Fujian Province at the southeast coast of China (latitude:
–
, longitude:
–
E); the total land area of the study site is 7065.22
. The area is typical of subtropical estuarine wetlands in China and is rich in aquatic resources. Maps and images of the reserve are presented in
Figure 1. A total of 193 species, 41 families, and 12 orders of fish have been identified in the reserve; most species are perch flatfish. The ecological types include migratory, estuarine, and bottom-dwelling fishes. The Chinese sturgeon is listed as a national class I wildlife protected animal and is listed as vulnerable in the
China Red Book of Animals: Fish. Moreover, finless porpoises and branchiostoma belcheri are listed as national class II wildlife protected animals. The park also has 10 types of marine mammals; among these, the Chinese white dolphin, finless porpoise, grey dolphin, and sperm whale are listed as national wildlife protected animals.
Aquatic animals leave the protected areas due to their natural migration behavior and human factors, such as the presence of cross-sea bridges and conduct of marine economic operations. Therefore, protecting the rare aquatic animals in existing protected areas is critical. In this study, the ecological function area of Puganqiangcheng estuary, in which the Jinjiang River and Luoyang River converge into the sea, as the primary research area to develop a better protection plan and to improve the protective effect for rare aquatic animals.
The government of Quanzhou city has increased the ecological restoration of the QBEWPNR, but the survival of rare marine animals has still been affected by human activities. Thus, the original functional area lacks sufficient protection [
34].
Figure 2 presents the Puguncheng estuary wetland ecological function area, which is an essential functional zone of the QBEWNR. The figure reveals that only a few rare aquatic animals appeared in the core or buffer areas; most of their activities are outside the functional zone. The results indicated that the existing ecological function zone of Pugangcheng does not include the primary living habitats of rare aquatic animals and does not sufficiently protect these rare aquatic animals.
2.2. Data
We took satellite images with a resolution of 8 m × 8 m captured by the Ziyuan-3 (ZY-3) satellite in 2016 as input and performed preprocessing, including geometric registration, atmospheric correction, and crop. Moreover, we selected the factors of the natural environment and human activity affecting the distribution of rare aquatic animals, including elevation, land use type, and nearest distance to the reserve. Land use types were divided into eight types: forest, mangrove, dry field, paddy field, building, bare land, road, and water using supervised classification and visual interpretation. The nearest distance to infrastructure that significantly affects the activities of rare aquatic organisms, namely paddy fields, roads, and buildings, was determined using the multiring buffer analysis function of ArcGIS 10.2. Grid sections were used as the base unit, and Quanzhou Bay was divided into 200 m × 200 m square grids using ArcGIS 10.2; a total of 5969 grid sections were produced.
The data of rare aquatic animals were derived from a field scientific investigations performed in 2016–2017 (
Figure 3). These data were collected by irregular fishing and observation at historical animal occurrence points provided by the Quanzhou Marine Bureau, the nature reserve, and local fishers. The distribution data described the occurrence of rare aquatic animals, such as the time and place of the observation, the species, and the number of animals observed. After data screening and elimination, we obtained 17 rare aquatic animal distribution sites. We assumed that if rare aquatic animals were successfully protected, other aquatic animals would also be protected because rare aquatic animals have strict habitat requirements. Therefore, we only analyzed rare aquatic animals in this study.
We randomly selected 75% of the distribution sites as the training set and 25% as the testing set based on the distribution data of the rare aquatic animals. Moreover, we analyzed the relationship between rare aquatic animals and environmental factors using the Maxent model and then predict the potential distribution probability of the rare aquatic animals; this number was used as the ecological suitability value (
Figure 4). The included environmental factors were land use type, a digital elevation model, and the nearest distance. The test area under the receiver operating characteristic curve (AUC) value generated from the Maxent model was 0.856, indicating that the model has high prediction accuracy. The potential distribution probability, as shown in
Figure 4, was used as the ecological suitability value in the first period.
3. Methods
An MDP algorithm for nature reserves based on cellular automata (CA), network flow, and other fundamental theories is proposed and developed to address the weaknesses of existing algorithms.
3.1. CA
In protected areas, the ecological benefits of a grid are affected by the surrounding grids. That is, the ecology of the central grid is also optimized if the surrounding grids have good ecology. Moreover, the execution of this process is similar to the theory of the CA model in which each cell is modified based on the state of neighboring cells according to certain system-level rules in the cellular space with discrete units and finite states. Therefore, if the ecological suitability values of adjacent grids are all higher than that of the central grid, the ecological suitability values of the plot is assumed to correspondingly increase. This dynamic updating system for the CA is designed based on this assumption and is expressed as follows:
Here,
S is the ecological suitability value of the grids, as presented in
Figure 4,
is the ecological suitability value of grid
i in the period
t (such that
is the ecological suitability value of grid
i in period 1), and
is the ecological suitability value between grid
i and its adjacent grid
j. If the ecological suitability value grid
i is less than the adjacent grid
j,
; otherwise,
.
k is a random number that increases the ecological suitability value. The ecological suitability values in the rest periods are produced using the CA dynamic update mechanism in Equation (
2). Specifically, the ecological suitability value
is the potential distribution probability assigned to each grid, as shown in
Figure 4. We then analyze the relationship between grid
i and its neighborhood using Equation (
2) to adjust the ecological suitability value of grid
i and to generate the updated ecological suitability value,
, in period 2. These steps are repeated until the final period
t to obtain
.
3.2. Optimized Objective Function
The objective function of the MPD algorithm has three goals: (1) the maximization of the protection benefits, (2) the minimization of the total perimeter, and (3) the minimization of the number of buffer zones. Finally, the overall protection level was maximized using the multiperiod procedure. The relevant equations are as follows:
We consider the maximization of the ecological value as the objective function of protection benefits in the designed multiperiod programming; this maximization is expressed as follows:
Here, E is the total ecological suitability value of the nature reserve. indicates whether grid i is selected as a core area in period t and is 1 if selected and 0 otherwise. T is the total number of periods. I is the total number of grids.
Moreover, the minimization of the perimeter of the selected grids of core and buffer areas is an additional goal of the objective function. The objective function for minimizing the perimeter of the selected core grids is as follows:
Here, M is the total perimeter of the nature reserve, and is a binary variable indicating whether grid i is selected as part of the core area; it is 1 if selected and 0 otherwise. indicates whether grid i and its adjacent grid j are selected as core areas and form a connection; if so, ; otherwise, .
Furthermore, we minimize the number of grids selected as buffer areas to limit the total number of selected grids; the minimization function is as follows:
Here, N is the total number of buffer grids, and is a binary variable. If grid i is selected as a buffer area, ; otherwise, . We only consider the minimization of the perimeter of the selected core and buffer areas (grids) in the final period because we must first ensure the overall spatial compactness of the core areas, and the buffer areas must surround the core areas.
Finally, we integrate these subobjectives into an objective maximization function formulated as follows:
Here,
is an empirically determined parameter. Equation (
6) is a modification of Equations (
4) and (
5) from minimization to maximization functions with a multiplication by
.
3.3. Basic Constraints
Two basic constraints were included: cost and type constraints.
3.3.1. Cost Constraints
Grid selection: Each grid can only be selected as a core area once per period as follows:
Grid Aggregation: The construction cost of nature reserves was used to limit the number of grids selected in each period; this limitation is expressed as follows:
Here, is the total number of grids selected in each period.
Minimum Protection of Core Area: The total ecological suitability value of grids selected as core areas in each period must be higher than the minimum protection requirement. The minimum protection requirements are expressed as follows:
Here, is the conservation proportion of the target species in the core area.
3.3.2. Type Constraints
Statistics of Core Area: The number of times a grid
i was selected as a core area over all periods dynamically can be formulated as follows:
Here,
is a binary variable of 0 or 1, as defined in Equation (
4).
Continuity of the Core Area: Core areas should be selected to achieve both spatial continuity and compactness to provide stable breeding sites for species. If grid
i and its adjacent grid
j are both selected as core areas, they must have a spatial connection. However, if one of these grids is not a core area, this connection may not exist. Core area continuity is expressed as follows:
Here,
C and
L are the binary variables defined in Equation (
4).
Partition Attribute: The grid
i cannot be selected for both the core and the buffer areas simultaneously.
Here, C and B are binary variables indicating inclusion in the core and buffer areas, respectively.
Partition Connection: A key requirement in nature reserve programming is ensuring that the buffer area surrounds the core area. If grid
i is selected as a core area, its adjacent grid
j must either be a core area or buffer area. This condition is formulated as follows:
Here, C and B are binary variables indicating inclusion in the core and buffer areas, respectively.
3.4. Multiperiod Constraints Based on Virtual Point
The selected grids must be continuous in all periods to achieve overall spatial continuity of the nature reserve MDP.
Figure 5 presents an example of grid selection over three periods indicated by red, green, and blue. Each grid square in
Figure 5 is labeled at the upper left corner, and its ecological suitability value is centered. In this example, three grids in each period are selected: blue in period 1, green in period 2, and red in period 3. Although the selected grids in each period are continuous both within and between periods (i.e., local and global continuity, respectively), grid 8 (light grey, ecological suitability value of 0.90) is not included to ensure continuity. This example reveals that the balance between spatial continuity and ecological suitability is critical when selecting grids in reserve programming.
Therefore, we introduce the concept of the virtual point in the MDP to balance the conflict between spatial continuity and ecological suitability value when selecting grids. A virtual point is a virtual grid independent of the grids but connected to each grid already selected in the nature reserve, as indicated by the orange square in
Figure 6. Grids in each period are selected beginning at the virtual point; this enables selecting grids that are locally discontinuous but globally continuous. For example, the grids selected in periods 2 (green) and 3 (red) are not locally continuous, but the final selection is globally continuous, and the total ecological suitability value is 7—0.20 greater than the result of
Figure 5. Thus, the protection benefits can be maximized by setting virtual points to ensure global continuity in each period. We set the relevant constraints of the virtual points as follows:
Number of Starting Points: The Netflow algorithm specifies the number of starting points to limit the final selection of set regions; the constraint is formulated as follows:
Here, is the flow from the virtual point to grid i in period t. indicates that only one grid is connected to the virtual point in period one to limit the number of reserve areas to 1.
Limitation of Flow Traffic: If grid
i is not selected in period
t, the inflow number of grid
i is 0; otherwise, the maximum number of the inflow from the adjacent grid
j to grid
i is (
M − 1). The constraint formula is expressed as follows:
Here, is the total number of flows from core grid j to its adjacent core grid i in period t. is a binary variable; indicates that grids i and j are connected; otherwise, . is the set of grid numbers adjacent to grid i. M is a constant and can be set to an arbitrary value but should be larger than the total number of grids.
Quantity Relationship between Inflow and Outflow: The difference between the outflow and inflow of grid
i in period
t must be greater than or equal to
. Its constraint is formulated as follows:
Here, is the total number of flows from core grid i to an adjacent core grid j in period t. If grid i in period t is not connected to the virtual point, its outflow is greater than or equal to the inflow + 1.
Continuity constraint between the virtual point and the core grid: If grid
i in period
t is connected to the virtual point, it must be selected as a core grid. The formulation is expressed as follows:
Multiperiod grid continuity: If grid
i connects to the virtual point, it must connect to grid
j, selected as a core grid in the previous period, to ensure the grid’s continuity between each period. We express the constraint as the following formulation:
Here, indicates whether grid j adjacent to grid i is selected as a core area in period . If grid j is selected in period t − 1, ; otherwise, .
6. Conclusions
The distribution of rare aquatic animals simulated by Maxent was used to determine ecological suitability values in the QBEWNR. We used CA to simulate variations in ecological suitability values during different periods of partitioning the QBEWNR. Finally, we proposed a novel MDP algorithm based on network flow theory to design the nature reserve. The MDP algorithm satisfies the primary optimization constraints, such as cost, partitioning, and continuity.
Moreover, we introduced virtual points in the MDP algorithm to design multiperiod constraints, such as numerous starting points, networked traffic, and continuity between virtual points and core sites. The MDP algorithm maximized the ecological suitability value, minimized the perimeter, and minimized the number of sites for the reserve and thus rationally selected a protection scheme for the QBEWNR over an extended period.
The results reveal that the total ecological suitability value (212.92) of the reserve generated by the MDP algorithm was substantially higher than that of the original reserve (111.35). Moreover, the number of selected sites (372) selected by MDP was substantially lower than the original reserve (547), demonstrating the superior performance of the MDP algorithm.