Development of a Composite Model for Simulating Landscape Pattern Optimization Allocation: A Case Study in the Longquanyi District of Chengdu City, Sichuan Province, China

Abstract

The simulation of landscape pattern optimization allocation (LPOA) to achieve ecological security is an important issue when constructing regional ecological security patterns. In this study, an LPOA model was developed by integrating a binary logistic regression model and a nonlinear programming model with a particle swarm optimization algorithm in order to consider the complexity of landscape pattern optimization in terms of the quantitative structure and spatial layout optimization, integrating the landscape suitability and factors that influence landscape patterns, and under constraints to maximize the economic, ecological, and comprehensive benefits of landscape patterns. The model was employed to simulate the LPOA in the Longquanyi District of Chengdu City, Sichuan Province, China. The model successfully obtained an appropriate combination of the landscape quantitative structure and spatial layout, as well as effectively integrating the landscape suitability and factors that influence the landscape pattern. Thus, the model addressed the problems of previous studies, such as neglecting the coupling between quantitative structure optimization and spatial layout optimization, ignoring the macrofactors that affect landscape patterns during optimization modeling, and initializing particles without considering the suitability of the landscape. Furthermore, we assessed and analyzed the accuracy and feasibility of the landscape pattern spatial layout optimization results, where the results showed that the overall accuracy of the optimization results was 84.98% with a Kappa coefficient of 0.7587, thereby indicating the good performance of the model. Moreover, the simulated optimization allocation scheme for the landscape pattern was consistent with the actual situation. Therefore, this model can provide support and a scientific basis for regional landscape pattern planning, land use planning, urban planning, and other related spatial planning processes.

1. Introduction

During the rapid economic development process worldwide, the lack of effective regulation for urban expansion has led to a disorderly urban sprawl, where the blind occupation of ecological spaces and severe destruction of the natural ecological environment have caused many global or regional ecological environmental problems, such as global warming, ozone layer depletion, loss of biodiversity, spread of acid rain, sharp decline of forests, land desertification, air pollution, water pollution, and marine pollution. These problems have severely restricted the sustainable development of human society. Since the start of economic reform in 1978, China has achieved remarkable economic progress. At present, China is the second largest economy and it is among the fastest growing economies in the world. However, similar to other countries in the midst of rapid economic development, disorderly urban sprawl, rapid industrial development, and inappropriate farmland use in China have led to problems, including the blind development of land resources and casual occupation of ecological space. These problems have resulted in inappropriate land use and poor ecological protection, as well as adversely affecting the structure and function of terrestrial ecosystems, thereby resulting in environmental issues that restrict sustainable development practices in China, such as farmland loss and land degradation, increased pollution with waste from industry, severe soil pollution of farmland, shortages of water and soil resources, loss of biodiversity, and fragmentation of the countryside landscape [1,2,3,4,5]. The local ecological environment has been improved by controlling the environment in China with tough measures, such as central government environmental supervision and ecological conservation redline planning, but the overall trend toward deterioration has not been halted. This is mainly because the protection of natural resources has not been planned scientifically from a spatial perspective. Thus, the current ecological protection processes were often implemented blindly and in an inefficient manner. Moreover, the conflict between the requirements for economic development and ecological protection has intensified [6,7,8,9]. Thus, coordinating the contradictory relationship between economic development and ecological protection in order to maximize the value of natural resources without affecting ecological protection has become a key issue that affects the sustainable development of China. Many studies have confirmed that constructing a reasonable landscape (land resource) spatial pattern by optimizing the landscape (land resource) allocation from a spatial perspective is effective for mitigating conflicts between economic development and ecological protection, and to achieve the goal of sustainable development [6,9,10,11,12,13,14,15,16].

Landscape pattern optimization allocation (LPOA) is a complex and challenging spatial resource optimization allocation problem, which involves resolving conflicts of interest between economic development and ecological protection, considering attribute characteristics (e.g., ecological, environmental, economic, and cultural factors) and spatial characteristics (e.g., morphology and compactness), and solving multiobjective optimization problems and complex models and algorithms [7,17]. In order to achieve LPOA, researchers in China and other countries have launched extensive studies and initiatives regarding this problem. Four main types of optimization models have been applied comprising quantitative structure optimization, spatial evolution simulation, spatial layout optimization, and composite optimization models.

The methods used for optimizing the quantitative structures of landscape patterns mainly comprise mathematical programming (MP)-based models (e.g., linear programming, nonlinear programming (NP), goal programming, integer programming, and uncertain programming), system dynamics (SD)-based models, and heuristic algorithm (HA)-based models, which have been employed widely for optimizing land use in terms of the quantitative structure and land area [16,18,19,20,21,22,23]. MP-based models can rapidly determine the optimal land use structure according to specific objectives and constraints, but they cannot change the land use for parcels and allow spatial optimization [24]. However, compared with other types of models, MP-based models can be generated relatively conveniently from the optimization toolbox in software such as MATLAB, LINGO, and LINDO. Therefore, MP methods still have many possible applications in the optimization of the quantitative structure of landscape patterns. SD-based models can explain the driving factors and the trends in land use changes in a region fairly well, as well as dynamically allocating the future land use quantitative structure in a region [23]. Thus, SD-based models mainly focus on dynamic simulations of the components of the system as well as the relationships among the components, and they are generally used to simulate how land use demands are influenced by the economy, technology, population, policy, and their interactions at macroscales [25]. Hence, it is impossible to optimize the spatial layout of the landscape pattern directly with these models.

Due to advances in computers and GIS technology, numerous models have been developed such as cellular automata [26,27], the conversion of land use and its effect on small regional extent (CLUE-S) [28,29], agent-based models [30,31], scenario analysis [32,33], future land use simulation model [27,34], and hybrid simulation models [35,36,37], which can be employed for simulating the spatial dynamics of changes in landscape patterns. However, these simulation models aim to predict future landscape patterns or land use, rather than optimizing the spatial layout of landscape patterns or land use. Moreover, these models are limited to optimizing landscape patterns or land use types to generate only a few landscape patterns or land use optimization allocation schemes. Heuristic methods can generate many more landscape pattern (land use) optimization allocation schemes to search for a better solution [14]. Therefore, heuristic algorithms such as genetic algorithm (GA), simulated annealing algorithm, particle swarm optimization (PSO), artificial immune system algorithm, artificial bee colony algorithm, ant colony optimization, and tabu search methods have been employed frequently for land use spatial layout optimization with the support of GIS technology [17,19,38,39,40,41,42], whereas they have been applied rarely to landscape pattern spatial layout optimization (LPSLO). Compared with other heuristic algorithms, the main advantages of PSO algorithms are the flexibility and simplicity of its operators, improved space search capability and adaptability, and rapid convergence rate. Thus, PSO algorithms have been widely applied to solve the land-use spatial allocation optimization problem [9,40,43], and they are effective tools for solving complex spatial layout optimization problems for landscape patterns [44]. However, these heuristic-based models mainly focus on spatial pattern optimization whereas they ignore the quantitative structures and beneficial optimization of land use.

The aforementioned models, including quantitative structure optimization, spatial evolution simulation, and spatial layout optimization models, have made positive contributions to studies of LPOA and land-use optimization allocation, but they only have specific advantages. For example, MP-based models are good at optimizing quantitative structures, SD-based models are good at integrating macrofactors in resource allocation processes, and HA-based models are good at optimizing spatial layouts, and thus these models have difficulty with LPOA, which is a complicated spatial optimization decision-making problem. To address these challenges, various composite models for making spatial optimization decisions have been applied to solve the problems of optimizing the landscape and land resource spatial allocation because composite models are capable of integrating the advantages of various models, including quantitative structure optimization, spatial evolution simulation, and spatial layout optimization models. These composite optimization models include a loosely coupled model based on GA and game theory [41], spatially explicit genetic algorithm that integrates land use planning knowledge with GA [45], differential evolution-cellular automata model [46], multiobjective land use optimization allocation model that integrates a multiagent system with PSO [47], mathematical-spatial optimum utilization model using fuzzy goal programming and multiobjective land allocation [48], and a method for supporting land use planning by combining the GA method, CLUE-S model, and water assessment tool model [49]. In recent years, some composite optimization models that integrate PSO algorithms with other optimization methods have been applied widely in land resource spatial optimization configuration research. For example, Liu et al. [6] presented a PSO model combined with multiobjective optimization techniques for the spatial optimization of rural land-use allocation. They first obtained the optimal land-use quantitative structure by linear programming, and then conducted land-use spatial layout optimization according to the initial particles generating by the optimal land type area. Liu et al. [25] proposed a novel model that integrated SD and hybrid PSO for solving land-use allocation problems in a large area. They first used the SD module to project land use demands under various scenarios, and further modified the PSO by incorporating genetic operators to allocate land use in various scenarios. The composite optimization model is an effective method for solving the problem of complicated spatial optimization decisions, and it provides a reference for the application of LPOA.

However, the aforementioned models for making spatial optimization decisions are still affected by problems that need to be solved urgently. First, most of the optimization models and methods employed for making spatial optimization decisions are focused on the quantitative structure or the spatial layout, whereas few consider both. They usually perform spatial layout optimization according to the land demand areas, but most ignore the coupling between quantitative structure optimization and spatial layout optimization. Second, some models and methods have failed to integrate macrofactors that influence landscape (land use) patterns, such as social, economic, ecological, policy, and institutional factor, thereby resulting in landscape pattern (land use) optimization allocation schemes that are obviously not appropriate for real situations. Third, when making spatial optimization decisions by integrating PSO algorithms with other optimization methods, most studies employed random functions to generate the initial particles, while a few studies initialized the particles according to the land requirement areas, whereas few studies generated the initial particles based on the actual suitability of the landscape or land use spatial layouts, thereby affecting the rationality of the optimization results to some extent.

To address the problems mentioned above, the composite model developed in this study aimed to simultaneously optimize the quantitative structure and spatial layout of the landscape pattern, as well as effectively integrating the landscape suitability and macrofactors that affect landscape pattern. Therefore, in terms of the optimization approaches employed, this composite model mainly involved landscape suitability evaluation methods, landscape pattern quantitative structure optimization (LPQSO) methods, and LPSLO methods. Early efforts at land-use spatial optimization mainly focused on allocating the most feasible land with the highest suitability to a specified land use unit [50,51]. Thus, landscape suitability evaluation was an important basis for LPOA. Binary logistic regression (BLR) models are highly appropriate for landscape suitability evaluations (LSEs) because each landscape type can have two statuses, i.e., “present” and “not present,” within a certain spatial range. Optimizing the quantitative structure of landscape pattern is a complex constrained optimization problem because this problem includes many constraints, such as land use planning and the ecological environment. NP is an effective method for solving constrained optimization problems and the underlying theory of the algorithms employed is mature (e.g., reduced gradient method and penalty function method), while they can be readily implemented using MATLAB. Therefore, there have been many successful applications of NP to practical quantitative structure optimization [45,52,53,54]. In addition, PSO is a stochastic optimization algorithm based on swarm intelligence, which guides particles to determine the optimum feasibility region in complex search spaces by simulating the social behavior of bird flocking and fish schooling. Compared with other traditional evolutionary algorithms, PSO has a superior capacity for searching the space as well as adaptability, in addition to a rapid convergence rate [25,43,47,55]. Furthermore, the successful applications of PSO to various problems have demonstrated its potential, such as the spatial optimization of land-use allocation [6,9,25,40,43,47], facilities location selection [56], optimal allocation of earthquake emergency shelters [57], model parameter optimization [58], feature selection in classification [59], and other problems [60,61]. Many studies have demonstrated that PSO is highly robust and it can obtain more different routes through the problem hyperspace than other evolutionary algorithms [62]. Therefore, a composite optimization method that integrates BLR and NP with PSO for LPOA may be an effective approach for addressing the problems caused by neglecting the coupling between quantitative structure optimization and spatial layout optimization, ignoring the macrofactors that affect landscape patterns when optimizing modeling, and initializing particles without considering the suitability of the landscape, as well as enhancing the practical utility of the LPOA results. Thus, in this study, we developed an LPOA composite model by integrating BLR and NP with PSO, which we employed to simulate the LPOA for Longquanyi District of Chengdu City, Sichuan Province, China. The results of this study may provide a useful reference for formulating landscape pattern planning, land use planning, urban planning, and other related spatial planning, as well as providing a valuable basis for implementing similar studies in other areas.

2. Materials and Methods

2.1. Study Area and Data

2.1.1. Study Area

Longquanyi District is the main eastern urban development area in Chengdu City and the location of the National Economic and Technological Development Zone. We selected this district as the study area (Figure 1) because of its high representativeness and research value for LPOA. This area has also been affected by series of ecological and environmental problems during the process of industrialization and urbanization. The area is located in the central section of Sichuan Province, to the east edge of Chengdu Plain and the west side of Longquan Mountain. Geographically, this district is located between 104°08′19″ E–104°27′09″ E and 30°27′52″ N–30°43′23″ N, at 12 km from the center of Chengdu City and it has the typical characteristics of a suburban area in a metropolis. The total area of the entire district is 556.4 km². The geological structure comprises a structural fault block between the Chengdu fault depression belt and the Longquan Mountain uplift fold belt. The characteristic terrain is high in the southeast but low in the northwest. The highest elevation is 1037 m, the lowest elevation is 407 m, and the relative elevation is 630 m. The geomorphic types mainly include plains, hills, and mountains, where the plains are piedmont alluvial dam distributed in the mid-west of the district, mountains are present with a northeast-southwest alignment in the southeast of the district, and the hills are distributed in the west of the district to the east and west sides of the middle section of Longquan Mountain. The regional climate type is subtropical humid in the Sichuan Basin. The mean annual temperature = 16.5 °C, precipitation = 852.4 mm, sunshine duration = 1021 h, and evaporation capacity = 984.7 mm. The soil types mainly comprise rice soil, yellow soil, purple alluvial soil, and purple soil. The forest vegetation type is subtropical evergreen broad-leaved forest and the major vegetation is natural secondary forest and plantation.

2.1.2. Data Sources

The original raw data included GIS data, remote sensing images, socioeconomic statistics, environmental data, as summarized in

Table 1. Data used for landscape pattern optimization allocation.

Types	Data	Source
GIS data	GDEM V2	National Aeronautics and Space Administration (NASA)
	Current land use status map	Local Bureau of Land and Resources
Remote sensing images	Landsat OLI image	United States Geological Survey (USGS)
Socioeconomic statistics	Total population, permanent populations in city and countryside, fruit yields, total production by major agriculture, gross products of primary, secondary, and tertiary industries, GDP	Local Bureau of Statistics
Environmental data	Soil organic matter content	Local Bureau of Rural Development and Forestry
	Annual rainfall, annual temperature	Local Bureau of Meteorological
Other data	Interpretation of remote sensing images, amounts of fertilizer applied to the main fruits and crops	Field survey
	General plan of Chengdu National Economic-Technological Development Zone	Local Bureau of Planning and Administration

. The GIS data included Global Digital Elevation Model Version 2 (GDEM V2) at 30 m resolution and the current land-use status map (1:100,000) for the year 2009. The remote sensing image comprised a Landsat Operational Land Imager (OLI) image taken in August 2014, which included eight multispectral bands (with a spatial resolution of 30 m) and one panchromatic band (with a spatial resolution of 15 m) with the WGS84 coordinate system (WGS84 ellipsoid and datum, universal transverse Mercator projection, and central meridian 105° E) and the world-wide reference system (path/row, 129/039). The socioeconomic statistics included the total population for the years 1978–2014, the permanent populations in the city and countryside during 2001–2014, fruit yields in 2000–2012, total production by major agriculture (agriculture, forestry, animal husbandry, side-lines, and fisheries) in 2014, the gross products of the primary, secondary, and tertiary industries, and the gross domestic product (GDP) in 2000–2014. The environmental data included the soil organic matter contents of 1842 soil samples in 2013, and annual rainfall and annual temperature data from nine automatic meteorological observation stations in Longquan District during 2004–2014. In addition, we used other data related to field survey results for landscape types, field survey data for the amounts of fertilizer applied to the main fruits and crops, and the general plan of Chengdu National Economic and Technological Development Zone.

2.1.3. Data Processing

In this study, the data imported into the model had two formats: spatial data and nonspatial data. During spatial data processing, the landscape type was classified into six categories (farmland, orchard, forest, urban-rural residential and industrial-mining, and waters) using the quick unbiased efficient statistical tree (QUEST) classification based on the Landsat OLI image from 2014, GDEM V2 data, the current land use status map, and the landscape field survey results, where the overall classification accuracy was 95.94% and it was sufficient for the purposes of this study [63]. We then obtained the landscape map of the study area (Figure 2a) and landscape maps of the priority planning area (Figure 2b) in the study area. Furthermore, derivative data such as the slope, aspect, and hypsography degree were obtained based on calculations using GDEM V2 data. Environmental data related to the soil organic matter content, average annual rainfall, and average annual temperature were obtained by spatial interpolation using ESRI ArcGIS 10.0. Neighborhood factors including the nearest distances to the city center, town center, major roads, and waters were calculated using a Euclidean distance tool based on related data such as the city center, administrative center of each town, major roads, and waters derived from the landscape map in 2014 with ESRI ArcGIS. In addition, socioeconomic data such as the population density and per capita GDP were calculated according to relevant sources, including the administrative district area, population, and GDP statistical data, before they were converted into grid data with the vector-raster conversion function in ESRI ArcGIS. In order to implement the models, all of the spatial data were converted into grid data format at 30 m resolution using the Xi’an 1980 coordinate system and the Gauss–Kruger projection. During nonspatial data processing, nonspatial data were collected for the study area, including natural status data and socioeconomic statistical data, which were used to establish the NP models to conduct the LPQSO.

2.2. LPOA Model

The LPOA model comprises three parts: LSE, LPQSO, and LPSLO. The implementation of the LPOA model can be divided into three parts, as follows. First, based on the values of the LSE indexes, the BLR model is used to conduct the LSEs. Second, according to the actual regional situation, where the conflict between economic development and ecological protection has become increasingly intense, we consider three scenarios comprising economic development, ecological protection, and overall considerations. Next, the NP model was used to conduct the LPQSO by selecting the maximum benefit in each development scenario as the optimization objective, as well as using each landscape type area as a decision variable, and the landscape area, ecological service value, nonpoint source pollution, and industrial structure as constraints. Finally, based on the spatial distribution of the landscape suitability obtained in the first step, optimal area of each landscape type obtained in the second step, landscape type conversion rules, and landscape type raster graphics in the base year for particle initialization, as well as the LPQSO results produced in the second step as the quantitative constraint conditions for the particle fitness function, an LPSLO model was established according to the principle of the PSO algorithm in order to optimize the landscape pattern spatial layout for each scenario. The framework of the LPOA model is shown as Figure 3.

2.2.1. LSE Model Based on BLR

An LSE measures the adaptability category for a certain landscape spatial distribution and it forms the basis of LPOA. The BLR model is highly appropriate for conducting an LSE because each landscape type has two possible statuses, i.e., “present” and “not present,” within a certain spatial range. Given that the probability of an event occurring in the BLR model represents the landscape suitability degree, then in the two-dimensional grids of $M\times N$ landscape pattern spatial distribution cells, a greater probability value of some landscape type appearing in a grid domain indicates that it is more suitable for allocating this landscape type. Thus, a higher probability for a landscape type represents a higher degree of suitability for the landscape spatial distribution. Hence, the LSE model based on BLR can be expressed as

$$p_{\mathsf{\theta }}=\frac{\exp (a+{\beta }_1{X}_1+{\beta }_2{X}_2+\dots +{\beta }_{\mathsf{\theta }}{X}_{\mathsf{\theta }})}{1+\exp (a+{\beta }_1{X}_1+{\beta }_2{X}_2+\dots +{\beta }_{\mathsf{\theta }}{X}_{\mathsf{\theta }})}$$

(1)

where $p_{\mathsf{\theta }}$ within the [0, 1] region is the probability of some landscape type appearing in a grid domain in the study area, i.e., it denotes the degree of suitability for a landscape type spatial distribution in a grid domain in the study area; $X_{\kappa }(\kappa =1,2,\dots ,\theta )$ is a factor that influences the landscape spatial distribution; $a$ is the regression equation constant; and ${\beta }_{\kappa }(\kappa =1,2,\dots ,\theta )$ is the regression coefficient. In the present study, the BLR model was implemented by Python programming using the ESRI ArcGIS platform. The interpretability of this model was checked using the receiver operating characteristic (ROC) method according to Pontius and Schneider [64].

2.2.2. LPQSO Model Based on NP

We established an LPQSO model to optimize the landscape quantitative structure under different scenarios in the target years of 2021 and 2028 by using the year 2014 as the base year, and then considering the areas with landscape types including farmland, orchard, forest, urban-rural residential and industrial-mining, and waters as the decision variables, as well as the landscape area, ecological service value, nonpoint source pollution, and industrial structure as constraints, and by selecting the maximum economic benefit, ecological security degree, and comprehensive benefit as optimization objectives.

Objective Functions

The objective functions in the LPQSO model are as follows.

(1) Objective functions in the economic development scenario: Producing more products and providing more services through the rational use of limited landscape resources were the main goals of economic development scenario. Therefore, in this study, the maximum summed economic outputs from all landscape types was selected in the model as the optimization objective for the economic development scenario, i.e.,

$$\mathrm{Maximize}: f(z)={\displaystyle \sum _{k=1}^K{c}_k{z}_k}$$

(2)

where $f(z)$ is the summed economic gross product value for all landscapes types, $c_k$ is the product value coefficient for landscape type $k$, $z_k$ is the area of landscape type $k$, and $K$ is the number of landscape types.

(2) Objective functions in the ecological protection scenario: Obtaining the maximum ecological security degree through the rational distribution of landscape resources was the main goal of ecological protection scenario. Therefore, in this study, the maximum summed ecological security indexes for all landscapes types was selected in the model as the optimization objective for the ecological protection scenario, i.e.,

$$\mathrm{Maximize}: g(z)={\displaystyle \sum _{k=1}^K{a}_k{z}_k}$$

(3)

where $g(z)$ is the summed ecological security indexes for all landscapes types, $a_k$ is the ecological security coefficient of landscape type $k$, and $z_k$ has the meaning stated above.

(3) Objective functions in the overall consideration scenario: Overall, the aim was to ensure the security of the ecological environment and to promote the continuous steady growth of the economy. Obtaining the maximum comprehensive benefit based on the overall layout of various landscape resources in the study area was the main goal of the overall consideration scenario. Thus, the maximum summed economic outputs and ecological security indexes for all landscape types were selected in the model as the optimization objectives for the overall consideration scenario, i.e.,

$$\mathrm{Maximize}: F(z)={\omega }_1{\displaystyle \sum _{k=1}^K{v}_k{z}_k}+{\omega }_2{\displaystyle \sum _{k=1}^K{v}_k^{\prime }{z}_k}$$

(4)

where $F(z)$ are the comprehensive benefits of various landscapes in the study area; ${\omega }_1$ and ${\omega }_2$ are the weights of the objective functions for the economic development scenario and ecological protection scenario, respectively; $v_k$ and $v_k^{\prime }$ are the standardized product value coefficient and standardized ecological security coefficient for landscape type $k$, respectively; and $z_k$ has the meaning stated above.

Constraints

The constraints in the LPQSO model are as follows.

(1) Landscape area: In this study, there were six landscape area constraints, i.e., total landscape area, farmland area, orchard area, forest area, urban-rural residential and industrial-mining area, and waters area. These constraint conditions can be expressed as:

$${\displaystyle \sum _{k=1}^K{z}_k=A}$$

(5)

where $A$ is the total area of each landscape type and $z_k$ and $K$ have the meanings stated above;

$$A_L\le z_1<A$$

(6)

where $A_L$ is the minimum demand area for farmland in the target year, $z_1$ is the optimal area of farmland in the target year, and $A$ has the meaning stated above;

$$A_G\le z_2<A$$

(7)

where $A_G$ is the minimum demand area for orchard in the target year, $z_2$ is the optimal area of orchard in the target year, and $A$ has the meaning stated above;

$$A_S\le z_3<A$$

(8)

where $A_S$ is the minimum demand area for forest in the target year, $z_3$ is the optimal area of forest in the target year, and $A$ has the meaning stated above;

$$A_C\le z_4<m_1{p}_1+m_2{p}_2+\phi +\delta$$

(9)

where $A_C$ is the minimum demand area for urban-rural residential and industrial-mining in the target year; $z_4$ is the optimal area of urban-rural residential and industrial-mining in the target year; $m_1$ and $m_2$ are the highest per capita land use standard in urban and countryside areas, respectively; $p_1$ and $p_2$ are the urban and rural resident populations in the target year, respectively; $\phi$ is the planning area for industrial land; and $\delta$ is the maximum planning area for transportation land;

$$A_W\le z_5<A$$

(10)

where $A_W$ is the minimum demand area for waters in the target year, $z_5$ is the optimal area of waters in the target year, and $A$ has the meaning stated above.

(2) Ecological service value: The ecological service value refers to the economic value of ecological services provided by landscape functions. This value can be estimated according to the ecosystem ecological services value equivalence factor [65], the grain yield per unit area, and the market price. In this study, the ecological service value is expressed as:

$${\xi }_{\min }\le {\displaystyle \sum _{k=1}^K{\epsilon }_k{z}_k<{\xi }_{\max }}$$

(11)

where ${\xi }_{\min }$ and ${\xi }_{\max }$ are the lowest and highest ecological service values in the study area, respectively; ${\epsilon }_k$ is the ecological service value per unit area for landscape type $k$; and $z_k$ has the meaning stated above.

(3) Nonpoint source pollution: In this study, we considered three constraints on nonpoint source pollution, i.e., annual load of chemical oxygen demand (COD), annual load of total nitrogen (TN), and annual load of total phosphorus (TP), which can be expressed as:

$$CO{D}_{\min }\le {10}^4{\displaystyle \sum _{k=1}^K{u}_{k}hL{C}_k{z}_k}<CO{D}_{\max }$$

(12)

where $CO{D}_{\min }$ and $CO{D}_{\max }$ are the minimum and maximum annual loads of COD in the target year, respectively; $u_k$ is the runoff coefficient for landscape type $k$; $h$ is the average rainfall; $L{C}_k$ is the COD concentration in surface runoff for landscape type $k$; and $z_k$ has the meaning stated above;

$$T{N}_{\min }\le {10}^4{\displaystyle \sum _{k=1}^K{u}_{k}hL{N}_k{z}_k}<T{N}_{\max }$$

(13)

where $T{N}_{\min }$ and $T{N}_{\max }$ are the minimum and maximum annual loads of TN in the target year, respectively; $L{N}_k$ is the TN concentration in surface runoff for landscape type $k$; $h$ and $z_k$ have the meanings stated above;

$$T{P}_{\min }\le {10}^4{\displaystyle \sum _{k=1}^K{u}_{k}hL{P}_k{z}_k}<T{P}_{\max }$$

(14)

where $T{P}_{\min }$ and $T{P}_{\max }$ are the minimum and maximum annual loads of TP in the target year, respectively; $L{P}_k$ is the TP concentration in surface runoff for landscape type $k$; and $u_k$, $h$, and $z_k$ have the meanings stated above.

(4) Industrial structure: In this study, the industrial structure is expressed as:

$${\zeta }_{\min }\le \frac{c_4{z}_4}{c_1{z}_1+c_2{z}_2+c_3{z}_3+c_5{z}_5}\le {\zeta }_{\max }$$

(15)

where ${\zeta }_{\min }$ and ${\zeta }_{\max }$ are the minimum and maximum secondary industry product value divided by the summed product values for primary industry and service industry, respectively; and $c_{k(k=1,2,\dots ,5)}$ and $z_{k(k=1,2,\dots ,5)}$ have the meanings stated above.

LPQSO Model Solution

In this study, the LPQSO models for the economic development, ecological protection, and overall consideration scenarios were implemented by MATLAB programming using the solution function FMINCON for the constraint optimization problem in the MATLAB optimization toolbox. We then determined the optimal areas for various landscape types under each scenario in the target years.

2.2.3. LPSLO Model Based on the PSO Algorithm

Spatial Optimization Principle for the LPSLO Model

Adjusting the position of a pixel according to a specific optimization objective is the basis of spatial layout optimization using landscape type raster data. Therefore, simulating the spatial distribution of pixels in a landscape pattern raster map using particle positions is important for conducting LPSLO with the PSO evolution algorithm. The raster map can be regarded as a real matrix, and thus its pixels correspond to the matrix elements, the pixel locations correspond to the row and column numbers of the matrix elements, and the pixel attribute values denoting landscape type codes correspond to the values of the matrix elements. Hence, adjusting the pixel locations and attribute values in the raster map is equivalent to adjusting the element row and column numbers and values in the matrix. Therefore, we assumed that the matrix ${\mathit{A}}_{m\times n}$ denoted the landscape type raster map and the matrix element value denoted the landscape type code, while the matrix denoted particles, matrix element values denote particle elements, and the row and column numbers of the matrix elements denote the locations of the particle elements. Thus, according to the principle of the PSO algorithm [66], regardless of how much the particle element space position changes, the particle element itself, i.e., the landscape type code, was not changed, and thus the element values comprising the matrix never changed. Therefore, a matrix element moved from one location to another by optimizing its row and column numbers using the PSO algorithm. Subsequently, a new matrix was obtained containing the optimized matrix element. Thus, the LPSLO was obtained successfully when the new matrix corresponding to the landscape pattern spatial layout obtained the maximum optimization objective. The spatial optimization principle for the LPSLO model is shown in Figure 4.

Key Modeling Techniques of the LPSLO Model

Considering that MATLAB is an advanced programming language where its basic programming unit is matrix, and that it has a powerful capacity for calculating matrices and image processing [67], we designed a PSO model and algorithm for implementing the LPSLO using MATLAB. We employed the following four key modeling techniques.

(1) Spatial mapping of the relationship between landscape type raster graphics and particles: We assume that $B$ is a vector quantity $({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_N)$ for the valid elements of matrix ${\mathit{A}}_{m\times n}$ corresponding to the landscape type raster graphics in the base year. A valid matrix element is the landscape type code, which has a corresponding grid value, but we exclude an element with a null value or one that is not within the range of encoding values for the landscape type. The vector quantity $B$ denotes a particle and it corresponds to a landscape pattern spatial layout scheme. $\alpha \in B$ is the value of a valid element and it denotes a particle element, where its row and column number also denote the position of the particle element. Therefore, the spatial mapping relationship between landscape type raster graphics and particles during the model optimization process is summarized in Figure 4.

(2) Optimization objective and constraints in the LPSLO model: Considering that the LPSLO is a multiobjective optimization problem, the maximum landscape suitability and spatial aggregation degree were selected as the optimization objective functions in the LPSLO model. The summed products for various landscape type spatial suitability values and their weights were selected to denote the landscape suitability, and the patch shape index was employed to denote the landscape spatial aggregation degree. In addition, three types of constraints were considered: landscape area, landscape type, and landscape type conversion rules. The optimization objective function and constraints can be expressed as

$$\begin{array}{cc}\mathrm{Maximize}:& F(\mathit{A})=w_1{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{l=1}^Q{\lambda }_k{p}_{kl}}}+w_2{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{r=1}^R{c}_{kr}/\sqrt{s_{kr}}}}\\ \mathrm{Subject} \mathrm{to}:& \{\begin{array}{c}{d}_k(\mathit{A})=D_k^{*},{\displaystyle \sum _{k=1}^K{d}_k(\mathit{A})}={\displaystyle \sum _{k=1}^K{D}_k^{*}}\\ {\alpha }_j\in [1,2,\dots ,K]\\ 0\le T_j(k\to k^{\prime })\le 1\end{array}\end{array}$$

(16)

where ${\lambda }_k$ is the weight of landscape type $k$; $p_{kl}$ is the spatial suitability value for pixel $l$ of landscape type $k$; $c_{kr}$ and $s_{kr}$ are the perimeter and area for patch $r$ of landscape type $k$, respectively; $w_1$ and $w_2$ are the weights for the landscape suitability and spatial aggregation degree during the objective function calculation process, respectively; $k$, $l$, and $r$ are the serial numbers for the landscape type, pixel, and patch, respectively; $K$, $Q$, and $R$ are the number of landscape types, grids of landscape raster graphics, and landscape patches, respectively; $d_k(\mathit{A})$ is the grid number for landscape type $k$ in LPSLO scheme $\mathit{A}$; $D_k^{\ast }=z_k^{\ast }/e^2$ ($z_k^{\ast }$ is the optimal area of landscape type $k$ and $e$ is the grid size) is the grid number in LPQSO result for landscape type $k$ in a designed scenario; ${\alpha }_j\in [1,2,\dots ,K](j=1,2,\dots ,N)$ is the encoding range of landscape types corresponding to grids; and $T_j(k\to k^{\prime })$ is the conversion possibility that a landscape type corresponding to grid ${\alpha }_j$ can be transformed from landscape type $k$ to $k^{\prime }$.

(3) Particle initialization: According to the landscape type conversion rules and landscape suitability, the spatial redistribution of the grids corresponding to the LPQSO results is achieved based on the landscape type raster graphics in the base year in order to generate an initial particle with an element number that is equal to the grid number for quantitative structure optimization of the corresponding landscape type. We assumed that $V$ (the structure of which is shown in Equation (17) was the cell array [67] used to store the row and column values of valid elements as well as the suitability values of the corresponding landscape types; $\Delta d_k=\left|d_k(\mathit{A})-D_k^{*}\right|$ is the absolute difference between the grid number $d_k(\mathit{A})$ for landscape type $k$ in the base year and the grid number for quantitative structure optimization for landscape type $k$ in the designed scenario; $sort(Val)$ and $sort(Val, ‘\mathrm{descend}’)$ denote the ascending and descending sorting of variable $Val$, respectively; and $\mathit{A}(x,y)=\alpha$ indicates that the grid of a landscape type, where a row is $x$ and a column is $y$, in LPSLO scheme $\mathit{A}$ is converted into a pixel of the landscape type corresponding to that encoding $\alpha$. Thus, the particle initialization process is illustrated in Figure 5:

$$V=\left\{\left(\begin{array}{ccccc}{x}_1& y_1& p_1(x_1,y_1)& \dots & p_K(x_1,y_1)\\ x_2& y_2& p_1(x_2,y_2)& \dots & p_K(x_2,y_2)\\ ⋮& ⋮& ⋮& & ⋮\\ x_L& y_L& p_1(x_L,y_L)& \dots & p_K(x_L,y_L)\end{array}\right),\dots ,\left(\begin{array}{ccccc}{x}_1& y_1& p_1(x_1,y_1)& \dots & p_K(x_1,y_1)\\ x_2& y_2& p_1(x_2,y_2)& \dots & p_K(x_2,y_2)\\ ⋮& ⋮& ⋮& & ⋮\\ x_W& y_W& p_1(x_W,y_W)& \dots & p_K(x_W,y_W)\end{array}\right)\right\}$$

(17)

where, $V\left\{1\right\},V\left\{2\right\},\dots ,V\left\{K\right\}$ are the row and column values, and the corresponding suitability values of valid elements for farmland, orchard, forest, urban-rural residential and industrial-mining, and waters landscapes, respectively; $(x_L,y_L),(x_G,y_G),(x_S,y_S),(x_C,y_C),(x_W,y_W)$ are the row and column numbers for farmland, orchard, forest, urban-rural residential and industrial-mining, and waters landscapes in matrix ${\mathit{A}}_{m\times n}$, respectively; and $p_K(x_L,y_L),p_K(x_G,y_G),p_K(x_S,y_S),p_K(x_C,y_C),p_K(x_W,y_W)$ are the suitability values for landscape type $K$ corresponding to the elements of farmland, orchard, forest, urban-rural residential and industrial-mining, and waters landscapes, respectively.

(4) Updating the particle velocity and location: We assume that ${\mathit{P}}_{M\times (4N+1)}$ is the storage matrix for the locations, velocities, and current fitness (the objective function values) of $M$ particles, where in matrix $\mathit{P}$, the column $1~2N$ is the particle location variable $(x_{i1},x_{i2},\dots ,x_{iN},y_{i1},y_{i2},\dots ,y_{iN})(i=1,2,\dots ,M)$, the column $2N+1~4N$ is the particle velocity variable $(\delta x_{i1},\delta x_{i2},\dots ,\delta x_{iN},\delta y_{i1},\delta y_{i2},\dots ,\delta y_{iN})$, and the column $4N+1$ is the current fitness $p{F}_i$; ${\mathit{O}}_{(2M+1)\times 2N}$ is the storage matrix for the historical best values, local best values, and global best values of $M$ particles, where in matrix $\mathit{O}$, row $1~M$ is the historical best value variable $(x_{i1},x_{i2},\dots ,x_{iN},y_{i1},y_{i2},\dots ,y_{iN})$ of a particle, row $M+1~2M$ is the local best value variable $(l{x}_{(M+i)1},l{x}_{(M+i)2},\dots ,l{x}_{(M+i)N},l{y}_{(M+i)1},l{y}_{(M+i)2},\dots ,l{y}_{(M+i)N})$ of a particle, and the row $2M+1$ is the global best value variable $(g{x}_1,g{x}_2,\dots ,g{x}_N,g{y}_1,g{y}_2,\dots ,g{y}_N)$ of a particle. The velocity and location are constantly updated for a particle according to its historical best value and local best value during the flight process. The operation equation for the particle velocity is expressed as:

$$\{\begin{array}{l}{\mathit{P}}_t(i,2N+j)=\omega (t){\mathit{P}}_{t-1}(i,2N+j)+c_1\cdot {\mu }_1\cdot [{\mathit{O}}_{t-1}(i,j)-{\mathit{P}}_{t-1}(i,j)]+c_2\cdot {\mu }_2\cdot [{\mathit{O}}_{t-1}(M+i,j)-{\mathit{P}}_{t-1}(i,j)]\\ {\mathit{P}}_t(i,3N+j)=\omega (t){\mathit{P}}_{t-1}(i,3N+j)+c_1\cdot {\mu }_1\cdot [{\mathit{O}}_{t-1}(i,N+j)-{\mathit{P}}_{t-1}(i,N+j)]+c_2\cdot {\mu }_2\cdot [{\mathit{O}}_{t-1}(M+i,N+j)-{\mathit{P}}_{t-1}(i,N+j)]\end{array}$$

(18)

where ${\mathit{P}}_t(i,2N+j)$ and ${\mathit{P}}_t(i,3N+j)$ are the longitudinal and transverse velocities of a particle in iterations $t$, respectively; $t=1,2,\dots I_{\max }$ is the number of iterations, where $I_{\max }$ is the maximum number of iterations; $\omega (t)$ is the inertia weight for iterations $t$; $c_1$ and $c_2$ are the different acceleration weights; ${\mu }_1$ and ${\mu }_2$ are random numbers within the [0, 1] range; ${\mathit{O}}_{t-1}(i,j),{\mathit{O}}_{t-1}(i,N+j)$ is the historical best value for a particle in iterations $t-1$; ${\mathit{O}}_{t-1}(M+i,j),{\mathit{O}}_{t-1}(M+i,N+j)$ is the local best value for a particle in iterations $t-1$; and ${\mathit{P}}_{t-1}(i,j),{\mathit{P}}_{t-1}(i,N+j)$ is the current best position of the particle in iterations $t-1$. In this study, we consider two constraints on a particle’s longitudinal velocity, i.e., when ${\mathit{P}}_t(i,2N+j)>\upsilon x_{\max }$, ${\mathit{P}}_t(i,2N+j)=\upsilon x_{\max }$; and when ${\mathit{P}}_t(i,2N+j)<\upsilon x_{\min }$, ${\mathit{P}}_t(i,2N+j)=\upsilon x_{\min }$. The constraints on a particle’s transverse velocity are similar to those on the longitudinal velocity. The operation equation for a particle’s location can be expressed as:

$$\{\begin{array}{c}{\mathit{P}}_t(i,j)=int[{\mathit{P}}_{t-1}(i,j)+{\mathit{P}}_t(i,2N+j)]\\ {\mathit{P}}_t(i,N+j)=int[{\mathit{P}}_{t-1}(i,N+j)+{\mathit{P}}_t(i,3N+j)]\end{array}$$

(19)

where the constraints on the particle’s longitudinal location can be expressed as: if ${\mathit{P}}_t(i,j)>m$, then ${\mathit{P}}_t(i,j)=m$; and if ${\mathit{P}}_t(i,j)<1$, then ${\mathit{P}}_t(i,j)=1$. Thus, the constraints on a particle’s transverse location are similar to those on its longitudinal location. The constraints given above can ensure that a particle’s movement is restricted to a rectangular range, the length of which is the longest transverse span of the study area and the width is the longest longitudinal span of the study area, but some of the particles will still fly throughout the optimization range because the range is an irregular polygon. Therefore, new constraints are formulated as follows: when the value of the element with a row value of ${\mathit{P}}_t(i,j)$ and column value of ${\mathit{P}}_t(i,N+j)$ in matrix $\mathit{A}$ satisfies $\mathit{A}[{\mathit{P}}_t(i,j),{\mathit{P}}_t(i,N+j)]\le 0$, the particle’s location $[{\mathit{P}}_t(i,j),{\mathit{P}}_t(i,N+j)]$ is changed back to the original location.

LPSLO Model Solution

In the LPSLO model, the solution steps are as follows.

Step 1: The particle locations are initialized as the row and column values of the valid elements in matrix ${\mathit{E}}_{m\times n}^{*}$, and the particle velocities are then initialized as random integers within the region of the particle velocity constraint variable $rg$, where they are generated using the interior function RANDINT in MATLAB.

Step 2: The current fitness of a particle is calculated according to Equation (16) and the local best value of the particle within the region of the first iterations is determined by the annular topology structure method [68], before the global best value of the particle is determined according to the maximum current fitness.

Step 3: The inertia weight of a particle is calculated based on a linearly decreasing formula according to Ma et al. [40].

Step 4: The velocity and location of a particle are updated according to Equations (18) and (19).

Step 5: The current fitness $p{F}_i(t)$ and historical fitness $h{F}_i(t)$ of a particle in iterations $t$ are calculated according to Equation (16), before comparing the fitness of these particles, and the historical best value of a particle is then updated when the particle’s fitness satisfies $p{F}_i(t)>h{F}_i(t)$.

Step 6: The local best value and the corresponding local fitness $l{F}_i(t)$ for particle $i$ within the region of iterations $t$ are extracted by the annular topology structure method [68]. The local fitness $l{F}_i(t-1)$ of a particle in iterations $t-1$ is then calculated according to Equation (16), and when the fitness of this particle satisfies $l{F}_i(t)>l{F}_i(t-1)$, the local best value for the particle is updated.

Step 7: The maximum current fitness $\max p{F}_i(t)$ and the corresponding location of the particle in iterations $t$ are extracted. The global fitness $gF(t-1)$ of a particle in iterations $t-1$ is calculated according to Equation (16) and when the fitness of this particle satisfies $\max p{F}_i(t)>gF(t-1)$, the global best value for the particle is updated.

Step 8: When the iterations $t$ satisfy $t\le I_{\max }$, the algorithm executes $t=t+1$ and returns to steps 3–7; otherwise, it returns to step 9.

Step 9: The global best value for a particle in iterations $I_{\max }$ is exported for mapping according to the spatial mapping relationship shown in Figure 4, where the LPSLO map is finally produced.

2.3. Implementation of the LPOA Model

2.3.1. Implementation of the LSE Model

As described previously [69,70,71] and according to the basic principles of LSEs [72], we selected various indexes in terms of natural factors, neighborhood factors, and socioeconomic factors to construct the index system for LSEs (

Table 2. Indexes used for landscape suitability evaluations (LSEs) in Longquanyi District.

Index Categories	Evaluation Indexes
Natural factors	Elevation, slope, aspect, hypsography degree, average annual rainfall, average annual temperature, soil organic matter content
Neighborhood factors	Nearest distance from city center, nearest distance from town center, nearest distance from major roads, nearest distance from major waters
Socioeconomic factors	Population density, per capita GDP

). In the study area, the landscape spatial distribution is influenced greatly by topographic differences due to the three common terrain types comprising plains, mountainous regions, and hills, so the elevation, slope, aspect, and hypsography degree were used to denote the topography impact factors for the landscape pattern [69,71]. Climate is the major factor that determines the landscape distribution, and the spatial variations in temperature and rainfall in the study area are significant due to differences in the surface morphology. Thus, we selected the average annual rainfall and temperature as climatic factors that affected the landscape distribution. Landscape spatial patterns are affected by the soil distribution to some extent, while the soil organic matter content is the main index that reflects soil fertility as well as forming the basis of the soil classification [73], and thus it was selected to denote the effects of soil factors on the landscape distribution. According to the first law of geography [74], landscape patterns are affected by relevant neighborhood factor, as confirmed in many similar studies [25,47,70,71]. Hence, considering the actual situation in the region, the nearest distances from the city center, town center, major roads, and major waters were employed to denote the effects of neighborhood factors on the landscape pattern. At small spatio-temporal scales, the effects of socioeconomic factors on landscape patterns are more intense than those of natural factors [72]. The study area comprised a suburban area of a metropolis and the population distribution was imbalanced because people migrate continually to cities as a consequence of rapid development in terms of urbanization and industrialization. Furthermore, regional differences in economic development were obvious. Thus, the effects of population and economics on the landscape pattern had great spatial differences. Therefore, the population density and per capita GDP were used to denote the effects of socioeconomic factors on the landscape distribution.

Using the indexes mentioned above and ESRI ArcGIS, we implemented the LSE model of the study area as follows. First, we extracted farmland, orchard, forest, urban-rural residential and industrial-mining, and waters landscape maps from the landscape map of the study area for the year 2014. The locations in the grids where landscape was present were assigned values of 1; otherwise, the grid locations were assigned values of 0. Second, we extracted the grid values from each landscape type map and the corresponding spatial distribution raster graphics for various evaluation indexes, and we then used these extracted grid values to determine each landscape pattern influence index and the corresponding regression coefficients by using stepwise regression for BLR analysis with SPSS software. Finally, we calculated the probability map for each landscape type spatial distribution, i.e., a suitability map for each landscape type spatial distribution, according to Equation (1) by Python programming based on the effect indexes and the corresponding regression coefficients for each landscape type.

2.3.2. Implementation of the LPQSO Model

Implementation of the Objective Functions

In the objective functions of the LPQSO model, the values of the optimization objectives (Equations (2)–(4)) were obtained by calculating the summed economic gross product value for all landscapes types, the summed ecological security index for all landscapes types, and the comprehensive benefits of various landscapes types. According to Equations (2)–(4), calculating the target values required the determination of the product value coefficient for each landscape type, the ecological security coefficient for each landscape type, the standardized product value coefficient for each landscape type, the standardized ecological security coefficient for each landscape type, the weight of the economic development scenario, and the weight of the ecological protection scenario. The product value coefficient for each landscape type was expressed as the total major agricultural and industrial production divided by the corresponding landscape type area. The ecological security index was obtained by calculating the summed grid values of various index raster data using a weighted sum method, which was performed by embedding the comprehensive evaluation index model into GIS based on the ecological security evaluation index system [75]. This index system comprised 25 indicators, including the population density, natural population growth rate, per capita water area, per capita grain output, quantity of chemical fertilizer applied per unit area of cultivated land, urbanization level, industrial added value, regional development index, elevation, slope, grain yield per unit area of cultivated land, per capita construction land area, annual average rainfall, annual average temperature, soil organic matter content, soil type, landscape fragmentation, area-weighted average patch fractal dimension, contagion index, NDVI index, population mortality, agricultural mechanization level, tertiary industry proportion, per capita GDP, and landscape types [75]. The summed ecological security index for each landscapes type was calculated using the SUM function in the Zonal Statistics Tool in ArcGIS, where the input feature zone data in this tool were set as the landscape type vector data and the input value raster in this tool was set as the raster data for the ecological security index in the base year. After obtaining the summed ecological security index for each landscapes type, the ecological security coefficient for each landscape type was expressed as the summed ecological security index for each landscapes type divided by the corresponding landscape type area. The standardized product value coefficient and the standardized ecological security coefficient were calculated for each landscape type using the linear normalization method. Furthermore, considering that economic development is as important as ecological protection, the objective function weights for the economic development scenario and ecological protection scenario were both set to 0.5.

Implementation of the Constraints

The LPQSO model includes many constraints in terms of the landscape area, ecological service value, nonpoint source pollution, and industrial structure. Thus, the constraints were implemented as follows.

First, according to Equations (5)–(10), implementing the landscape area constraints required the determination of the minimum demand areas and maximum areas for farmland, orchard, forest, urban-rural residential and industrial-mining, and waters in each target year. The minimum demand area for farmland was determined based on the minimum per capita farmland area and anticipated registered population in the target year. The minimum per capita farmland area in Longquanyi District during 1988–2014 was 0.0106 hm², and the anticipated registered population in the target year was calculated by a first order linear regression model based on the demographic data for 1978–2014. After obtaining the minimum per capita farmland area and the anticipated registered population, the minimum demand area for farmland in the target year was equal to the product of both. According to the requirement that the basic fruit demand for residents in the study area was met, the minimum demand area for orchards was determined by calculating the annual mean yield per unit area for orchards and the anticipated total yield for orchards in the target year. The annual mean yield per unit area for orchards was determined according to the fruit yield during 2000–2012, and the anticipated total yield for orchards in the target year was calculated based on the mean level of fruit consumption by residents and the resident population in the study area. After obtaining the anticipated total yield for orchards and the annual mean yield per unit area for orchards, the minimum demand area for orchards in the target year was equal to the former divided by the latter. In order to consider the establishment of the national ecological demonstration district in Longquanyi, the minimum demand area for forest in the target year was set as the current area. Given that it would be almost impossible to reduce the urban-rural residential and industrial-mining area because the study area was located in suburban district of Chengdu City, the minimum demand area in the target year was set as the current area. The eastern section of the study area comprises arid mountains and water shortages are very severe, so the minimum demand area for waters in the target year was set as the current area to ensure the safety of the regional water resources. The maximum area for urban-rural residential and industrial-mining was obtained by determining the construction land areas in the city and countryside, the industrial land area, and the traffic land area in the target year. The construction land areas in the city and countryside were calculated based on the maximum residential land area per capita, where both were determined as described previously [76], and the anticipated residential population of both areas in the target year, which were calculated by a first order linear regression model using the resident populations in the city and countryside during 2001–2014. The industrial land area was set as the planning area of the National Chengdu Economic and Technological Development Zone. The traffic land area was set as 15% of the total area of city construction land, countryside construction land, and industrial land [77]. After obtaining the maximum residential land area per capita and the anticipated resident populations of the city and countryside areas, the industrial land area, and traffic land area, we calculated the maximum area of urban-rural residential and industrial-mining in the target year according to Equation (9). In addition to urban-rural residential and industrial-mining, the maximum areas of other landscape types were all set as the total areas of the various landscape types in the study area.

Next, according to Equation (11), implementing the ecological service value constraints required the determination of the ecological service value per unit area for each landscape type, and the minimum and maximum ecological service values for the study area. The ecological service value per unit area was determined according to the ecosystem ecological service value equivalence factor, regional average grain yield per unit area, and the market price using the ecosystem service value calculation methods proposed by Xie et al. [65]. The ecosystem ecological service value equivalence factor was determined according to the table of the “ecosystem service values per unit area for the Chinese terrestrial ecosystem” [65]. This table provides the ecological service value equivalence factors for gas regulation, climate regulation, water conservation, soil formation and protection, waste disposal, biodiversity protection, food production, raw materials, and entertainment culture functions for five types of terrestrial ecosystems comprising forests, orchards, farmland, wetlands, waters, and deserts. The ecological service value equivalent factors for each function of the farmland, forest, and waters landscapes were determined by referring to this table. Thus, due to their absence from this table, the ecological service value equivalence factors for various functions of orchards, urban-rural residential and industrial-mining were estimated indirectly according to the ecological service value equivalence factors for the various functions of each ecosystem provided in this table. The ecological service value equivalence factor for each function of orchards was approximately equal to the mean of the ecological service value equivalence factors for the corresponding functions of forest and farmland in this table. The ecological service value equivalence factor for the entertainment culture function of urban-rural residential and industrial-mining was equal to the mean equivalence factors for the corresponding functions of forest, orchard, farmland, and wetland, and the equivalence factors for the other service functions of urban-rural residential and industrial-mining were assumed to be zero. After determining the ecological service value equivalence factors for various functions of each landscape type, the regional average grain yield per unit area and the market price in the base year, the ecological service value per unit area of landscape type could be expressed as the sum of 1/7th of the product of the ecological service value equivalence factors, grain yield per unit area, and the market price for each landscape function. Subsequently, the minimum and maximum ecological service values for the study area were equal to the sum of the product of the ecological service value per unit area for each landscape type, and the corresponding minimum and maximum constraints as mentioned above.

According to Equations (12)–(14), implementing the nonpoint source pollution constraints requires the calculation of the minimum and maximum annual loads for COD, TN, and TP. Previous studies of nonpoint source pollution [78,79] indicate that calculating the values of the constraints requires the determination of the mean annual rainfall, the runoff coefficient for each landscape type, the concentrations of COD, TN, and TP in the surface runoff from each landscape type, and the amounts of fertilizer applied to the main fruits and crops. The mean annual rainfall was calculated based on rainfall monitoring data from the study area, the amounts of fertilizer applied to the main fruits and crops were obtained from field surveys, and the values of other variables were obtained from previous studies [65,78]. After obtaining these variable values, the annual loads per unit area were estimated for COD, TN, and TP according to Equations (12)–(14), and the minimum and maximum annual loads of COD, TN, and TP were then calculated according to the corresponding minimum and maximum landscape constraints for the area as mentioned above.

Finally, according to Equation (15), implementing the industrial structure constraints required the determination of the minimum and maximum secondary industry product values divided by the summed product values for primary industry and service industry. Based on the macroeconomic situation and industrial development status in Longquanyi District, the growth rates for primary industry, secondary industry, and service industry were set as 2%, 12%, and 8%, respectively. Subsequently, the minimum and maximum secondary industry product values divided by the summed product values for primary industry and service industry were estimated based on industrial development statistical data and the corresponding growth rates given above.

In summary, the values of the main variables for each constraint on the LPQSO model of Longquanyi District in this study are shown in

Table 3. Values of the main variable for each constraint on the landscape pattern quantitative structure optimization (LPQSO) model of Longquanyi District.

Constraint Variable		By 2021	By 2028
Landscape area (hm2)	Total area of each landscape type	55,569	55,569
	Minimum demand area for farmland	7235	7901
	Minimum demand area for orchard	4777	5216
	Minimum demand area for forest	5167	5167
	Minimum demand area for urban-rural residential and industrial-mining	16,207	16,207
	Maximum demand area for urban-rural residential and industrial-mining	19,754	21,586
	Minimum demand area for waters	1610	1610
Ecological service value (million CNY)	Lowest ecological service value in the studied area	669.72	691.80
	Highest ecological service value in the studied area	10,081.07	10,086.88
Nonpoint source pollution (kg)	Minimum annual load of COD	9,412,702	9,481,387
	Maximum annual load of COD	17,895,612	18,864,542
	Minimum annual load of TN	682,316	715,596
	Maximum annual load of TN	3,863,575	3,897,836
	Minimum annual load of TP	138,262	141,448
	Maximum annual load of TP	470,088	481,631
Industrial structure	Minimum secondary industry product value divided by the summed product values for primary industry and service industry	55,569	55,569
	Maximum secondary industry product value divided by the summed product values for primary industry and service industry	7235	7901

.

2.3.3. Implementation of the LPSLO Model

According to the LPSLO model and solution algorithm, five main parameters had to be determined, i.e., the landscape type weight for each scenario, optimal grid number for each landscape type, landscape suitability, landscape type encoding, and landscape type conversion rules. The landscape type weight was calculated for each scenario using the analytic hierarchy process method. We invited 15 experts in related fields (four experts in ecology, four experts in economics, four experts in land, and three experts in agriculture) to construct the judgment matrix based on three rounds of scoring, before calculating the weights for each landscape type in each scenario. The calculated results are shown in

Table 4. Weights for each landscape type in each scenario in Longquanyi District.

Scenario	Weight
	Farmland	Orchard	Forest	Urban-Rural Residential and Industrial-Mining	Waters
Economic development	0.1378	0.2107	0.0606	0.5514	0.0395
Ecological protection	0.0706	0.1366	0.5071	0.0350	0.2507
Overall consideration	0.0930	0.1613	0.3582	0.2072	0.1803

. The optimal grid number for each landscape type was equal to the landscape type area, which was equal to the landscape type optimal area calculated by the LPQSO model minus the corresponding landscape type area in the priority planning area, divided by the grid size setting as 60 m. The landscape suitability was obtained using the LSE model. Next, based on the landscape types in the study area, the landscape type codes for farmland, orchard, forest, urban-rural residential and industrial-mining, and waters in the areas, except for priority planning areas, were set as 1, 2, 3, 4, and 5, respectively, and the landscape type codes in the priority planning areas were set as a negative value. Finally, according to the characteristic changes in the landscape patterns in the study area [80], four landscape type conversion rules were constructed, as follows:

(1)

Urban-rural residential and industrial-mining, and farmland or orchard could be exchanged with each other fully;

(2)

Waters and farmland could be exchanged with each other fully;

(3)

Forest and farmland or orchard could be exchanged with each other fully;

(4)

Farmland and orchard could be exchanged with each other fully.

In addition, according to the actual situation in the study area and based on trial-and-error experiments, the configuration parameters for the LPSLO model and the algorithm were set as follows: maximum iterations $I_{\max }=100$, accelerated weight $c_1=c_2=100$, minimum longitudinal and transverse velocity of a particle $u{x}_{\min }=-6$, maximum longitudinal and transverse velocity of a particle $u{x}_{\max }=6$, maximum longitudinal span $m=478$, maximum transverse span $n=502$, and the inertia weight $\omega (t)$ was determined using a linearly decreasing formula, as described previously [40].

3. Results

3.1. Evaluation Results Obtained from the LSE Model

Using the LSE model, we determined the evaluation indexes and corresponding regression coefficients for each landscape type, before obtaining the spatial suitability maps for each landscape type, as shown in

Table 5. Regression coefficients and significant test results obtained by the binary logistic regression (BLR) models.

Evaluation Indicator	Regression Coefficient a
	Farmland	Orchard	Forest	Urban-Rural Residential and Industrial-Mining	Waters
Elevation	−6.37 × 10−3	2.34 × 10−3	−1.49 × 10−3	——	−8.56 × 10−3
Slope	——	——	——	——	——
Aspect	——	——	——	——	——
Hypsography degree	−1.19 × 10−2	1.17 × 10−3	2.69 × 10−3	−2.86 × 10−3	4.83 × 10−3
Nearest distance from city center	6.99 × 10−5	9.93 × 10−5	1.38 × 10−4	−1.90 × 10−4	——
Nearest distance from town center	2.72 × 10−4	——	——	−1.08 × 10−4	−1.57 × 10−4
Nearest distance from major roads	−9.97 × 10−5	——	−1.63 × 10−4	——	——
Nearest distance from major waters	——	——	−1.59 × 10−4	3.45 × 10−4	−6.65 × 10−4
Mean annual rainfall	−3.19 × 10−3	6.44 × 10−3	7.14 × 10−3	−1.10 × 10−2	——
Mean annual temperature	——	——	−1.21 × 10	3.72 × 10−1	——
Soil organic matter content	2.21 × 10−2	——	——	−3.08 × 10−2	——
Population density	——	−1.19 × 10−4	——	——	——
Per capita GDP	−1.72 × 10−2	4.25 × 10−2	4.42 × 10−2	−6.04 × 10−2	——
Constant	2.49 × 100	−7.29 × 100	1.06 × 101	3.82 × 100	1.33 × 100
ROC	0.8040	0.7125	0.7892	0.7858	0.7000

^a In the BLR analysis process, if the significance level of the regression coefficient between a landscape type and evaluation indicator was less than 0.02, the indicator was retained and the regression coefficient is shown in the table. However, if the significance level was more than 0.02, the indicator was deleted and this regression coefficient is not shown in the table.

and Figure 6.

Among the BLR models for all the landscape types (Table 5), the predictive accuracies of the modeling and test data using the BLR model were 85.1% and 85.3% for farmland landscape, respectively, both accuracies were 65.6% for orchard landscape, 90.1% and 89.5% for forest landscape, 82.7% and 81.7% for urban-rural residential and industrial-mining landscape, and 97.2% and 97.6% for waters landscape. Thus, the predictive accuracy was high for each landscape type and the results predicted by the model were highly reliable. In addition, the ROC values obtained by the BLR model for each landscape type were all higher than 0.7 (Table 5). Therefore, the evaluation indicators selected for the regression equation had high explanatory power for the effects of landscape patterns, and thus these indicators were used for suitability evaluations of landscape patterns.

In addition, Figure 6 shows the spatial characteristics of the suitability evaluation results for each landscape type in the study area. The higher suitability areas for farmland were mainly distributed in the towns of Shiling, Huangtu, Hongan, Luodai, Xihe, and east of Wanxing, surrounding Longquan Lake in the town of Chadian. The higher suitability areas for orchard were mainly distributed in the towns of Shanquan and Tongan, east of Baihe, and west of Chadian and Wanxing. The higher suitability areas for forest were mainly distributed in the towns of Shanquan, Wanxing, and Baihe, east of Tongan, and west of Chadian. The higher suitability areas for urban-rural residential and industrial-mining were mainly distributed in the towns on plains in the study area, such as Shiling, Damian, Xihe, Longquan, and Baihe. The higher suitability areas for waters were mainly distributed in the towns on plains in the study area, such as Shiling, Xihe, Huangtu, Hongan, Luodai, Damian, Tongan, and Baihe, as well as Longquan Lake and its surrounding areas to the south of Chadian in the mountains in the study area. According to these characteristic distributions, there was some overlapping in the higher suitability areas for orchard and forest, where the overlapping distribution areas were distributed centrally to the east of Shanquan, Baihe, and Tongan, and to the west of Chadian.

3.2. Optimal Results Obtained by the LPQSO Model

Using the LPQSO model, we obtained the optimal areas for farmland, orchard, forest, urban-rural residential and industrial-mining, and waters under the economic development, ecological protection, and overall consideration scenarios in the target years of 2021 and 2028, as shown in

Table 6. Quantitative structure optimization results for each landscape type area under each scenario for the target year in Longquanyi District (hm²).

Target Year	Scenario	Farmland	Orchard	Forest	Urban-Rural Residential and Industrial-Mining	Waters
2021	Economic development	7235.00	21,803.00	5167.00	19,754.00	1610.00
	Ecological protection	7235.00	7380.75	23,136.25	16,207.00	1610.00
	Overall consideration	7235.00	16,231.31	10,738.69	19,754.00	1610.00
2028	Economic development	7901.00	19,305.00	5167.00	21,586.00	1610.00
	Ecological protection	7901.00	5216.00	24,635.00	16,207.00	1610.00
	Overall consideration	7901.00	5216.00	19,256.00	21,586.00	1610.00
Base year 2014	Current landscape pattern	6723.00	25862.00	5167.00	16,207.00	1610.00

. Compared with the area for each landscape type in the year 2014, the optimal areas for each landscape type under each scenario in the target year (Table 6) exhibited specific characteristics, as follows.

Under the economic development scenario, the areas for urban-rural residential and industrial-mining and farmland increased, whereas the area for orchard decreased, while the areas for forest and waters were unchanged. Moreover, the area for urban-rural residential and industrial-mining increased the most, whereas the area for orchard decreased the most. Thus, the scheme can achieve the maximum objective in terms of economic benefit because the economic benefit of urban-rural residential and industrial-mining is higher than that of orchard. Hence, the results obtained by the scheme agreed with the actual economic development scenario.

Under the ecological protection scenario, the areas for forest and farmland increased, whereas the area for orchard decreased, while the areas for urban-rural residential and industrial-mining and waters were unchanged. Moreover, the area for forest increased the most, whereas the area for orchard decreased the most. Thus, orchard was changed to forest by the scheme to achieve the maximum ecological security degree because the ecological security associated with forest is higher than that linked with orchard. Hence, the results obtained by the scheme agreed with the actual ecological protection scenario.

Under the overall consideration scenario, the areas for forest, urban-rural residential and industrial-mining, and farmland increased, whereas the area for orchard decreased, while the area for waters was unchanged. Moreover, the area for forest increased the most, followed by the area for urban-rural residential and industrial-mining, and the area for orchard decreased the most. Thus, orchard was changed to forest and urban-rural residential and industrial-mining, thereby significantly improving the maximum comprehensive benefit, including economic and ecological benefits, because the ecological security degree associated with forest and the economic benefit of urban-rural residential and industrial-mining are higher than those linked with orchard. Hence, the results obtained by the scheme agreed with the actual overall consideration scenario.

3.3. Optimal Results Obtained by the LPSLO Model

According to the basic data, such as the landscape map for the base year 2014 (Figure 2a), landscape maps of the priority planning areas (Figure 2b), spatial suitability evaluation maps for each landscape type (Figure 6), and quantitative structure optimization areas for each landscape under each scenario in the target year (Table 6), we obtained the landscape pattern spatial layout schemes for each scenario in the target years of 2021 and 2028 by using the LPSLO model and solution algorithm implemented in MATLAB, as shown in

Table 7. Solutions obtained by the LPSLO model under each scenario for the target year in Longquanyi District.

Target Year	Designed Scenario	Farmland	Orchard	Forest	Urban-Rural Residential and Industrial-Mining	Waters
2021	Economic development	19,870 (0.98%)	60,677 (0.34%)	14,600 (1.88%)	54,541 (0.45%)	4435 (0.67%)
	Ecological protection	19,796 (1.35%)	20,950 (2.34%)	64,261 (0.14%)	44,664 (0.64%)	4452 (0.29%)
	Overall consideration	19,890 (0.88%)	45,121 (0.23%)	30,165 (1.28%)	54,539 (0.46%)	4408 (1.28%)
2028	Economic development	21,754 (0.73%)	53,271 (0.51%)	14,829 (3.47%)	59,825 (0.08%)	4444 (0.47%)
	Ecological protection	21,810 (0.47%)	14,734 (1.85%)	68,246 (0.12%)	44,905 (0.10%)	4428 (0.83%)
	Overall consideration	21,840 (0.34%)	14,533 (0.46%)	53,570 (0.31%)	59,705 (0.28%)	4475 (0.22%)
Base year 2014	Current landscape pattern	22,136	82,917	14,268	30,439	4363

The statistical analysis of the area of each landscape type was conducted according to the number of grid units measuring 60 m. The values in brackets show the relative error between the optimal areas obtained by LPSLO and LPQSO.

and Figure 7. Figure 7 shows that the optimization results obtained for the landscape pattern spatial layout under each scenario in the target year had specific characteristics, as follows.

Under the economic development scenario, in the target year 2021, farmland was distributed throughout towns on the plains, including Huangtu, Hongan, and Xihe, as well as east of towns in mountainous areas such as Wanxing and Chadian; orchard was distributed in towns in mountainous areas, such as Shanquan, Chadian, and Wanxing, as well as east of Luodai, Tongan, and Baihe; forest was distributed in the towns of Wanxing and Chadian, as well as east of the towns of Luodai, Tongan, and Baihe; urban-rural residential and industrial-mining was distributed centrally in the towns of Shiling, Xihe, Damian, Longquan, Huangtu, and Hongan, and west of Luodai, Tongan, and Baihe; waters were mainly distributed in lakes in the mountainous region to the east of the study area, and in ponds and canals in the plains to the west of the study area. Compared with the landscape spatial layout under the economic development scenario in the target year of 2021, the areas covered by urban-rural residential and industrial-mining and farmland increased in the target year of 2028, whereas the area covered by orchard decreased.

Under the ecological protection scenario, in the target year of 2021, farmland was distributed throughout towns on the plains, including Huangtu, Hongan, Xihe, and Longquan, as well as the areas surrounding Longquan Lake in the town of Chadian in the mountainous region of the study area; orchard was distributed in towns on the plains such as Huangtu, Hongan, and Luodai, as well as southeast of Chadian in the mountainous region of the study area; forest covered the whole areas of the towns of Shanquan and Wanxing, as well as east of Baihe, Tongan, and Luodai, and northwest of Chadian; urban-rural residential and industrial-mining was distributed centrally in the towns of Shiling, Xihe, Damian, and Longquan, as well as west of the towns of Tongan and Baihe; while the characteristic distribution of the waters landscape spatial layout was similar to that under the economic development scenario. Compared with the landscape spatial layout under the ecological protection scenario in the target year of 2021, the areas covered by forest and farmland increased in the target year of 2028, whereas the area covered by orchard decreased.

Under the overall consideration scenario, in the target year 2021, farmland was distributed throughout the towns on the plain, including Huangtu, Hongan, and Xihe, as well as east of Wanxing and Chadian in the mountainous region of the study area; orchard was distributed mainly in the towns of Luodai and Chadian, as well as east of Tongan and Baihe, but less widely in the towns of Shanquan, Wanxing, Huangtu, Xihe, and Longquan; forest was distributed mainly to the east of Tongan and Baihe, and in the towns on the plains, such as Wanxing and Shanquan; urban-rural residential and industrial-mining was distributed mainly in the towns of Shiling, Xihe, Damian, and Longquan, and west of Tongan, Baihe, and Luodai, but less widely in the towns of Huangtu and Hongan; while the distribution characteristics of the landscape spatial layout for waters was similar to that under the economic development scenario as well as the ecological protection scenario. Compared with the landscape spatial layout under the overall consideration scenario in the target year of 2021, the area covered by forest increased in the target year of 2028 to cover the whole of the towns of Wanxing and Chadian, as well as east of the towns of Luodai, Tongan, and Baihe, and northwest of Chadian, whereas the area covered by orchard decreased greatly in the target year of 2028, where it was mainly distributed southeast of Chadian and on the piedmont to the east of Luodai.

4. Discussion

4.1. Analysis of the LPOA Model for Longquanyi

4.1.1. Accuracy of the Spatial Optimization Results

There were differences in the optimal areas obtained by LPSLO and LPQSO for each scenario in the target year (Table 5 and Table 6). The maximum error obtained after spatial optimization was for the forest area under the economic development scenario in the target year of 2028, where the relative error before and after optimization with the LPSLO model was 3.47%. The minimum error obtained after spatial optimization was for the urban-rural residential and industrial-mining area under the economic development scenario in the target year of 2028, where the relative error before and after optimization with the LPSLO model was 0.08%. Thus, the solutions obtained by the LPSLO model could not satisfy the equality constraints because of the effects of some particle elements during the flight process. Therefore, the solutions obtained by the model could not satisfy the equality constraints because the locations of particle elements corresponding to certain landscape types were replaced by particle elements corresponding to other landscape types during the flight process (i.e., the first particle element moved to a location while the second particle element moved to the same location occupied previously), thereby decreasing the previous landscape type area but increasing the subsequent landscape type area. According to the principles of the LPSLO model and the algorithm, the initial particles fully satisfy the equality constraints at initialization. However, the equality constraints are broken by some particles because the objective function values of some particles that break the equality constraints are higher than those of the initial particles. Thus, the landscape pattern spatial layout undergoes further optimization. Moreover, breaking the equality constraints in a specific error range may be a strategy for resolving the coupling problem between the objective functions for LPQSO and LPSLO. Thus, the LPSLO model can be used to optimize the landscape pattern spatial layout.

Furthermore, in order to validate the reliability and feasibility of the LPSLO model proposed in this study, we quantitatively assessed the accuracy of the landscape pattern simulation results using the widely employed confusion assessment method [46], where we compared the simulated and actual landscape pattern maps, and the accuracy was computed based on point-by-point comparisons using a confusion matrix. First, considering that Longquanyi District has been undergoing a period of rapid economic development and urbanization in the past ten years, we set the current landscape pattern status in the base year of 2014 as the actual landscape pattern for the economic development scenario. Next, according to the grid number for each landscape type in the current landscape pattern status in 2014, we obtained the simulated spatial allocation of the landscape pattern under the economic development scenario by using the LPSLO model with the parameters given above. Using ENVI 5.1, the confusion matrix shown in

Table 8. Confusion matrix for the simulated landscape pattern obtained using the LPSLO model and the actual landscape pattern in the base year of 2014.

Ground Truth Classification	Farmland	Orchard	Forest	Urban-Rural Residential and Industrial-Mining	Waters
Farmland	15,280	5218	71	405	83
Orchard	3770	75,499	3354	766	167
Forest	198	6505	7277	39	58
Urban-rural residential and industrial-mining	468	1046	51	29,291	30
Waters	269	575	30	43	3574

was computed by comparing the simulated results obtained by the LPSLO model with the current landscape pattern status map in 2014 for each cell. The results showed that the overall accuracy of the simulated results was 84.98% and the Kappa coefficient was 0.7587, thereby indicating that the LPSLO model could simulate the landscape pattern spatial layout objectively and accurately. Therefore, it is feasible to employ this model to simulate the spatial allocation of future landscape patterns under each scenario based on the LPQSO model and LSE model, and thus the simulated results obtained were acceptable.

4.1.2. Rationality of the Spatial Optimization Results

In order to analyze the rationality of the spatial optimization results, we compared the characteristics of the current landscape pattern spatial status in the base year of 2014 with the landscape pattern optimization spatial layouts for the target years of 2021 and 2028, where the changes were analyzed according to the algorithm for the LPSLO model. As shown in Figure 2 and Figure 7, Table 6 and Table 7, the landscape pattern spatial optimization layout for each scenario in the target year exhibited specific characteristics, as follows.

Under the economic development scenario, the landscape pattern spatial layout obtained by the optimal schemes for the two periods indicated large increases in the urban-rural residential and industrial-mining areas but large decreases in orchard, while the forest and waters areas increased slightly and farmland decreased slightly. According to the algorithm in the LPSLO model, parts of some landscape types with lower suitability would be adjusted for other landscape types with higher suitability, while the optimized quantitative structure area for the former landscape type was decreased in the target years compared with the area of the corresponding landscape types in the base year, whereas the area for the latter landscape type in the target year was better than the area for the corresponding landscape types in the base year. Thus, the urban-rural residential and industrial-mining areas increased greatly compared with the corresponding area in the base year because the optimized quantitative structure areas in the target years were greater than the corresponding area in the base year. In addition, the orchard area decreased greatly because the optimized quantitative structure area was smaller in the target years than the corresponding area in the base year. However, the optimized quantitative structure areas for farmland in the target years were greater than those in the base year, and the optimized quantitative structure areas for forest and waters were the same as those in the base year, while the areas for farmland decreased slightly and the areas for forest and waters increased slightly. This was probably because the objective function value was improved by reducing the area of farmland and by increasing the areas for forest and waters according to the algorithm in the LPSLO model, as well as due to the error caused by particles colliding during the flight process.

Under the ecological protection scenario, the optimized landscape pattern spatial layouts obtained in the two periods indicated increases in the forest and urban-rural residential and industrial-mining areas, but decreases in the farmland and orchard areas, while the waters areas were unchanged. The optimized spatial areas for forest were increased compared with the corresponding area in the base year. This was mainly because the optimized quantitative structure areas in the target year were larger than the corresponding area in the base year, and thus other landscape types that were highly suitable for conversion into forest landscape were converted by the algorithm in the LPSLO model. Thus, the optimized spatial areas for orchard decreased because their optimized quantitative structure areas in the target year were smaller than the corresponding area in the base year, and thus the low suitability regions were converted into other landscape types by the algorithm in the LPSLO model. However, similar to the farmland, forest, and waters types under the economic development scenario, although the optimized quantitative structure areas for farmland were larger in the target years than the base year and the areas for urban-rural residential and industrial-mining were also the same in the target years as the base year, the optimized spatial areas for farmland decreased and the areas for urban-rural residential and industrial-mining increased.

Under the overall consideration scenario, the optimized landscape pattern spatial layouts in the two periods indicated that the areas for urban-rural residential and industrial-mining and forest increased greatly, whereas the areas for orchard and farmland decreased, while the waters areas were unchanged, in the same manner as the ecological protection scenario. Given that the human demand for landscape resources is unlimited, game selection among the various landscape types was necessary during the spatial optimization process because of the finite amount of land to sustain landscape resources. Therefore, the weights for the urban-rural residential and industrial-mining and forest types shown in Table 4 were larger than those for other landscape types, and thus they were preferentially increased in the specified area. Furthermore, waters had a low weight but there was little spatial change in the areas for water areas because a protection strategy for water resources was considered in the LPSLO model (Figure 2). In addition, the difference between the weights was slight for orchard and farmland, but the orchard type decreased greatly, whereas the farmland type only decreased slightly. This was mainly because the high suitability areas for orchard overlapped with those for forest in the study area (Figure 6), and thus orchard was favored less by the game selection where forest had a higher weight than orchard, so the spatial expansion of orchard decreased greatly.

4.1.3. Comparative Analysis of the LPOA Scheme Under Each Scenario

In order to conduct comparative analyses of the LPOA schemes under each scenario, we analyzed the feasibility of the landscape pattern spatial allocation scheme for each scenario in the target years by considering the actual status of socioeconomic development and ecological construction. The Chinese economy has slowed down, but the overall economy and population of Longquanyi District will continue to increase in the future because this area is located in a special location for economic growth, i.e., a suburban area of Chengdu City, and it is the location of the high-end manufacturing industry in Sichuan province. Population growth will inevitably lead to increases in the demand for land for residential and public service uses as well as for economic development, which will increase the demand for industrial land and ancillary land. Hence, the potential growth of urban-rural residential and industrial-mining areas will be very high in this area in the future. Substantial increases in urban-rural residential and industrial-mining areas will decrease the availability of other landscape resources, including farmland, orchard, forest, and waters landscape types because of the finite availability of land resources. However, the growth of forest and waters landscape types has greater ecological benefits and Longquanyi is a National Ecological Demonstration Zone, so a decrease in ecological construction efforts is unlikely, and thus potential decreases in the farmland and orchard landscape types are very likely in the future. Under the economic development scenario for landscape patterns, we found that the urban-rural residential and industrial-mining area increased greatly whereas the orchard area decreased considerably in order to achieve the maximum economic benefit. However, the areas of the forest, waters, and farmland landscape types changed little, especially forest, which has higher ecological value and it increased slightly. This is inconsistent with the actual situation because the growth of forest in this area has high potential. Thus, the feasibility of the schemes obtained under this scenario may be lower in the future. Under the ecological protection scenario for landscape patterns, the area of forest increased greatly when emphasizing the maximum ecological security level. However, in this scheme, there was not an adequate supply of urban-rural residential and industrial-mining land, which maintains steady economic growth, and this is inconsistent with the actual situation where a substantial expansion of land allocated to urban-rural residential and industrial-mining is very likely. Thus, the feasibility of the scheme obtained under this scenario is also less likely in the future. Under the overall consideration scenario for landscape patterns, the scheme considered the supply of land resources required to maintain economic growth as well as the ecological resources required to improve the quality of the environment, so there were adequate allocations to the urban-rural residential and industrial-mining landscape type to maintain economic growth and to the forest landscape type to preserve ecological security. This scheme agrees with the likely changes in the regional landscape as well as providing appropriate allocations to meet the actual needs of economic development and ecological construction in the study area. Thus, the scheme obtained under this scenario is the most likely compared with the other two schemes and it is the optimal scheme for determining the future landscape pattern spatial layout in the study area.

4.2. Discussion of the Potential Usefulness of the LPOA Model

4.2.1. Useful Results obtained from the Application of the LPOA Model

Three useful results were obtained from the application of the LPOA model, i.e., a planning concept for balancing economic and ecological benefits, a planning strategy that considers both quantitative structure optimization and spatial layout optimization, and a planning strategy that closely connects the theoretical model with the actual region, and these results provide a useful reference for formulating related spatial planning activities, including landscape pattern planning, land use planning, and urban planning in other regions.

First, the conflict between economic development and ecological protection has become the principal contradiction that affects sustainable development in human society. During the modernization process throughout the world, many disasters have been caused by ignoring ecological protection while focusing on economic construction, such as nuclear pollution, air pollution, ocean pollution, rainforest felling, and industrial mining. However, the desire for a better life is the tireless pursuit of humanity. We cannot stop economic development to protect the ecological environment or seek economic development at the expense of the environment. Thus, balancing the conflict between the needs of economic development and ecological protection is a major problem which must be solved urgently to allow the sustainable development of human society. The LPOA model optimization results confirmed that LPOA, by balancing economic and ecological benefits, is an effective approach for mitigating conflicts between economic development and ecological protection from the perspective of spatial planning. This planning concept can provide a useful reference for other regions when conducting landscape pattern planning, land use planning, urban planning, and other related spatial planning processes.

Second, in the past, when we formulated spatial planning processes such as land use planning and urban planning, we attached great importance to planning and forecasting the resource quantity in the target year, whereas we often ignored the rational spatial allocation of resources, thereby making some plans difficult to implement because spatial optimization allocation planning was ignored. The LPOA model achieves a successful planning strategy by effectively coupling the landscape quantity structure with the spatial layout during the LPOA process. Other regions can employ this planning strategy to guide the preparation of landscape pattern planning, land use planning, urban planning, and other related forms of spatial planning.

Third, many spatial optimization models were mentioned above, but most are difficult to apply directly for practical planning. This is mainly because some optimization models fail to consider macrofactors, including social, economic, policy, system, and other factors, and thus the optimization results are far from reality and difficult to implement. During the LPOA process, the LPOA model first formulates large urban built-up areas and waters as priority planning areas, before applying the optimization model to spatially optimize the allocation of the landscape pattern in other areas, thereby allowing the effective combination of macrofactors and theoretical models. This approach is universal because when formulating similar types of spatial planning in other regions, specified unadjusted areas (e.g., ecological protection and basic farmland protection zones) can be planned as priority planning zones and the other areas can be subjected to spatial layout optimization allocation using the optimization model, thereby ensuring a close connection between the theoretical model and the reality in the region.

4.2.2. Beneficial Contributions of the LPOA Model

The LPOA model addresses three problems related to the current spatial optimization decision-making methods, such as neglecting the coupling between quantitative structure optimization and spatial layout optimization, ignoring the macrofactors that affect landscape patterns when optimizing modeling, and initializing particles without considering the suitability of the landscape. Our proposed LPOA model is beneficial because it can complement the simulation of spatial optimization allocation processes, such as LPOA, land use spatial allocation, urban space optimization allocation.

First, the key to establishing the LPOA model is integrating the landscape optimal area, landscape suitability, and macrofactors (e.g., social, economic, policy, and system factors) in order to conduct the LPOA to maximize the landscape benefits. We established a new LPOA model for optimizing the allocation of the landscape spatial pattern by integrating the landscape optimal area, landscape suitability, and macrofactors that influence landscape patterns in order to maximize the economic, ecological, and comprehensive benefits of landscape patterns. The optimal landscape areas in various landscape types are optimized with the LPQSO model based on NP. The suitability of each landscape type is then calculated with the LSE model based on BLR. The macrofactors that affect landscape patterns are combined with the LPOA models by setting model constraint conditions and priority planning areas. A PSO evolutionary algorithm is then embedded into the LPSLO model in order to integrate the results obtained by the LPQSO and LSEs, thereby establishing an LPOA model and algorithm to address the problems caused by ignoring the coupling between quantitative structure optimization and spatial layout optimization, the macrofactors that affect landscape patterns when optimizing modeling, and the suitability of landscape types, as found in previous studies.

Second, based on raster data, the focus when optimizing the landscape pattern spatial layout using PSO is how to use the positions of particle elements to simulate the spatial positions of cells for landscape type raster graphics, and establishing the spatial mapping of the relationships between landscape type raster graphics and particles. We assume that the matrix denotes the landscape type raster map and the matrix element value denotes the pixel attribute value of the landscape type raster map. The matrix is abstracted as a particle that corresponds to a landscape pattern spatial layout scheme. The matrix element value is abstracted as particle element that corresponds to a landscape type code. According to the principle of the PSO algorithm, we can program using the powerful matrix operation capability in MATLAB software to simulate the particle flight process to search for the best scheme for the landscape pattern spatial layout, i.e., a matrix element value moves from one location to another by adjusting its row and column numbers to form a new matrix (a new landscape pattern spatial layout scheme). Thus, the LPSLO is obtained successfully when the new matrix corresponding to the landscape pattern spatial layout achieves the maximum optimization objective. The type and number of the particle element (landscape type code) comprising the particle (landscape pattern spatial layout scheme) will never change in terms of its spatial position. Therefore, according to the landscape suitability and landscape type conversion rules, we generate an initial particle with an element number equal to the grid number for the optimal area of the corresponding landscape type based on the landscape type raster map in the base year, which effectively combines landscape pattern quantity structure optimization and spatial layout optimization, as well as skillfully integrating the landscape suitability evaluation results during the optimization of the landscape spatial layout. Therefore, this approach successfully overcomes the problem of neglecting the coupling between quantitative structure optimization and spatial layout optimization for the landscape pattern, as well as addressing the problem of ignoring the landscape suitability when initializing particles.

Third, during the establishment of the spatial mapping relationship between the landscape type raster graphics and particles, we first set the code value of landscape type in the priority planning area as negative values (e.g., setting the code value for urban-rural residential and industrial-mining as −4 or setting the code value for waters as −5), and encode the landscape type code values in other areas as positive values (e.g., setting the code values for farmland, orchard, forest, urban-rural residential and industrial-mining, and waters as 1, 2, 3, 4, and 5, respectively). We then define the elements with codes that have positive values as the valid elements and store them as vector quantity data structures denoting a particle, whereas we define the elements with codes that have negative and null values as invalid elements that do not map to particles, i.e., the invalid elements will not participate in the optimization allocation for the spatial locations. Thus, after each iterative computation of the particles, we amalgamate the landscape pattern spatial layout corresponding to the invalid elements in the priority planning area with the landscape pattern spatial layout searched for by the flight of the particles, and then calculate the particle fitness based on the merged spatial layout to obtain the optimal scheme of LPOA for the entire area. Therefore, this approach allows the effective combination of the spatial optimization model and macroscopic factors, such as society, economy, policy, and system factors, as well as successfully addressing the problem of coupling macroscopic factors that affect the landscape pattern with the spatial optimization model.

5. Conclusions

In this study, we proposed a new composite model called the LPOA model for optimizing the spatial allocation of a landscape pattern, which aims to solve a set of optimal problems in LPOA, such as neglecting the coupling between quantitative structure optimization and spatial layout optimization, ignoring the macrofactors that affects landscape patterns during optimization modeling, and initializing particles without considering the suitability of the landscape. The LPOA model mainly comprises LSEs, LPQSO, and LPSLO, where it successfully integrates BLR and NP with PSO, thereby overcoming the problem of optimizing either the quantitative structure or the spatial layout of the landscape pattern, as found in previous studies, as well as addressing the problem of ignoring the landscape suitability and macrofactors that influence the landscape pattern during the LPOA process. The model proposed in this study is a beneficial and useful complement to methods for simulating the spatial optimization allocation, such as LPOA, land use spatial allocation, and urban space optimization allocation, thereby providing a useful reference for formulating related spatial planning processes, including landscape pattern planning, land use planning, and urban planning, in other regions.

We employed the LPOA model to optimize the landscape pattern for the target years of 2021 and 2028 in Longquanyi District. We found that the LPOA model could simultaneously optimize the quantitative structure and spatial landscape pattern, as well as effectively integrating the landscape suitability and relevant factors that influence landscape patterns, such as social, economic, policy, and system factors. Moreover, the LPOA model could optimize the landscape pattern in terms of its quantitative structure, spatial layout, and benefits, where we established optimized landscape pattern schemes under economic development, ecological protection, and overall consideration scenarios. The model significantly improved the overall economical, ecological, and comprehensive benefits of the landscape pattern. In addition, we assessed and analyzed the accuracy and rationality of the spatial optimization results, where we found that the overall accuracy of the spatial optimal results was 84.98% with a Kappa coefficient of 0.7587. This indicates that performance of the LPSLO model was good and the application of this model can satisfy the demands for LPOA under multiple constraints. Furthermore, the results obtained by the simulated scheme were consistent with the actual situation. The proposed model can provide support and a scientific basis for regional landscape pattern planning, land use planning, urban planning, and other related spatial planning.

LPOA is a complex and multiobjective decision-making process. In this study, we present a new LPOA model that integrates LSE and LPQSO models with an LPSLO model based on the grid units in a landscape type raster map. Our results demonstrated that this model can achieve LPOA by simultaneously combining the quantitative structure optimized using the LPQSO model, the spatial layout optimized using the LPSLO model, the landscape suitability obtained using the LSE model, as well as the macrofactors that affect landscape patterns, including social, economic, policy, and system factors. However, the number of calculations required by this model and the runtime of the algorithm increase with the resolution of the raster images or the size of the study area, so it is necessary to select a raster map with an appropriate resolution according to the study area size and research scale before conducting spatial optimal decision making with this model. Furthermore, by considering the coupling of the objective functions for both the LPQSO and LPSLO models, as well as the efficiency of the algorithm, the proposed model does not include particle collision constraints. The relative error in the calculation results is acceptable for macroplanning in terms of regional landscape security patterns. Nevertheless, it is necessary to specify particle collision constraints after resolving the coupling problem between the objective functions for LPQSO and LPSLO, as well as improving the efficiency of the algorithm for the LPOA model when conducting detailed planning and design for landscape patterns and other related spatial planning.