Optimization of Land Use Based on the Source and Sink Landscape of Ecosystem Services: A Case Study of Fengdu County in the Three Gorges Reservoir Area, China

: Promoting the preservation and appreciation of ecosystem services is an important value guide for land use optimization. In this research, Fengdu County in the Three Gorges Reservoir Area was selected as the focus of a case study. From the perspective of the source and sink landscape of ecosystem services, a MOP model and FLUS model were used to optimize the areas of various land use types and the spatial conﬁgurations of those land use types in the study area in 2035 under a strict ecological constraint (SEC) scenario, a moderate ecological constraint (MEC) scenario, and a relaxed ecological constraint (REC) scenario. We also superimposed and adjusted the results of land use optimization under the three ecological constraint scenarios, and obtained land use regionalization results that integrated multiple scenarios. The results indicated that (1) there were large differences in the areas and spatial distributions of the source and sink landscapes under the three scenarios. Under the SEC scenario, the important source landscapes (ISLs), common source landscapes (CSLs), and sink landscapes (SLs) areas covered 1676.62 km 2 , 1190.43 km 2 , and 33.81 km 2 , respectively. A large area of the CSLs and a small area of the SLs were transformed into ISLs area, and the degree of fragmentation of the landscape was low. Under the MEC scenario, the ISLs, CSLs, and SLs areas covered 1609.22 km 2 , 1241.60 km 2 , and 49.74 km 2 , respectively. The development of the source landscapes and sink landscapes was similar, and the degree of fragmentation was moderate. Under the REC scenario, the ISLs, CSLs, and SLs areas covered 1603.96 km 2 , 1243.32 km 2 , and 53.58 km 2 , respectively. A large area of CSLs was transformed into SLs area, and the degree of fragmentation was high. (2) Fengdu County was divided into seven types of areas: ecological conservation area; agricultural production area; construction optimization area; construction-ecological area; ecological-agricultural area; agricultural-construction area; and integrated development area. The results of this study can provide references for the territorial spatial planning and management of ecological barrier zones.


Introduction
Ecosystem services (ESs) are the sum of life-sustaining products and services that humans acquire from ecosystems-these products and services are closely related to human well-being and sustainable development [1]. Ecosystems support and maintain balance in areas occupied by humans by regulating the air and water quality, and maintaining biodiversity [2]. Ecosystems also provide the food and raw materials needed to sustain life and Land 2021, 10, 1242 3 of 24 source pollution [31], and improving water resource utilization efficiency and economic income [32]. Land use optimization models include quantity optimization models and spatial configuration optimization models. At present, there are three main types of quantity optimization models. The first category is traditional mathematical optimization models, such as gray linear programming models [33] and multiobjective programming (MOP) models [25,34]. The second category is heuristic algorithms, such as genetic algorithms [35], simulated annealing algorithms [36], ant colony algorithms [37], and particle swarm algorithms [38]. The third category includes other methods, such as Markov chains [39], system dynamics models [40], and neural network models [41]. Among these models, the MOP model has the advantages of flexibility, practicability, and high credibility, and can address problems associated with land use quantity optimization in dynamic situations and for multiple objectives and plans. Spatial configuration optimization models include the cellular automata (CA) model [42], multiagent model [43], conversion of land use and its effects at small regional extents (CLUE-S) model [27], future land use simulation (FLUS) model [44], and so on. Moreover, heuristic algorithms, such as the simulated annealing algorithm [45], particle swarm algorithm [46], and ant colony algorithm [47], are commonly used. The FLUS model adopts an adaptive inertial competition mechanism based on roulette selection: this mechanism can effectively cope with the uncertainty and complexity of land use change under the interaction of the natural environment and human activities, and the simulation accuracy achieved with this model is high [44]. Generally, although the existing research methods have gradually matured, leading to more robust results, the following shortcomings still exist: (1) Scholars mainly focus on the application of source-sink landscape theory in the research of specific ecological processes, but no scholars have applied this theory to the research of land use optimization. (2) Existing research has mainly focused on technical methods and optimization goals, and insufficient attention has been given to the guiding theory of land use optimization. Moreover, the constraint conditions and optimization rules considered in previous land use optimization studies have been focused on mostly macro constraints, such as the numerical values associated with land use planning, whereas policy-oriented constraint conditions, such as those associated with ecological civilization construction and high-quality development, have not been considered frequently enough. (3) Existing research often regards improvements in the functions of ESs as the objective of land use quantity optimization, and there is a lack of research on how to promote the preservation and appreciation of ESs through the spatial optimization of land use.
In this research, we selected Fengdu County as a study area, and combined the sourcesink landscape theory with ES theory to construct a land use optimization framework. Moreover, we comprehensively considered the constraints imposed by the ecological environment and intensive development in the region. A MOP-FLUS coupling model, which combined the MOP model and FLUS model, was used to optimize the areas and configurations of land use types in 2035 under a strict ecological constraint (SEC) scenario, a moderate ecological constraint (MEC) scenario, and a relaxed ecological constraint (REC) scenario. We also superimposed and adjusted the land use optimization results for these three ecological constraint scenarios, and developed a land use regionalization plan that integrated multiple scenarios. Finally, we explored the role of the results in the functional optimization of ESs and space governance. The key scientific question of this research was how to scientifically and reasonably optimize the county-level land use pattern through the combination of ESs and source-sink landscape theory. The results of this research enrich the theory, perspective, and model of land use optimization research. At the same time, it also provides references for future territorial spatial planning in ecological barrier zones.

Study Area and Data Sources
The Three Gorges Reservoir Area (TGRA) is an important ecological barrier in the Yangtze River Basin (Figure 1). Its ecological environment is directly related not only to the long-term safe operation of the Three Gorges Dam and the provision of stable support for Land 2021, 10, 1242 4 of 24 millions of immigrants, but also to the ecological security and sustainable development of the whole Yangtze River Basin [48]. Fengdu County, ranging from 107 • 28 -108 • 12 E and 29 • 33 -30 • 16 N, is a typical mountainous county in the ecological barrier zone of the TGRA, Southwest China. The predominant landform in this region is mountains, followed by hills, and only a few flat mountain and river valleys-generally, the terrain is high in the southern parts of the county and low in the northern parts. The main types of land use are cultivated land and forest, with abundant forest resources and diverse ecosystems. The county covers an area of 2900.86 km 2 , with 2 subdistricts, 23 towns, and 5 townships. In 2018, the per capita GDP of Fengdu County was CNY 40,400, and the population urbanization rate was 46.48%. Since the start of the 21st century, the immigrant population in Fengdu County has become larger than the nonimmigrant population, and the total permanent population decreased from 671,100 in 2000 to 585,200 in 2018. The departing population mainly originates from rural areas. Between 2000 and 2018, the rural population decreased by 210,600. Fengdu County is faced with multiple problems, such as immigration relocation, degradation of the countryside, rapid urban expansion, and deterioration of the ecological environment. Therefore, there is an urgent need to allocate and optimize the land resources of Fengdu County according to scientific rationale.

Study Area and Data Sources
The Three Gorges Reservoir Area (TGRA) is an important ecological barrier in the Yangtze River Basin (Figure 1). Its ecological environment is directly related not only to the long-term safe operation of the Three Gorges Dam and the provision of stable support for millions of immigrants, but also to the ecological security and sustainable development of the whole Yangtze River Basin [48]. Fengdu County, ranging from 107°28′-108°12′ E and 29°33′-30°16′ N, is a typical mountainous county in the ecological barrier zone of the TGRA, Southwest China. The predominant landform in this region is mountains, followed by hills, and only a few flat mountain and river valleys-generally, the terrain is high in the southern parts of the county and low in the northern parts. The main types of land use are cultivated land and forest, with abundant forest resources and diverse ecosystems. The county covers an area of 2900.86 km 2 , with 2 subdistricts, 23 towns, and 5 townships. In 2018, the per capita GDP of Fengdu County was CNY 40,400, and the population urbanization rate was 46.48%. Since the start of the 21st century, the immigrant population in Fengdu County has become larger than the nonimmigrant population, and the total permanent population decreased from 671,100 in 2000 to 585,200 in 2018. The departing population mainly originates from rural areas. Between 2000 and 2018, the rural population decreased by 210,600. Fengdu County is faced with multiple problems, such as immigration relocation, degradation of the countryside, rapid urban expansion, and deterioration of the ecological environment. Therefore, there is an urgent need to allocate and optimize the land resources of Fengdu County according to scientific rationale. In this study, the data needed to optimize the land use in Fengdu County included spatial data and statistical data. The data sources are listed in Table 1.  In this study, the data needed to optimize the land use in Fengdu County included spatial data and statistical data. The data sources are listed in Table 1.

Types Datasets Sources Resolution
Statistical data Socio-economic indicators (gross domestic product, urban resident population) and agricultural production indicators (output of farming, forestry, animal husbandry, fishery, and agricultural service industries) "Chongqing Statistical Yearbook", and "Fengdu Yearbook" -

Methodology
A flow chart of the methods applied in this study is shown in Figure 2: the approach involved three main parts. After data collection and preprocessing, the MOP-FLUS coupling model was used to optimize the areas and configurations of land use types in 2035 under the three ecological constraint scenarios. Then, we used overlay analysis to develop a land use regionalization plan that integrated the results from multiple ecological constraint scenarios. Socio-economic indicators (gross domestic product, urban resident population) and agricultural production indicators (output of farming, forestry, animal husbandry, fishery, and agricultural service industries) "Chongqing Statistical Yearbook", and "Fengdu Yearbook" -

Methodology
A flow chart of the methods applied in this study is shown in Figure 2: the approach involved three main parts. After data collection and preprocessing, the MOP-FLUS coupling model was used to optimize the areas and configurations of land use types in 2035 under the three ecological constraint scenarios. Then, we used overlay analysis to develop a land use regionalization plan that integrated the results from multiple ecological constraint scenarios.

Land Use Optimization Framework Based on the Source-Sink Landscape of ESs
ESs are the sum of various material products and nonmaterial services provided by an ecosystem [4]. The processes by which ecosystems provide various services to human Land 2021, 10, 1242 6 of 24 society are continually undergoing material and energy changes, and these changes can be regarded as an ecological process that includes several subprocesses [1,4]. With respect to supporting ESs, some landscapes act as "source", while others act as "sink", and the attribution of landscapes to these categories can be performed according to their roles in the ES supply process [49,50]. ESV is an important indicator to measure ES supply capacity [51]. Referring to the research of Zhang et al. [52], we investigated land use type and the ESV per unit area as the bases to identify the source and sink landscape of ESs in the study area. The results show that paddy fields, dry land, forest, grassland, water, and unused land are source landscapes, whereas urban land and rural residential areas are sink landscapes (SLs). In addition, the contribution to ES supply differs by landscape type. Forests, grasslands, and water provide the main ecological products, and their ESVs per unit area are much higher than those of dry land and paddy fields, which mainly provide agricultural products. Therefore, forests, grasslands, and water are further classified as important source landscapes (ISLs), whereas dry land and paddy fields are classified as common source landscapes (CSLs). In addition, although the ESV coefficient of unused land is low, it is also classified as a CSL.
With the rapid development of industrialization and urbanization, competition for space between sources and sink landscapes continues, promoting the change of land use patterns [53]. However, source landscapes are often at a disadvantage, which has allowed sink landscapes, such as urban land and rural residential areas, to continue to overtake other landscapes, resulting in a continuous decrease in source landscapes [54]. Therefore, it is necessary to balance the relationship between source and sink landscapes, and optimize regional land use patterns.
Land use optimization is the process of optimizing the areas and configurations of different land use types [55,56]. The optimization of land use type area allocates land resources to land use types with high efficiency under some constraint conditions to improve the overall benefits of the land [25]. Spatial land use optimization involves rearranging land resources in space, and allocating those resources to spatial units with higher suitability. This process is performed at two scales: land use spatial configurations optimization on a micro scale, and land use zoning on a macro scale [57]. The relationship between these optimization processes indicates that the optimization of the areas of land use types is determined by the spatial constraint of the spatial land use optimization, whereas the spatial land use optimization is determined by the quantity constraint imposed by the optimization of land use type area [55] (Figure 3).

Land Use Optimization Framework Based on the Source-Sink Landscape of ESs
ESs are the sum of various material products and nonmaterial services provided by an ecosystem [4]. The processes by which ecosystems provide various services to human society are continually undergoing material and energy changes, and these changes can be regarded as an ecological process that includes several subprocesses [1,4]. With respect to supporting ESs, some landscapes act as "source", while others act as "sink", and the attribution of landscapes to these categories can be performed according to their roles in the ES supply process [49,50]. ESV is an important indicator to measure ES supply capacity [51]. Referring to the research of Zhang et al. [52], we investigated land use type and the ESV per unit area as the bases to identify the source and sink landscape of ESs in the study area. The results show that paddy fields, dry land, forest, grassland, water, and unused land are source landscapes, whereas urban land and rural residential areas are sink landscapes (SLs). In addition, the contribution to ES supply differs by landscape type. Forests, grasslands, and water provide the main ecological products, and their ESVs per unit area are much higher than those of dry land and paddy fields, which mainly provide agricultural products. Therefore, forests, grasslands, and water are further classified as important source landscapes (ISLs), whereas dry land and paddy fields are classified as common source landscapes (CSLs). In addition, although the ESV coefficient of unused land is low, it is also classified as a CSL.
With the rapid development of industrialization and urbanization, competition for space between sources and sink landscapes continues, promoting the change of land use patterns [53]. However, source landscapes are often at a disadvantage, which has allowed sink landscapes, such as urban land and rural residential areas, to continue to overtake other landscapes, resulting in a continuous decrease in source landscapes [54]. Therefore, it is necessary to balance the relationship between source and sink landscapes, and optimize regional land use patterns.
Land use optimization is the process of optimizing the areas and configurations of different land use types [55,56]. The optimization of land use type area allocates land resources to land use types with high efficiency under some constraint conditions to improve the overall benefits of the land [25]. Spatial land use optimization involves rearranging land resources in space, and allocating those resources to spatial units with higher suitability. This process is performed at two scales: land use spatial configurations optimization on a micro scale, and land use zoning on a macro scale [57]. The relationship between these optimization processes indicates that the optimization of the areas of land use types is determined by the spatial constraint of the spatial land use optimization, whereas the spatial land use optimization is determined by the quantity constraint imposed by the optimization of land use type area [55] (Figure 3).  The areal proportion of various source and sink landscapes directly determines the total regional ESV [52,56]. Therefore, to optimize land use for ESs, and according to the areas of various land use types, regional land resources need to be allocated to land use Land 2021, 10, 1242 7 of 24 types with higher ESVs per unit area under constraint conditions. Moreover, the spatial configuration of source and sink landscapes also determines the spatial characteristics of landscape connectivity and fragmentation, which have large impacts on regional ecological processes and ES functions [58,59]. For example, when a source landscape is invaded by a sink landscape, there is an increase in the resistance to the migration and flow of species and energy between landscape types-this resistance is not conducive to enhancing ecological processes, leading to a reduction in the ES supply capacity of surrounding areas. That is, the expansion of sink landscapes into source landscapes has a negative spatial effect on ESs [60][61][62]. This finding indicates that controlling the expansion and infiltration of sink landscapes into source landscapes, and maintaining the integrity and continuity of ecosystems are effective ways to improve regional ES functions. Therefore, land use types need to be allocated to more suitable spatial units through spatial land use optimization and land use zoning, thus, forming a spatial configuration more supportive of ESVs.

Optimization of the Areas of Land Use Types Based on the MOP Model
The MOP model is a scientific method for solving multiobjective optimization problems based on constraints and objectives [34]. The model includes three components: decision variables; objective functions; and constraint conditions. The formula of the objective function F(x) is as follows: where x j is the decision variable j; c j is the benefit coefficient of different land use types; s.t. represents the constraint conditions; a ij is the coefficient corresponding to variable j in the constraint i; b j is the constraint value; and n is the number of decision variables. The following land use types (eight categories) are used as the decision variables: x 1 , paddy field; x 2 , dry land; x 3 , forest; x 4 , grassland; x 5 , water; x 6 , urban land; x 7 , rural settlement; and x 8 , unused land.

Scenarios and Objective Function Setting
The implementation period of the new territorial spatial plan is 2021-2035, so we take 2035 as the optimization year. Three development scenarios are established according to the strictness of the ecological constraints: (1) The SEC scenario represents comprehensive implementation of the following development concept: "to step up conservation of the Yangtze River and stop its over development". The goal of this scenario is to maximize ecological benefits by managing land use. (2) The MEC scenario represents coordination and unity between development and protection. In this scenario, the aim of land use is to balance ecological benefits and economic benefits. (3) The REC scenario represents the maximization of economic benefits in the context of rapid urbanization.
Ecological benefits are measured by the ESV per unit area of different land types, and the specific coefficients are determined according to the findings of Zhang et al. [52] ( Table 2). Economic benefits are measured by the economic output per unit area of different land use types ( Table 2). In particular, the output value of secondary and tertiary industries is used to measure the economic output of urban land; the agricultural output value after deducting the output value of mulberry, tea, and fruit is used to measure the economic output of paddy fields and dry land; and the output values of the agricultural service industries are used to measure the economic outputs of rural residential land. To determine the economic outputs of forests, grasslands, and water, we consider the output of the forestry, animal husbandry, and fishery industries. The economic benefits per unit area for these land use types are calculated as follows: We first derive the economic output data from 2000 to 2018, and then, we predict their values in 2035 by using a GM (1, 1) model. Finally, the mean value of the three predicted outputs for each land use type is divided by the corresponding area, yielding the economic benefit per unit area. However, because unused land has not yet been exploited for economic use, we assume that it does not generate any economic value, and its coefficient is set at 0.01. Based on the ecological benefit coefficients and economic benefit coefficients of each land use type, the objective functions are established as follows: where Z represents the total benefit of a land use type; x i is the area of land use type i; K i and P i are the ecological benefit coefficient and economic benefit coefficient of land use type i; and a and b are the weights of the ecological benefit objective and the economic benefit objective, respectively. It was determined from related studies [63,64] that a is 0.8, 0.5, and 0.2, and b is 0.2, 0.5, and 0.8 in the SEC scenario, MEC scenario, and REC scenario, respectively.

Constraint Conditions Setting
(1) Total area constraint: The total area of each land use type is equal to the total area of the study area (Table 3). (2) Planning objective constraints: According to the study area's land use plan and related standards, the areas of cultivated land, forest, urban land, and unused land are limited in 2035. (3) Ecological environment constraints: To ensure ecological security and promote the continuous optimization of the ecosystem, constraint conditions are set for four factors (the total area of ecological land, forest area, water area, and habitat quality). Grasslands and water are important for maintaining the diversity of species in ecosystems. From 1990 to 2018, the grassland area in the study area decreased, whereas the water area increased, but the areas of these two types land cannot be reduced or expanded indefinitely. We limit the areas of grassland and water to between the value in 2018 and the predicted value in 2035. In addition, according to the method for calculating the habitat quality index E m described in the "Technical Criterion for Eco-environmental Status Evaluation", the E m in 2035 must be greater than the E m in 2018. (4) Intensive development constraints: To promote the intensive use of land and improve the sustainability of social and economic development, we set constraint conditions from two factors (urban land area and rural residential area). According to the "Code for Classification of Urban Land Use and Planning Standards of Development Land", Fengdu County belongs to architectural climatic zone III, and the planned upper limit of its urban land area per capita is between 85 m 2 and 105 m 2 [65]. Therefore, the upper limit of the urban land area per capita in 2035 is set at 105 m 2 . In addition, although a large number of people from rural areas have left the region for work, the rural residential area has increased rather than decreased. It is very important to promote the withdrawal of idle homesteads-however, this withdrawal is affected by policies, funds, and farmers' willingness, and is a relatively slow process. According to related research, the theoretical potential of rural residential land reclamation is approximately 40% of the actual area [66]. Therefore, we limit the area of rural settlements in 2035 to between 60% and 100% of the area of rural settlements in 2018. The predicted values of the above indicators in 2035 were calculated by the GM (1, 1) model based on data from 2000 to 2018.

Constraint Types. Constraint Variables Constraint Conditions
Total area constraint The total area of each land use type x 1 + x 2 +, · · · , +x 8 = The total area of the study area Planning objective constraints Cultivated land area x 1 + x 2 ≥ planning value in 2020 Forest area x 3 ≥ planning value in 2020 Urban land area x 6 ≥ planning value in 2020 Unused land area x 3 ≤ planning value in 2020 Ecological environment constraints Total area of ecological land

Optimization of the Land Use Spatial Configuration Based on the FLUS Model
The FLUS model is a scientific method suitable for the simulation of future land use changes [67]. This model is composed of two main modules: (1) an artificial neural network (ANN)-based probability of occurrence estimation module; and (2) a self-adaptive inertial competition mechanism CA module with roulette selection.

ANN-Based Probability of Occurrence Estimation Module
An ANN is a machine learning model based on biological neural network simulation. It is usually used to simulate and calculate nonlinear functions with many variables. It can continuously fit the complex relationship between input data and training targets through a large number of learning and recall iterations to ensure the generation of a higher suitability probability distribution and establish the relationship between the probability of each land use type and the factors that determine land use type [44,67]. The ANN can be expressed as follows: where SP(p, k, t) represents the probability of suitability for land use type k on grid cell p at time t; ω j,k and sigmoid are the weight and activation function between the hidden layer and the output layer, respectively; and net j (p, t) is the signal received from grid cell p in the hidden layer at time t.
Based on the results of relevant studies [68,69], 14 factors considered driving forces were selected to build the ANN to determine the probability of suitability for each landscape type in this study. The natural factors include elevation, slope, annual average temperature, annual average precipitation, soil texture, and soil organic matter content. The socioeconomic factors include population and GDP. The distance factors include the distances to rivers, railways, highways, roads, residential areas, and town areas. Because of the uneven distribution of land use types in the study area, random sampling was performed to obtain training samples. Due to the large number of pixels in the study area, the sampling ratio was set to 1% of the effective number of pixels in the study area. Since the unit of the sampling ratio was one thousandth, the sampling ratio parameter was set to 20. Moreover, the number of hidden layers in the ANN was set to 12. Finally, the ANN model was applied to generate accurate images of the probability of suitability for each land use type ( Figure 4). The deeper the red is, the more suitable the landscape is.
the uneven distribution of land use types in the study area, random sampling was performed to obtain training samples. Due to the large number of pixels in the study area, the sampling ratio was set to 1% of the effective number of pixels in the study area. Since the unit of the sampling ratio was one thousandth, the sampling ratio parameter was set to 20. Moreover, the number of hidden layers in the ANN was set to 12. Finally, the ANN model was applied to generate accurate images of the probability of suitability for each land use type (Figure 4). The deeper the red is, the more suitable the landscape is.

Self-Adaptive Inertial Competition Mechanism CA Module
The FLUS model introduces a self-adaptive inertial competition mechanism into the CA model, and combines the probability of occurrence, the influence of the neighborhood between cells, the conversion cost, and the inertia coefficient to calculate the total probability of the conversion of each land use type in each cell. Then, the roulette mechanism is used to reflect the competition among various land use types and determine the conversion and distribution of land use types in a cell to achieve more accurate land use change simulation [70,71]. The equation is expressed as follows: where , indicates the probability of grid cell being converted from the initial land use type to the target land use type at iteration time ; indicates the inertia coefficient of land use type at iteration time ; , is the neighborhood effect of land use type on grid cell at iteration time ; and → represents the conversion cost of the initial land use type to the target land use type .
Spatial constraints are areas where a land use type will not change, and three types of spatial constraints are included in this research: ecological conservation redline areas; basic farmland conservation areas; and ecological safety network areas ( Figure A1). The ecological safety network is constructed based on landscape security pattern theory [72]. First, the Getis-Ord Gi* statistics are calculated by ArcGIS software version 10.7 to describe the spatial locations of cold spots and hot spots of ESVs in Fengdu County, and the hot spots with a confidence level of more than 90% are considered ecological sources and extracted [73]. Second, a minimum cumulative resistance model is used to construct an ecological resistance surface of Fengdu County, and Linkage Mapper software is used to calculate the distance path with the lowest cumulative consumption and, then, identify ecological corridors [74]. Finally, we analyzed buffer zones around the corridors and found

Self-Adaptive Inertial Competition Mechanism CA Module
The FLUS model introduces a self-adaptive inertial competition mechanism into the CA model, and combines the probability of occurrence, the influence of the neighborhood between cells, the conversion cost, and the inertia coefficient to calculate the total probability of the conversion of each land use type in each cell. Then, the roulette mechanism is used to reflect the competition among various land use types and determine the conversion and distribution of land use types in a cell to achieve more accurate land use change simulation [70,71]. The equation is expressed as follows: where TProb t p,k indicates the probability of grid cell p being converted from the initial land use type to the target land use type k at iteration time t; Intertia t k indicates the inertia coefficient of land use type k at iteration time t; θ t p,k is the neighborhood effect of land use type k on grid cell p at iteration time t; and sc c→k represents the conversion cost of the initial land use type c to the target land use type k.
Spatial constraints are areas where a land use type will not change, and three types of spatial constraints are included in this research: ecological conservation redline areas; basic farmland conservation areas; and ecological safety network areas ( Figure A1). The ecological safety network is constructed based on landscape security pattern theory [72]. First, the Getis-Ord Gi* statistics are calculated by ArcGIS software version 10.7 to describe the spatial locations of cold spots and hot spots of ESVs in Fengdu County, and the hot spots with a confidence level of more than 90% are considered ecological sources and extracted [73]. Second, a minimum cumulative resistance model is used to construct an ecological resistance surface of Fengdu County, and Linkage Mapper software is used to calculate the distance path with the lowest cumulative consumption and, then, identify ecological corridors [74]. Finally, we analyzed buffer zones around the corridors and found that the corridor width does not include construction land within 300 m, and that the effects of human intervention are small. Considering that the extent to which land use types are excluded from ecological corridors differs by scenario, the widths of the ecological corridors in the SEC scenario, MEC scenario, and REC scenario were set to 300 m, 200 m, and 100 m, respectively [75]. Neighborhood factors refer to the difficulty of conversion between different land use types, with values ranging from 0 to 1: the closer the value is to 1, the stronger the expansion ability of the land use type is. The conversion cost is represented by a matrix, which is used to characterize the difficulty of converting from the current land use type to the demand type, where 1 represents a scenario in which the conversion between the two land types is allowed, and 0 represents a scenario in which there is no conversion. The neighborhood weight parameters (Tables A1-A5) and  conversion cost matrices (Tables A2-A5) for different scenarios included in this research were determined based on the results of related research [76,77]. Finally, the number of iterations was set to 1000, and the neighborhood size was set to 3 to simulate the spatial configuration of land use under the different ecological constraint scenarios.

Model Validation
To better simulate land use changes in the future, the figure of merit (FoM) coefficient and Kappa coefficient were introduced to verify the results of the FLUS model [78][79][80]. The larger the FoM coefficient is, the better the simulation effect and the higher the accuracyhowever, practical verification shows that FoM coefficients are mostly maintained within the range of 0-0.3, with values from 0.1 to 0.2 being the most common [80]. The closer the Kappa coefficient is to 1, the better the simulation accuracy. Kappa coefficients greater than 0.8 indicate that the simulation accuracy has reached an ideal state [79]. Based on land use data from 2010 and the FLUS model, the simulation result was obtained for 2018. Then, the simulation results for 2018 are compared with the observed results for 2018, and the results demonstrate that the overall accuracy of the model simulation is 93.20%, the FoM coefficient is 0.165, and the Kappa coefficient is 0.88. The validation results show that the FLUS model produces a highly accurate simulation of land use change in Fengdu County; therefore, this model can be applied for the following multiscenario simulation.

Results of the Land Use Type Area Optimization
The results of the land use type area optimization under the three scenarios were calculated by Lingo software version 12 ( Table 4). The results showed that there were large differences in the areas of various types of land use under the three scenarios. 1243.32 km 2 , and 53.58 km 2 , respectively. With the relaxation of the ecological constraints, the increases in the ISLs and the reductions in the CSLs were further reduced, whereas the increases in the SLs were further increased. The areas of forest and urban land increased to 1369.96 km 2 and 45.38 km 2 , and the areas of paddy fields and dry land decreased to 350.24 km 2 and 893.08 m 2 , respectively. The areas of rural residential sites, grassland, and water remained unchanged.

Optimization of the Spatial Configuration of Land Use
Based on the 2018 land use data, the corresponding quantitative structure, and spatial criteria, the spatial configuration of the land use types in the study area in 2035 under the three ecological constraint scenarios was optimized ( Figure 5). the source and sink landscapes was consistent, and the increases in the ISLs and the reductions in the CSLs were less than those in the SEC scenario, but the spatial change trend was similar to that in the SEC scenario. Among the land use types in the sink landscapes, whereas urban land continued to expand westward, there was also a trend of eastward expansion. The newly added urban land was mainly concentrated in the Mingshan subdistrict and Sanhe subdistrict, Shuanglu town, and Xingyi town. Moreover, a small number of rural settlements were transformed into forests, dry land, and paddy fields. Under this scenario, the MPFD of the study area was 1.16, the DIVISION was 0.96, the SHDI was 1.32, and the AI was 93.67. The configuration of the central town was relatively regular, the degree of fragmentation was moderate, and the spatial aggregation of different landscape types was moderate. (3) The development goal of the REC scenario is to maximize economic benefits. As shown in Figure 5c, under this scenario, the spatial scope of the sink landscapes was significantly expanded, the area occupied by sink landscapes exceeded that occupied by source landscapes, and the ES supply in the study area faced very large decreasing pres-  (1) The development goals of the SEC scenario are to ensure regional ecological security, reduce regional ecological problems through a series of ecological restoration projects, and comprehensively improve ES functions. As shown in Figure 5a, under this scenario, the spatial scope of the ISL expanded rapidly, and the area of the source landscapes was generally greater than that of the Sins-this balance provided support for improving the overall ES function of the study area. The advancement of the Grain for Green Program has transformed many dry land areas and paddy fields with large farming radii and steep slopes into forestland. The area covered by water also expanded significantly. The newly added water was mainly concentrated along the Yangtze River. The low-lying terrain of this area is conducive to the formation of water features. With the effective implementation of prevention and control measures for rocky desertification, all the unused land in the study area has been transformed into other land use types. Urban land mainly expanded to the west, but the expansion was not obvious. Most of the newly added urban land was transformed from dry land, and the new urban land area was mainly concentrated in the Mingshan subdistrict and Sanhe subdistrict. Moreover, some rural settlements located in remote mountainous areas, with underdeveloped transportation systems and extreme rates of population decrease, were transformed into forests, dry land, or paddy fields after the abandonment of the settlements. The mean patch fractal dimension (MPFD), landscape division index (DIVISION), Shannon's diversity index (SHDI), and aggregation index (AI) were calculated in Fragstats software version 4.2.
The results showed that the MPFD was 1.09, DIVISION was 0.94, SHDI was 1.58, and AI was 93.37. The configuration of the central town was relatively regular, the degree of fragmentation was low, and the spatial aggregation of different landscape types was relatively high. (2) The development goal of the MEC scenario is to consider both economic growth and ecological protection. As shown in Figure 5b, under this scenario, the development of the source and sink landscapes was consistent, and the increases in the ISLs and the reductions in the CSLs were less than those in the SEC scenario, but the spatial change trend was similar to that in the SEC scenario. Among the land use types in the sink landscapes, whereas urban land continued to expand westward, there was also a trend of eastward expansion. The newly added urban land was mainly concentrated in the Mingshan subdistrict and Sanhe subdistrict, Shuanglu town, and Xingyi town. Moreover, a small number of rural settlements were transformed into forests, dry land, and paddy fields. Under this scenario, the MPFD of the study area was 1.16, the DIVISION was 0.96, the SHDI was 1.32, and the AI was 93.67. The configuration of the central town was relatively regular, the degree of fragmentation was moderate, and the spatial aggregation of different landscape types was moderate. (3) The development goal of the REC scenario is to maximize economic benefits. As shown in Figure 5c, under this scenario, the spatial scope of the sink landscapes was significantly expanded, the area occupied by sink landscapes exceeded that occupied by source landscapes, and the ES supply in the study area faced very large decreasing pressure. Of the three scenarios, the increase in urban land use was highest in this scenario. In addition to Mingshan subdistrict, Sanhe subdistrict, Shuanglu town, and Xingyi town, urban land also increased significantly within Zhanpu town. The urban land area expanded from east to west along the Yangtze River. The main reason for this trend is that regions to the south and north are densely covered by basic farmland protection areas of existing towns, restricting the expansion of urban land. In addition, the western part of the county has flat terrain, dense road networks, excellent geographic locations, and functional water and land transportation facilities; therefore, it is suitable for urban development. The increases in the ISLs and the decreases in the CSLs in this scenario were further reduced compared with those in the SEC scenario, and the spatial change trend was similar to that in the SEC scenario.
In the REC scenario, the MPFD of the study area was 1.38, the DIVISION was 0.98, the SHDI was 1.23, and the AI was 92.08. The configuration of the central town was relatively scattered, the degree of fragmentation was high, and the degree of spatial aggregation of different landscape types was low.

Land Use Regionalization and Control Strategy in Fengdu County
We used ArcGIS software version 10.7 to superimpose the optimized spatial distributions of the source and sink landscapes under the three ecological constraint scenarios. With the objective of optimizing the ES function of the study area, the results of this superposition were amended, and Fengdu County was divided into seven types of areas ( Figure 6). The amendment rules are as follows: (1) Priority was given to ensuring ecological areas in Fengdu County and the grid cells that contained only ISL under all three ecological constraint scenarios were allocated to the ecological conservation area (ECA). Second, the grid cells that contained only CSLs and the grid cells that contained only SLs under the three scenarios were allocated to the agricultural production area (APA) and construction optimization area (COA). (2) The grid cells with different landscape types under the three scenarios were allocated to flexible control areas with multiple spatial functions, including the construction-ecological area (CEA), ecological-agricultural area (EAA), agriculturalconstruction area (ACA), and integrated development area (IDA). (3) Within the study area, the grid cells in the COA, CEA, and IDA that exceeded the boundaries of the permitted construction area were allocated to the ECA, and the grid cells that exceeded the boundaries of the permitted construction area in the ACA were allocated to the APA. (4) To ensure agricultural and ecological stability relative to the spatial configuration of the seven types of areas, a small number of COAs with heavy fragmentation were allocated to the ECAs or APAs based on their suitability. (2) The ECA was the most widely distributed type of area in Fengdu County, and the highest concentrations of the ECA were located in the area south of the Yangtze River. The ECA covered 1556.83 km 2 , accounting for approximately 53.67% of the total area of the study area. Urban and rural construction, and cultivated land reclamation activities, must be strictly prohibited in this area, and activities that destroy landscapes, vegetation, topography, and landforms, such as mining and borrowing are prohibited. Second, ecological restoration should be carried out in key areas that have been destroyed by human activities to optimize the vegetation coverage and the water environment. Moreover, it is possible to sustainably develop and utilize the ecological resources in this region, develop ecotourism, popularize science education and other projects, and enhance the cultural service values of the ecosystem. At the same time, it could be commented that better ecological corridors for the forest should be restored in the north of and along the Yangtze river. (3) The APA was widely distributed throughout the study area, and the highest concentration was in the area north of the Yangtze River. The area of the APA was 1123.23 km 2 , accounting for approximately 38.72% of the total area of the study area. First, it should be confirmed that the basic farmland in the study area is occupied and not abandoned. Second, the transformation of land should be encouraged to promote management of cultivated land at the appropriate scale. Moreover, to comprehensively improve the production environment and enhance agricultural modernization, high-quality farmland construction projects and rural land consolidation projects should be actively carried out, and appropriately planned ditches, mechanical tillage roads, breeding and seedling sites, greenhouses, hardening and drying yards, and other agricultural facilities should be allowed to occupy cultivated land. In addition, the demand for land associated with the development of agricultural science and technology research, ecological agricultural tourism, and other projects in the region should be supported. Moreover, it is necessary to strengthen the spatial matching between the permanent population and construction land, encourage the turnover of idle homesteads, and carry out cross-regional transactions and allocation of land use indicators to achieve a balance between the growth and decline in construction land between urban and rural areas. Second, the redevelopment of old urban areas and inefficient land should be actively carried out to develop emerging industries or green spaces. Furthermore, village renovation should be encouraged, and rural settlements should be guided to moderate concentration. The construction of transportation facilities, public service facilities, park squares, and other projects should also be carried out to comprehensively improve the living conditions of rural residents. (5) The CEA was mainly distributed around the COA, and the area of the CEA was only 2.60 km 2 . In the development and planning processes, the CEA should be prioritized for the protection and restoration of the ecological environment, and the development of eco-friendly spaces, such as country parks and ecologically friendly green spaces, should be encouraged. While adhering to the priority of ecological protection, the CEA can also be developed as construction land-however, the land in the area should not be transferred to high-polluting enterprises, such as the chemical and metal smelting industries. (6) The EAA was mainly distributed around the ECA, with an area of 163.56 km 2 , accounting for approximately 5.64% of the total area of the study area. In the future, this area should be given priority to be used as ecological space for the cultivation of forests. Moreover, it can be used for agricultural space to carry out production activities, such as grain and cash crop planting-however, green agriculture should be developed, and pollutants, such as pesticides and chemical fertilizers, should be completely prohibited. The EAA can also encourage the development of eco-agricultural tourism in accordance with local conditions to promote rural revitalization in Fengdu County. (7) The ACA, which was largely in the middle of the study area, was mainly distributed around the ECA. The area of the ACA was 21.96 km 2 , accounting for approximately 0.76% of the total area of the study area. This area can not only be used as agricultural land for the cultivation of agricultural products, but also can be developed to meet the land demand for economic development and residents' lives. Moreover, the ACA is close to urban and rural residential areas, giving it good geographic locations. Therefore, the cultivated land in the ACA should be planted mainly with high valueadded cash crops to improve the economic benefits of land use. (8) IDA was distributed mainly around the ECA, APA, and COA. This type of land had the smallest area in Fengdu County, with an area of only 1.17 km 2 . The IDA had the most diverse land use functions related to the process of regional development, and its use should be focused on the protection of ecological spaces and the establishment of a good foundation for the improvement of ES functions in the study area. After those objectives are achieved, then, this area can be used as construction space to increase the living area of residents, and as agricultural space to increase the amount of arable land. Therefore, a highly flexible land use control policy should be implemented in the IDA, allowing this land use type to provide space for various types of land according to the development needs of Fengdu County, and promoting the coordinated development of space in the territory.

Land Use Regionalization and Control Strategy in Fengdu County
We used ArcGIS software version 10.7 to superimpose the optimized spatial distributions of the source and sink landscapes under the three ecological constraint scenarios. With the objective of optimizing the ES function of the study area, the results of this superposition were amended, and Fengdu County was divided into seven types of areas ( Figure 6). The amendment rules are as follows: (1) Priority was given to ensuring ecological areas in Fengdu County and the grid cells that contained only ISL under all three ecological constraint scenarios were allocated to the ecological conservation area (ECA). Second, the grid cells that contained only CSLs and the grid cells that contained only SLs under the three scenarios were allocated to the agricultural production area (APA) and construction optimization area (COA). (2) The

Advantages of Land Use Optimization Models
On the basis of the observed situation in the study area and the source and sink landscapes associated with ESs, we proposed a method to optimize the areas of various land use types and their spatial configurations in Fengdu County by using a MOP-FLUS coupling model. When optimizing the areas of the land use types, we considered total area constraints, planning objective constraints, ecological environment constraints, and intensive development constraints, and established different objectives according to the varying emphasis on ecological and economic benefits in different scenarios. This optimization process is of great significance for improving the comprehensive benefits of space in a territory. We selected 18 driving factors of land use change, and 3 types of space-restricted areas to construct the FLUS model, thereby optimizing the spatial configuration of land use under different scenarios. The model validation results showed that the kappa coefficient was 0.88 and the FoM coefficient was 0.172, indicating that the model can effectively capture the trends in land use in Fengdu County, and that the results were more accurate than those of some studies that also use the FLUS model [64,69]. This finding may be because the land use data used in this study had a higher spatial resolution and fewer land use types than those used in other studies, and because the selection of driving factors was more comprehensive in this study than in other studies. In addition, spatial configuration optimization takes into account the areas occupied by ecological conservation redlines, basic farmland conservation areas, and ecological safety networks, which are of great significance for improving the overall ecological security of the region.
In summary, the method applied in this study addressed the problems that have arisen in recent land use optimization research, such as the lack of systematic constraints for determining the areas of various land use types, the lack of comprehensive spatial configuration optimization rules, and the unsatisfactory model accuracy, to a certain extent. Considering these issues, the predicted results had high credibility, and can provide valuable references for spatial planning in territories.

Feasibility of Land Use Regionalization
Due to the complexity of the natural environment, land use zoning of the TGRA is not only an effective way to plan and manage the land, but also a significant factor in the improvement of ES functions [64,81]. We divided Fengdu County into three areas with rigid control (the ECA, APA, and COA) and four areas with flexible control (the CEA, EAA, ACA, IDA). Among these areas, the COA was distributed mainly in the flat area along the Yangtze River, which is the core area for residential sites and economic development. The ECA was distributed mainly in high-altitude mountainous areas and along the main parts of the Yangtze River, which is the core area of ES supply. The APA was distributed mainly in the flat area between the two mountains in the study area. This is the main production area of grains in Fengdu County, and this area is of great significance for ensuring regional food security. The ECA and the ACA were distributed around the COA. These areas provide a buffer zone between the Sous and the Sins, which can effectively prevent the encroachment of important regional ecological spaces, and limit the pollution of soil, water, and air by various pollutants from construction land [82]. Moreover, the ECA and ACA can also alleviate the problem of overcrowded living spaces. The EAA, which establishes a natural barrier between agricultural space and ecological space, and can limit the damage of agricultural pollution to the ecological environment, was widely distributed around the APA. The IDA was distributed mainly around the COA, ECA, and ACA. The IDA has the most diverse spatial functions, and can provide space for the future development needs of Fengdu County.
In general, the land use regionalization plan proposed in this study is essentially a spatial functional regionalization plan. It does not only promote the preservation of values and appreciation of ESs, but also gives consideration to the coordinated development of urban areas, agriculture, and ecological functions in a space, which can provide the basis for delineating the "three zones and three lines" (three zones: ecological zone, agricultural zone, and urban zone; three lines: permanent basic farmland conservation redline, urban development boundary, and ecological conservation redline), and identifying key areas for ecological restoration in territorial spatial planning. Moreover, the land use regionalization plan reflects a moderately strict spatial management method for territories, which is in line with the observed status and future trends of the study area, and is highly feasible.

Implications for Spatial Management in the TGRA Territory
Since the United Nations issued the 2030 Sustainable Development Goals in 2015, territorial spatial governance has gradually become a hot topic in global change and sustainable development research [83]. The terrain of the TGRA is dominated by mountains and hills, and the water network is dense. Various ecological elements, including mountains, water, forestland, lakes, and grass, are integrated in the TGRA, and the natural background conditions are good. However, due to the very large natural and artificial disturbances associated with the construction of the Three Gorges Dam, and the rapid urbanization and industrialization that began when Chongqing became a municipality in 1997, the original ecosystem structure and functions of the TGRA have been seriously negatively affected, and a series of ecological and environmental problems have followed, such as forest vegetation degradation, water capacity decline, and frequent geological disasters [52]. Therefore, there is a need to strengthen the territorial spatial management in the region and build a strong ecological security barrier. Combining the source and sink landscape theory and the results of this research, we suggest that the spatial governance strategy of the TGRA should focus on the following two aspects: (1) At the macro level, in the TGRA, a more coordinated spatial development and protection plan should be developed, and the direct expansion of Sins into Sous should be slowed. First, a control policy related to the three zones and three lines must be strictly implemented in the TGRA, the system for delineating the main functional zones needs to be improved, and the ecological, agricultural, and urban spaces need to be rationally distributed. Moreover, the potential of the existing construction land needs to be fully exploited, and the renovation of rural residential sites needs to be vigorously carried out to continuously improve the compactness of the construction space. In addition, a market-oriented and diversified regional ecological compensation mechanism should be established to improve the enthusiasm and sustainability of ecological environmental protection.
(2) At the micro level, the TGRA should take ecological protection and restoration measures in response to the main problems faced by Sous. (a) The TRGA should continue to implement the Grain for Green project and Closing Hills for Afforestation project, especially in karst rocky desertification areas. Moreover, ecological restoration of mining areas should be implemented to repair damaged mountains, and continue to improve the vegetation coverage of the TGRA. In addition, the structure of natural forest vegetation should be optimized to prevent the formation of barren green deserts. (b) The prevention and management of urban and rural domestic sewage, industrial wastewater discharge and ship pollution, should be strengthened. Remediation projects in riparian zones should be implemented to prevent water-land cross-pollution. Comprehensive aquatic ecological restoration projects should also be launched to improve the quality of water sources and water biodiversity in reservoir areas. (c) It is necessary to optimize the configuration of agricultural lands, control agricultural nonpoint source pollution, and implement fallow fields on degraded cultivated land. Restoration projects should be implemented on sloping farmland in the TGRA to prevent further soil erosion.

Limitations
Due to the limited availability of basic data and the complexity of the model, the constraint conditions and optimization rules in this study are still relatively simple, and it is assumed that many driving factors of land use changes will continue as before, which may lead to results that are not robust or timely [63,64]. Therefore, although the results of this study provide a reference for regional spatial planning, they cannot provide an alternative proposal that can be implemented. In the future, more basic data and optimization rules should be added to the model, and the uncertainties in the driving factors of land use change should be considered more comprehensively and carefully to improve the accuracy of the land use optimization.
The demand for ESs is driven by the consumption of, or the hope to obtain, various services that support human society, and this demand is a potential determinant of optimal land use allocation [84]. As the strengthening of ES demand will inevitably cause stakeholders to balance, adjust, and control land use patterns, this will affect the direction, objectives, and planning of territorial spatial development. However, this study considered only the ES supply and did not consider the ES demand. In the future, the impacts of ES demand on territorial space should be fully considered, and the land use pattern should be optimized more accurately on the basis of the ES supply and demand relationship, spatial differentiation characteristics, and degree of matching.

Conclusions
(1) By integrating the MOP model and FLUS model, a study of county-level land use optimization was carried out according to the following flow: optimization of the areas of various land use types, followed by spatial configuration optimization and land use regionalization. This flow can effectively enhance the reliability and applicability of land use optimization results. (2) The results of the land use optimization process showed that there were large differences in the areas and spatial distributions of the source and sink landscapes under the three ecological constraint scenarios. Under the SEC scenario, the areas of the ISLs, CSLs, and SLs were 1676.62 km 2 , 1190.43 km 2 , and 33.81 km 2 , respectively. A large area of CSLs and a small area of sink landscapes were transformed into ISLs, the degree of landscape fragmentation was low, and the area covered by source landscapes was generally greater than that covered by sink landscapes. Under the MEC scenario, the areas of the ISLs, CSLs, and SLs were 1609.22 km 2 , 1241.60 km 2 , and 49.74 km 2 , respectively. The development of the source and sink landscapes was consistent, the degree of fragmentation was moderate, and the land use changes ranged from centralized development to balanced development. Under the MEC scenario, the areas of the ISLs, CSLs, and SLs were 1603.96 km 2 , 1243.32 km 2 , and 53.58 km 2 , respectively. A large amount of CSLs were transformed into sink landscapes, the degree of fragmentation was high, and the area occupied by Sins was generally higher than that occupied by source landscapes. (3) The results of land use optimization under different scenarios were superimposed, and the results indicated that Fengdu County is divided into seven types of areas: ecological conservation area; agricultural production area; construction optimization area; construction-ecological area; ecological-agricultural area; agricultural-construction area; and integrated development area. Different areas have different spatial functions, and differentiated land management policies should be adopted. (4) In the future, the spatial governance of all the counties in the TGRA territory should be strengthened, more coordinated land space development and protection patterns should be developed, and ecological protection and restoration projects should be carried out in mountains, rivers, forests, fields, lakes, and grasslands to comprehensively improve regional ES functions.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The abbreviations in this article:  Figure A1. Spatial constraints area.  Figure A1. Spatial constraints area.