Land Use Zoning at the County Level Based on a Multi-Objective Particle Swarm Optimization Algorithm: A Case Study from Yicheng, China

Comprehensive land-use planning (CLUP) at the county level in China must include land-use zoning. This is specifically stipulated by the China Land Management Law and aims to achieve strict control on the usages of land. The land-use zoning problem is treated as a multi-objective optimization problem (MOOP) in this article, which is different from the traditional treatment. A particle swarm optimization (PSO) based model is applied to the problem and is developed to maximize the attribute differences between land-use zones, the spatial compactness, the degree of spatial harmony and the ecological benefits of the land-use zones. This is subject to some constraints such as: the quantity limitations for varying land-use zones, regulations assigning land units to a certain land-use zone, and the stipulation of a minimum parcel area in a land-use zoning map. In addition, a crossover and mutation operator from a genetic algorithm is adopted to avoid the prematurity of PSO. The results obtained for Yicheng, a county in central China, using different objective weighting schemes, are compared and suggest that: (1) the fundamental demand for attribute difference between land-use zones leads to a mass of fragmentary land-use zones; (2) the spatial pattern of land-use zones is remarkably optimized when a weight is given to the sub-objectives of spatial compactness and the degree of spatial harmony, simultaneously, with a reduction of attribute difference between land-use zones; (3) when a weight is given to the sub-objective of ecological benefits of the land-use zones, the ecological benefits get a slight increase also at the expense of a reduction in attribute difference between land-use zones; (4) the pursuit of spatial harmony or spatial compactness may have a negative effect on each other; (5) an increase in the ecological benefits may improve the spatial compactness and spatial harmony of the land-use zones; (6) adjusting the weights assigned to each sub-objective can generate a corresponding optimal solution, with a different quantity structure and spatial pattern to satisfy the preference of the different decision makers; (7) the model proposed in this paper is capable of handling the land-use zoning problem, and the crossover and mutation operator can improve the performance of the model, but, nevertheless, leads to increased time consumption.


Introduction
Land-use zoning and regulation, which originated in the late 19th century and was ubiquitous in most major US cities from the 1920s, is an effective measure for optimizing the allocation of land resources [1,2] and mitigating the negative externalities caused by mixed land use, and has been implemented in many countries and regions [3][4][5]. As the most important method of land-use regulation undertaken by local governments, zoning divides an administrative division into geographically contiguous 'zones', where each zone has its ordinance which prescribes what may and what may not be done [6]. The structure and goals of land-use regulation vary with the economic development, social institution and traditional culture of the different countries. China, since the adoption of economic reform and the open-door policy in 1978, has experienced rapid urbanization and industrialization, with the resulting undesirable phenomena of urban sprawl and expansion, loss of farmland, damage to the ecological environment, etc. Therefore, the main purpose of the zoning policy of China is to restrict the urban encroachment on rural lands, slow the agricultural conversion process and provide protection for the preservation of farmland and forest land [7]. Xinji County of Hebei Province, as a pilot county, adopted the nation's first comprehensive zoning ordinance in 1987, but there was no specific legal provision for zoning until 1998 when Article 20 of the China Land Management Law introduced a regulation bringing in CLUP at the county level, which defines the land-use zones and the dominant land use of each zone. It can be seen, however, that the actual implementation of zoning policy occurred later in China, and the methods and laws for zoning and land-use regulation, in line with China's actual land-use conditions, are still in the process of being explored. So, the focus of this paper is to construct a rational model serving land-use zoning in China, in view of the fact that zoning is the centerpiece of local planning and land-use regulations [8].
Standards for CLUP at the county level note clearly that, in general, a county can be divided into eight different zones comprising basic farmland preservation areas (BFPA), general agricultural land areas (GALA), forestry land areas (FLA), animal agricultural land areas (AGLA), construction land for towns and villages (CLTV), independent industrial and mining land areas (IIMLA), tourism land areas (TLA), and natural and humanistic landscape preservation areas (NHLPA). Each land-use patch should be assigned to a certain zone, based on the actual land-use vector map [9,10]. A county can improve the land-use multi-efficiency and achieve the optimal allocation and sustainable utilization of land resources by strictly controlling the usages of land in each zone [11]. Traditional methods such as clustering analysis [12], the aggregative indicator method [13], and principal component analysis [14] are used for land-use zoning. All these methods are conducted based on the natural and socio-economic attributes of the land-use patches. However, land-use zoning is not only a spatial cluster problem, but also a multi-objective optimization problem, with the best contiguity between units belonging to the same zone and the highest economic returns, social effects and ecological benefits of land use [15]. Therefore, the aforementioned methods have obvious drawbacks in coping with land-use zoning with a MOOP. Furthermore, up to now, there's few research on solving land-use zoning problems with multi-objective optimization techniques.
Nevertheless, many multi-objective optimization models have been presented for solving the optimization of land-use structure and land-use allocation, which are also the significant components of land-use planning as land-use zoning in China. For example, Brookes [16] used a genetic algorithm (GA) to find the best configuration for patches subject to multiple criteria; Aerts [17] adopted linear integer programming (LIP) for multi-site land-use allocation, to minimize development costs and maximize compactness; Liu et al. [18] analyzed the application of a GA in the optimization of land-use structure, considering the economic benefits, ecological 'green' equivalent, and soil erosion; Santé et al. [19] applied simulated annealing (SA) in rural land-use allocation, to reconcile multiple conflicting interests such as land suitability and compactness; Sadeghi et al. [20] optimized land resource allocation by using linear programming, to minimize soil erosion and maximize economic returns; Gao et al. [21] generated an optimized land-use structure which is able to reduce soil erosion, enhance the utilization of water resources and raise agricultural productivity, based on a multi-objective programming (MOP) model and the geographic information system (GIS) technique; Eldrandaly [22] indicated that integrating GIS and gene expression programming (GEP) is a promising and efficient approach for solving multi-site land-use allocation (MLUA) problems; and Zhang et al. [23] applied a multi-agent system (MAS) to regional land-use optimization allocation under a multi-objective constraint. Multiple objectives defined in a spatial context leads to a complex spatial optimization problem and often exhibits substantial computational complexity, which makes exact optimization methods such as LIP infeasible when there are a large number of land units [17,24]; however, optimization techniques that incorporate heuristics have been found to be more effective in solving spatial optimization problems such as SA, GA [25]. Therefore, why not build a multi-objective optimization model for land-use zoning according to previous studies? The particle swarm optimization (PSO) method put forward for the first time by Kennedy and Eberhart [26], is a relatively recent heuristic search method whose mechanics are inspired by the swarming or collaborative behavior of biological populations, which is known as swarm intelligence [26]. As an evolutionary algorithm, it is based on a population of candidate solutions (swarm of particles), which have improved capability for solving complex problems, high convergence speed and good generality for different global optimization problems [27]. The most important advantages of PSO, compared to other optimization strategies, are that PSO is easy to implement and there are few parameters to adjust [28]. It therefore has wide applications in solving the MOPs in many fields [29][30][31][32]. Yicheng, a county-level city in Hubei Province, with rapid economic growth and a sensitive ecological environment, was chosen as the study area in this study. PSO was applied to the land-use zoning of Yicheng, employing an objective function that took into account the differences between zones, the spatial compactness of zones, the degree of spatial harmony of land-use zones, as well as the ecological benefits of land-use. The model obtained the Pareto-optimal solutions, with trade-off between the conflicting objectives, and was subject to the three typical constraints of containing a minimum parcel area in the land-use zoning map, limiting the area of land-use zones, and the general rules of dividing land units into various land-use zones. Additionally, considering the shortcomings of premature convergence and possibility of sinking into local optima of a basic PSO algorithm, a constriction factor and crossover and mutation operator were used in this model [33,34]. The algorithm, with several different sets of weights applied to these four objectives, was run, and the corresponding results were compared.
The remainder of the paper is organized as follows: Section 2 gives an introduction to the study area and relevant data. The application of PSO in a land-use zoning problem is described in Section 3. Thereafter, in Section 4, the solutions with different sets of weights are compared, and the improvements are discussed. Finally, the main conclusions are given in Section 5.

Study Area
The study area (31°26′-31°54′N, 111°57′-112°57′E) is located in Hubei Province, central China, with an area of 2,114.75 km 2 ( Figure 1). The area, whose average elevation is between 50 and 150 m, accounts for 76.8% of the total county area, and a dense river network and many ponds cover the whole area. The average annual rainfall in this area, which has a typical subtropical monsoon climate, is between 800 and 1,000 mm. For this area, the population was 5.15 × 10 5 and the average annual growth rate of total gross domestic product (GDP) in the 11th five-year period was more than 10%. The study area has abundant mineral resources, including silica, limestone, and the biggest deposits of bauxite in Hubei Province. Extensive industrial development of the area has led to a great threat to the local ecological environment.
In accordance with the CLUP (2006-2020) for Yicheng, the land-use zone of CLTV should continue to be divided into two zones comprising the construction land for towns (CLT) and construction land for villages (CLV). In addition, because there is little grassland for animal agriculture in this study area, the zone of LLA is excluded. In the end, the scenario of land-use zoning for Yicheng makes sure that eight land-use zones, including BFPA, GALA, FLA, CLV, CLT, IIMLA, TLA and NHLPA, are obtained.

Data Collection and Processing
The natural and social characteristics of each land-use patch should be considered when the work of land-use zoning is conducted, so we chose some indices that can reflect these characteristics and list them in Table 1. Table 1. Land-use zoning indices.

Target Criteria Indices
Land-use zoning soil quality surface soil type, soil organic matter content, pH value, soil depth, soil erosion intensity, soil profile construction topographic elevation, slope infrastructure for farmland irrigation guarantee rate, drainage conditions, rural road density land location distance from rural settlements, distance to center of towns, distance from roads land-use status land-use type Land-use type can be quantized as an integer number such as 1, 2…n.
The soil quality data such as pH value, soil depth, soil erosion intensity, etc., were obtained from the soil quality survey map (1:50,000). Slope and elevation were derived from a digital elevation model (DEM of 30 m resolution, provided by the International Scientific & Technical Data Mirror Site, Computer Network Information Center, Chinese Academy of Sciences). Some vector traffic maps (1:10,000) were used to provide the spatial information such as centers of towns, railways, highways, roads, and administrative boundaries (provided by local bureau of land and resources). The type of each land-use patch, the percentage area of each land-use type and the distribution of rural settlements were obtained from the actual vector land-use map (1:10,000, provided by the local Bureau of Land and Resources). The actual land-use map contains 47,851 patches, with each patch considered to be a land Yicheng unit which should be assigned to a certain land-use zone. Several fields were added to the actual land-use map to store the values of the aforementioned indices, after being standardized by Equation (1): where x i ' is the value of variable i before standardization, x i is the value after standardization, x min is the minimum value of the variable, and x max is the maximum. All the spatial data were processed and computed by the spatial analysis tool in ArcGIS 9.3. To ensure consistency of the data, all the images and maps were geometrically rectified with each other and subsequently referenced to the Gauss-Kruger projection Xian_1980_3_ Degree_GK_Zone_37 Projected Coordinate System. Figure 2 represents the main spatial data after GIS processing. Land-use planning based on ecological evaluation is of great significance to sustainable development [35]. The ecological benefits of land use can be evaluated by the ecosystem services value (ESV) [36], where different land-use types and zones have different ESVs. Xie et al. [37] established the ESV unit areas of Chinese terrestrial ecosystems, based on the partial global ESV evaluation results obtained by Costanza et al. [38]. This paper defines the level of ecological benefit and the corresponding coefficient of different land-use zones (see Table 2) in accordance with the foregoing results, integrated with the responses to ecological questionnaires from specialists of different departments in Yicheng City.

Objective Function
Optimized land-use zones can effectively improve the land-use multi-efficiency, achieve optimal allocation and guarantee sustainable utilization of land resources; therefore, each land-use patch should be assigned to a certain zone, not only considering the natural and social characteristics of each land-use patch but also the global spatial compactness, spatial harmony of the land-use zones and the ecological benefits of land use. Hence, optimal land-use zones should reach the following objectives: (1) The attribute difference between land-use zones maxf AD Land-use zoning needs to divide n patches x d (d = 1,2,…,m, m is the number of attributes) into C land-use zones G j (j = 1,2,…c). The difference between all the land-use zones is greatest when the sum of the distances from each land-use patch to the center of the corresponding land-use zone is smallest. So, the problem can be formulated as: where n is the number of land-use patches and c is the number of land-use zones when patch x k belongs to zone G j , then u kj is equal to 1, or otherwise equal to 0. v j is the center of zone G j that can be calculated by formula x is the value of attribute i for patch x N that belongs to zone G j , N j is the total number of patches that belongs to zone G j , and 2 k j x v  is the Euclidean distance between patch x k and v j that is given by formula (2) The spatial compactness of land-use zones maxf SC To meet the needs of land-use control, besides the optimization of the quantity and quality of land use, the optimization of the spatial pattern is also an important component of land-use zoning. Reducing the fragmentary zones benefits the spatial control. The problem can be formulated as Equation (3): If land-use patch x k is assigned to spatial cluster G j , b kj is the number of adjacent patches that also belong to G j . f SC achieves a maximum value when land-use patches with the same attributes are adjacent to each other, and f SC achieves a minimum value when each patch has no adjacent patches that belong to the same zone, a small area with several patches in Yicheng is chosen to describe this vividly ( Figure 3). Figure 3. The spatial pattern of land-use zoning: the spatial compactness of b is obviously better than a; b′ represents the spatial clusters.
(3) The degree of spatial harmony of land-use zones maxf HD While in the process of partitioning the n patches into C land-use zones and D spatial clusters (Figure 3b), the state and influence of the neighborhood of each zone (spatial cluster) should be considered. The proper adjacent distribution of two land-use zones will promote the convenience degree of social production activities of human. Therefore, the problem can be expressed as in Equation (4): where d is the number of spatial clusters, v ij is equal to 1 when cluster i and cluster j are adjacent to each other, otherwise v ij is equal to 0. hp ij represents the spatial harmony index between cluster i and cluster j, which belong to two different zones, respectively. The spatial harmony index between different land-use zones is obtained according to a combination of expert experience and the Delphi method, the results are listed in Table 3.
The ecological benefits of land-use zones maxf EB : Different land-use patterns provide a different ESV. Objective f EB aims to improve the global ESV, based on the optimal land-use zones from an ecological point of view. So, the problem can be expressed as follows: where ep k represents the ecological benefit coefficient (refer to Table 2) when patch x k is allotted to zone j. f EB achieves a maximum value when all the land-use patches are assigned to the zone of FLA, and f EB achieves a minimum value when all the patches are assigned to the zone of CLTV or IIMLA. Since these four objective functions (Equations (2-5)) represent conflicting land use demands, the scenario of land-use zoning with the optimal spatial pattern does not promise to give the maximum ecological benefits; therefore, there does not exist a single ideal solution which simultaneously satisfies the decision makers across all criteria. So, a compromise solution is sought based on a linear weighting method, which is a simple and practical way to convert the multi-objective problem into a single-objective one by using the following expression: where w 1 , w 2 , w 3 and w 4 are the weights of four objective functions, respectively, and subject to the condition w 1 + w 2 + w 3 + w 4 = 1. g norm normalizes the values of the objectives within the range [0, 1] to eliminate the impact of different dimensions by:

The stipulation of a minimum parcel area in a land-use zoning map
The standards for CLUP at the county level make clear a regulation about the minimum parcel area (MPA) of different land-use zone types in the land-use zoning map (see Table 4). The parcels (also spatial clusters) with a conflict with the regulation should be merged into an adjacent parcel which has a maximal area. In the process of land-use zoning, the quantity of land-use zones should meet the demands of the land-use plan, as well as the optimization of the spatial pattern. In accordance with the CLUP (2006-2020) for Yicheng, which states that the area of basic farmland should be not less than 62,462.14 hm 2 by 2020, and the amount of total construction land should not exceed 18,785.70 hm 2 , the forest coverage rate will simultaneously have to reach 26.36% (55,563.16 hm 2 ), so the area of BFPA, FLA, CLV, CLT and IIMLA should satisfy the inequality as follows: 62462.14 0 The general regulations of land-use zoning (1) Land-use zoning is conducted excluding the water areas. Therefore, the 527 land-use patches for the water areas will not participate in the land-use zoning and will be given a value of NONE in the attribute of land-use zone. The total area of water areas in the study area is 12,998.06 hm 2 (see Figure 4). (2) Nature reserve areas with important ecological value, special cultural relics, historical sites and ecological forest must be assigned to NHLPA as a priority, and the attribute values of 112 land-use patches for the above types remain to the very end of the process of land-use zoning, with an area of 10,111.58 hm 2 (see Figure 4).
(3) There is little likelihood of built-up areas being converted into other land-use types, based on the general laws of land-use change [39,40], so the land-use patches for urban areas and towns, and industry and mining should be only assigned to CLT and IIMLA, respectively. The same applies for tourism land, which can only be assigned to TLA.
(4) Cheng et al. [41] defined the general dominant land use, allowable land use, and prohibited land use in various land-use zones. This paper also provides the relevant regulations for the eight land-use zones of Yicheng, in line with the local conditions (see Table 5). The probability interval of a land-use patch of a certain type being assigned to a zone is determined based on the information from Table 5, and is listed in Table 6. The zone label of each land-use patch should be decided referring to the table of probability intervals, for the initialization of different scenarios of land-use zoning, as follows: where p is a random number and 0 ≤ p ≤ 1, PI kj [a,b] is the probability interval of land-use type k being allocated to land-use zone type j.
△ dominant land use, □ allowable land use，-prohibited land use.

The Particle Swarm Optimization Model
As an important branch of swarm intelligence, the PSO algorithm maintains a population of particles, where each particle represents a potential solution to an optimization problem. The set of particles, also known as a swarm, is flown through the D-dimensional search space of the problem. The position of each particle is changed, based on the experiences of the particle itself and those of its neighbors [42]. The particle in the model of land-use zoning, based on the multi-objective particle swarm optimization with constriction factor, and crossover and mutation operator (MOPSO-CCM), is seen as a potential scenario of land-use zoning. The velocity and position of the particle are updated under the guidance and restriction of the aforementioned objective functions and constraints. To implement the application of MOPSO-CCM in land-use zoning, it is necessary to define the coding scheme for the particle, design the crossover and mutation operator, and set up the updating mechanism for the particle's position and velocity, according to the particular land-use zoning problem.
(1) Coding scheme of the particle ( Figure 5) The eight land-use zones are encoded with a hybrid code system consisting of both decimal code and binary code (see Table 7). Decimal code is defined for the updating of the particle's position, and binary code is proposed for the crossover and mutation of the particle. The encoded information is stored in a relational table in the database, which is convenient for the operation of encoding and decoding.  The land-use patches are encoded with a serial integer and unique numbers in ascending order. The code of each land-use patch and corresponding land-use type are stored in a hash table named as Hash-UseType.
When the initialization of a population of random particles starts, each land-use patch obtains the type of land-use zone, based on the stochastic technique (ST) integrated with assigning the probability interval (API, refer to Table 6). The relevant information is stored in a hash table named as Hash-ZoneType.
The particle can be represented by the structure named as S_particle, which consists of a hash table (Hash-LandZoneType), the current fitness (fitness), the best fitness achieved by itself (pbest) and the best fitness obtained by any particle in the population (gbest). The array of S_particles is used to represent the swarm of particles.
(2) Designing the crossover and mutation operator ( Figure 6) In the process of crossover or mutation, the particle will be seen as a chromosome, and a land-use patch will be seen as a gene. A single-point crossover operator is adopted to reproduce offspring. The selection of parents for performing crossover is done at a specified probability and is displayed in the mating pool. The number and FID code of genes preparing for crossover are determined randomly. For the purpose of the crossover operation, the crossover point i (i∈{1,2,3}) is also chosen randomly, then the corresponding binary bits before and after the crossover point of the parents are interchanged, leading to two new individuals.
The mutation operator is designed as follows. A chromosome is chosen for mutation at a preset probability. Some genes and the corresponding position i (i {1 ∈ ,2,3}) for mutation then need to be chosen randomly. Thanks to the gene structure of the binary code, if the value at the mutation position is 1, then reset it to 0, otherwise to 1.
When crossover or mutation is proceeding, the constraint of API should be considered. For example, cropland is prohibited in the land-use zone of TLA, so if the transition of the gene offends against the constraint of API, the operation of crossover or mutation should be repeated until the transition is acceptable or reaches a preset try-number.
(3) Updating mechanism of the particle's position and velocity A new position for the particle is computed by the following formula in the basic PSO with inertial weight [43].
where x id and v id are the position and velocity of a particle, r 1 and r 2 are two distinct random values in [0, 1], c 1 and c 2 are acceleration constants, p id and p gd are the best previous position of the particle itself (pbest) and of all the particles of the swarm (gbest), respectively. Eberhart and Shi [44] found that constriction factor k introduced by Clerc [45], combined with constraints on V max (equal to x max ), significantly improved the PSO performance. The formula used to compute the velocity is then modified as: The land-use zoning problem can be regarded as an integer programming problem, which is described as follows: where Z n is the n-dimensional integer space, and S is an integer set. In this paper, S can be represented as {x|0≤ x≤ 7 and x∈Z n }. When PSO is used to solve the land-use zoning problem, the new position of the particle may be changed to a real number, even though the previous location and velocity are both integers, because the parameters of k, r 1 , r 2 , c 1 and c 2 are all real numbers. To make sure the particles are flown through the integer number space of the problem, v id (t + 1), including the three parts of kv id (t), kc 1 r 1 (p id -x id ), and kc 2 r 2 (p gd -x id ), must be an integer number in the evolutionary process. kv id (t) can be changed into an integer number by function int(). kc 1 r 1 (p id -x id ) is defined as μ 1 , when p id > x id , and μ 1 is a random number within [0, kc 1 (p id -x id )], otherwise within [-kc 1 (p id -x id ),0], so we can also define μ 1 as a random number within [a 1 , b 1 ] for simplification. On the assumption that there are n integer numbers in [a 1 , b 1 ], any integer number can be chosen to set to μ 1 , with an equivalent probability equal to 1/n. kc 2 r 2 (p gd -x id ) is defined as μ 2 with the same processing method. So, the formula to compute the velocity of the particle can be modified again as: As with the crossover and mutation, it may also offend against the constraint of API when the position of the particle is updating. Simultaneously, x id (t + 1) of the particle may be outside the range of S. So, the attribute of the land-use patch will be randomly decided, subject to the API when the phenomenon mentioned above occurred. Figure 7 presents a flow diagram of the MOPSO-CCM model, which is based on the objective functions, constraint condition, encoding and decoding scheme, and updating mechanism for the swarm. The performance of the model depends on several critical parameters [46], including the population size N, acceleration constants c 1 and c 2 , constriction factor k and maximum velocity V max . Carlisle and Doziert [47]

Comparison of the Different Scenarios
The MOPSO-CCM model adopted the crossover and mutation operator by setting the probability of mutation and crossover to 0.3 and 0.01, respectively, and the try-number to 3, and afforded several scenarios of land-use zoning by using different sets of weights applied to the four sub-objectives. The weights w n given to the various objectives were obtained according to the sub-objective weighting schemes in the study of Santé et al. [19], and are listed in Table 8.  The values of the four sub-objectives were obtained by MOPSO-CCM from the varied weighting scheme scenarios listed in Table 8 (see Table 9). The results in Table 9 show that whenever the objective of spatial compactness was included in the MOPSO-CCM objective function, along with the sub-objective of difference between land-use zones, the solution obtained displayed a significant expected increase ranging from 90.72% to 173.4% in f SC and a reduction ranging from 33.59% to 37.98% in f AD , with comparison Scenarios B-D, I-J and L-O with Scenario A. The same can be seen for the degree of spatial harmony and ecological benefit sub-objectives, with respect to Scenario A, Scenarios E-G, I and K-O increase f HD by 19.49% to 69.91%, in exchange for a reduction in f AD of 14.5% to 37.98%, and Scenarios H and J-O get a slight increase in f EB of 16.16% to 21.31%, at the expense of a reduction in f AD of 33.59% to 37.98%.
When the objective function of the MOPSO-CCM model only included the sub-objective of spatial compactness (Scenarios B-D), the value of f AD increased significantly by a factor of 2.5 as w 2 increased. The fundamental demand for attribute difference between land-use zones leads to a mass of fragmentary land-use zones; however, the sub-objective of spatial compactness can reduce the fragmentation by merging isolated land-use zones into a large contiguous land-use zone. For example, isolated CLV or GALA in a BFPA encourages merging into the BFPA. The garden land or rural settlements in BFPA will tend to be converted into farmland of high quality by consolidation of these lands in the future. Also, small BFPA or GALA showing up in a FLA encourages merging into the FLA, so the farmland with a slope greater than 25°, and a part of the farmland between 15° and 25°, should be transformed gradually into forest or grassland in the future. A typical small region is chosen from the villages of Qianfeng and Xinji to show how the sub-objective of spatial compactness optimizes the spatial pattern of land-use zones (see Figure 8). When only the sub-objective of the degree of spatial harmony degree was included, f HD had a gradual increase by a factor of more than 1.35 as w 3 increased, which can be seen from Scenarios E-G. The improvement of spatial harmony requires that the spatial pattern of BFPA is adjacent to GALA, NHLPA is contiguous to FLA, but IIMLA should be distant from BFPA, etc. A small area chosen from the village of Quanshui is presented in Figure 9 to show how the suitable neighbors are allocated to a certain land-use zone as w 3 increases. When w 1 is fixed as 0.5, the remainder of the weight is assigned equally to w 2 and w 3 (Scenario I), with a distinct decrease in f SC compared with the remainder of the weight being all assigned to w 2 (Scenario C). When it is equally assigned to w 2 and w 4 (Scenario I), in spite of f SC , it still had a slight reduction compared to Scenario C, but it was obviously greater than Scenario I. This may suggest that the pursuit of spatial harmony will bring about more fragmented land-use zones, which means worse spatial compactness, which can also be seen from the comparison between Scenarios E-G and graphically in Figure 9. An increase in the ecological benefits can also improve the spatial compactness. From Scenarios A-D, we can see that f HD reduced as f SC increased, which proved that the spatial compactness also has a negative effect on spatial harmony. Another comparison of Scenario A with Scenarios H and J shows that f HD increased by a factor of about 1.1 as f EB increased, which suggests that the increase in ecological benefits will boost the degree of spatial harmony of the land-use zones.
From the comparison of Scenario L with Scenarios M-O, it can be seen that assigning more weight to a certain sub-objective, and the same weight to other sub-objectives, had the expected effect of increasing the value of the corresponding sub-objective, indicating that the fitness values are positively correlated with the weights. That is to say, a decision maker for land-use planning can obtain the Pareto solution to satisfy their subjective preference by giving a higher weight to the corresponding sub-objective. Figure 10 presents the area of each land-use zone determined in the aforementioned four scenarios. It can be seen from Figure 10 that: (1) NHLPA in each scenario have the same area, because of the priority of being fixed, based on the general regulations for land-use zoning. The areas of BFPA and FLA are greater than 62,462.14 hm 2 and 55,563.16 hm 2 , respectively, in each scenario, and the sum of the areas of CLV, CLT and IIMLA is less than 18,785.7 hm 2 , which indicated the efficiency of the constraint of area limitation of the various land-use zones; (2) Scenario M got the largest area of BFPA, CLT and IIMLA and the smallest area of FLA, CLV and TLA amongst the three scenarios. With the drive for spatial compactness, the large spatial clusters annexed the small clusters surrounding them to generate an optimal spatial pattern, simultaneously resulting in either a larger land-use zone or a smaller one; (3) Scenario N had the largest land-use zone of CLV and the second largest area of BFPA and GALA. This is because the spatial harmony between these land-use zones is higher than with the others; (4) The largest FLA was obtained in Scenario O because the forests can provide the greatest ESV. Meanwhile, Scenario O got the smallest BFPA and GALA. This suggests that the enforcement of the policy of 'grain for green' can improve the ESV of the land-use system; (5) Scenario L, assigning the same weight to each objective function, provided a compromise solution with a moderate value for each objective function (see Table 9) and area for each land-use zone.
The global spatial distribution of the eight land-use zones in the four Scenarios of L, M, N, O are similar to each other, with only some microscopic differences. The reason for these disparities has been described above, with a vivid depiction in Figures 8 and 9. As far as Scenario L is concerned, the quantity structure and spatial pattern of the land-use zones are presented in Figure 11.
The BFPA is distributed mainly in the towns of Xiaohe, Zhengji, Kongwan, Wangji and Nanying, and accounts for 36.25% of the total area, excluding the water areas. A total of 85.62% of the fifth level and 63.25% of the fourth level of agricultural land in the study area are assigned to BFPA. (The quality of agricultural land in Yicheng is classified into five levels: first level-the worst quality, second level, third level, fourth level, and fifth level-the highest quality.) This zone is marked off to hold its ability to continue feeding a growing population and meet the food security goals. High-quality land with high productivity in this zone should be prohibited from converting to non-agricultural uses. The conversion of good-quality land with moderate productivity to non-agricultural uses is permitted by law under some circumstances, usually after a planned period of 5-10 years [48]. GALA is observed mainly in the towns of Wangji, Nanying, Banqiao and Liushui, and is contiguous with BFPA. This zone primarily contains the land for gardens, aquaculture, etc., with high ecological value, and the rural roads, infrastructural facilities, water conservancy facilities, protective forest, etc., serve the farmland. In this region, the land consolidation project is encouraged to make the 'dynamic balance of farmland' policy materialize. FLA is concentrated to the east and southwest of Yicheng and mainly consists of fast-growing and high-yield broad-leaf forests. In view of the functions of forests, such as soil and water conservation, and ecological protection, the conversion of forest land is severely restricted, except for approval by the competent forestry authorities. The spatial clusters of CLV are evenly dispersed across the maps, the surroundings of which are mostly of BFPA and GALA, which annex some isolated rural settlements. At the same time, some farmland patches surrounding rural settlements are merged into them. CLT is marked off, based on actual land use for urban areas and towns, including some agricultural land, which can be used as the preparatory land for the future expansion of the city. Small villages, 'hollow-village' and 'villages-in-the-city' [49] within CLV and CLT should be amalgamated so as to upgrade the scale and organization of villages and towns. The IIMLA only covers around 0.34% of the territory of the city, and mainly comprises some land for the mining and chemical industries. The land that was destroyed in the process of industrial production and construction should be reclaimed for new use, especially cultivated land, if the land is suitable for agriculture. TLA is located in the vicinity of the region's water areas, and NHLPA is distributed in the west, southeast and northeast corner of the city. In these regions, a small amount of agricultural land is permitted to improve the biodiversity and maintain the ecological balance, and the land for enterprises with industrial pollution is strictly prohibited.

Analysis of the Effect of Improvements for the Model
The focal point of this subsection is to analyze the effect of the adoption of the crossover and mutation operator on acquiring the optimal solutions. A comparative analysis of the optimal fitness values, convergence speed and run time is made to evaluate the performance of the MOPSO-CCM model and the basic particle swarm optimization with only a constriction factor (MOPSO-BC), in Scenarios I, J, K and L.
It is observed in Table 10 that the optimal solution obtained by the MOPSO-CCM model exhibits the expected higher fitness value for each sub-objective in each scenario, with f AD increasing by 15.6% to 17.48%, f SC increasing by 14.61% to 16.89%, f HD increasing by 9.6% to 14.18%, and f EB increasing by 1.26% to 2.38%. Figure 12 clearly illustrates that the weighted sum fitness obtained by MOPSO-CCM in each generation is much higher than that of MOPSO-BC. At the same time, the convergence speed of MOPSO-CCM, with a sudden rise in its curve, is distinctly higher than that of MOPSO-BC, even if the initial normalized fitness value of MOPSO-CCM is lower than MOPSO-BC in Scenario K. The sudden rises in the curves of MOPSO-CCM in the later period of evolution show that the crossover and mutation operator helps the model to escape the local optimum and achieve a better solution. In summary, the convergence and exploration ability of the MOPSO-CCM model outperforms MOPSO-BC, which proves the validity of improvement for the PSO algorithm with land-use zoning. It is necessary, however, to note that MOPSO-CCM requires some computational cost to perform the operation of crossover and mutation on particles, based on a preset probability, so it is more time consuming than the MOPSO-BC model.
A discussion of the stability of the improved PSO model is conducted next. We ran the MOPSO-CCM and MOPSO-BC model 100 times, respectively, with weighted Scenario L, and calculated the mean value and standard deviation of the four sub-objectives in each optimal solution (see Table 11). Although the mean value of MOPSO-CCM is obviously higher than MOPSO-BC, what further proves the improved exploration ability of MOPSO-CCM is that the standard deviation of MOPSO-BC is lower than MOPSO-CCM, which indicates that MOPSO-CCM is not as stable as MOPSO-BC. The optimal fitness values of each sub-objective are plotted in the curves. The curves of MOPSO-CCM in Figure 13 show a greater fluctuation than MOPSO-BC. The statistics table and graph reveal that the adoption of the crossover and mutation operator remarkably increases the probability of acquisition of an optimal solution by improving the diversity of the particle swarm, but also brings about some disturbance to the exploitation and exploration ability of the model.

Conclusions
The traditional approaches to the land-use zoning problem in China just take into consideration the characteristic attributes of the land units. In this article, the land-use zoning problem is treated as a MOP. Four objective functions were formulated, respectively, for attribute differences between land-use zones, spatial compactness, spatial harmony and the ecological benefits of land-use zones. A PSO-based model equipped with a crossover and mutation operator was constructed. The feasibility of the proposed approach in solving the land-use zoning problem was checked using a case study in Yicheng, China. The results obtained indicated that the integration of GIS and MOPSO-CCM is a promising and efficient approach for solving the land-use zoning problem. The inclusion of spatial compactness and spatial harmony in the objective functions allows the achievement of solutions with an optimal spatial pattern, and the inclusion of ecological benefits leads to a greater quantity of land-use zones with a high ESV. It is worth pointing out that the objective function can be extended in line with local conditions, such as the economic benefits for a county with a backward economy, spatial contiguity, etc. In this paper, a conventional weighted aggregation method is used for the MOPSO-CCM model to search for the Pareto-dominant solutions, the weakness of which is that an optimal solution can be obtained with only one run of the model. The next stage is to apply a dynamic weighted aggregation (DWA) technique to the MOPSO-CCM model to generate a series of Pareto-optimal solutions.
The incorporation of a crossover and mutation operator proved to be effective in improving the convergence and exploration ability, at the price of a relatively small increase in runtime and, simultaneously, a poorer stability than MOPSO-BC. As we know, the performance of PSO depends on a few different parameters [50]. Therefore, a further study will focus on choosing a group of parameters that will enable the MOPSO-CCM model to achieve a better performance with the land-use zoning problem.