Land Use Optimization for Coastal Urban Agglomerations Based on Economic and Ecological Gravitational Linkages and Accessibility

: Urban agglomerations (UA) are attracting increasing research attention as a global emergent phenomenon, whereby regional collaborative linkages between cities attracts and agglomerates development. However, these studies also acknowledge that ecological values may be negatively impacted by re-development, ecological fragmentation, and proximity or downstream impacts. Sus-tainable development, therefore, requires balancing forces from economic attraction and ecological repulsion. Forces similar to economic ones may also operate in attracting ecological enhancement towards higher-valued ecological regions; however, research regarding the role of the self-collaborative gravity-like forces shaping UA is limited in land use optimization. To assist planners, this study developed a new multi-objective land use optimization of UA that explored the intensity of economic ties and ecological gradients using the multi-objective NSGA-II algorithm. In this model, economic linkage intensity (ELI) and accessibility were used to calculate a modiﬁed GDP (gross domestic product), while the NDVI (normalized difference vegetation index) was used for the modiﬁed ESV (ecosystem services value). Spatial allocation with implicit economic accessibility relationships was enhanced through a two-step mutation operator, including a “gravity ﬂip” spatial orientation factor. Compared to the standard NSGA-II algorithm, models of future land use of the Guangdong-Hong Kong-Macao Greater Bay Area (GBA) in 2030 have shown that the modiﬁed GDP value in our model increased by 7.41%, while the conversion rate of high-density vegetation reduced by 7.92%. The results highlighted the importance of linkage and accessibility factors in enhancing the clustering of cities. In tandem, the modiﬁed ESV also enhances ecosystem services contributions of higher value vegetated land through decentralized built-up developments. The proposed model provides managers with a comprehensive and efﬁcient land use solution model that accounts for intrinsic linkage factors shaping the development of compact urban agglomerations.


Introduction
The emergence of urban agglomerations (UA) globally is thought to improve economic cooperation, livability, and transportation, especially along coastal areas [1][2][3][4].There are more than 20,000 UA in the 151 countries, and research has proved that the urban growth of UA in the 21st century will produce increasingly concentrated cities with high population density [5].Expansion of these UA is driven not only by the speed of economic development and population growth, but also through collaborative economic linkages between cities, and regional collaborative development has become even more important results according to the planner's needs [39].The non-dominated sorting genetic algorithm II (NSGA-II) is a classic metaheuristic algorithm based on GA, and it can obtain Pareto set by non-dominated sorting; it has been widely used for optimizations, due to its good robustness and global search capability [40][41][42].For example, it has been successfully applied to land use planning and integrated optimization using system dynamics models [43][44][45].However, the complexity of NSGA-II problems often leads to long calculation times, especially in land use optimization.Significant time reductions can be achieved with computationally parallelized NSGA-II [46,47].Even small computational speedup across numerous repetitive calculations in NSGA-II can provide substantial overall benefits, and parallel processing is a reliable method for achieving computational improvements.
Most previous urban land use optimization studies considered the conflicting relationship between economy and ecology but ignored the role of ecological gradients and economic linkage between cities [13].This paper assumes that the economic linkages and ecological gradients can play a positive role in the future regional optimization.To simplify the model and emphasize the collaborative development demands of UA, we do not consider the effects of factors such as food security [16,48], environmental deterioration [49], or climate risk [50] on land use optimization and propose a new multi-objective land use optimization of UA that explored the intensity of economic ties and ecological gradients using the multi-objective NSGA-II algorithm.In the rest of this paper, a description of the land use optimization model, consisting of three conflicting objectives, is described first, followed by details of the parallel NSGA-II algorithm used.Subsequently, the model is applied to the Guangdong-Hong Kong-Macao Greater Bay Area (GBA) to assess its agglomeration capabilities, and the results and conclusions are discussed.

Materials
Our case study is the Guangdong-Hong Kong-Macao Greater Bay Area (GBA), situated in the Pearl River Delta in China, which includes nine mainland cities and two special administrative regions, namely Guangzhou (GZ), Shenzhen (SZ), Zhuhai (ZH), Zhaoqing (ZQ), Huizhou (HZ), Jiangmen (JM), Dongguan (DG), Foshan (FS), Zhongshan (ZS), Hongkong (HK), and Macao respectively (Figure 1).Covering a total area of 64,823 km 2 , and with a total population of over 66 million at the end of 2015, it is one of the most economically dynamic regions in China.As of 2015, the total GDP reached 6745 billion CNY, where the economic volume and population density are mainly concentrated in HK, SZ, and GZ, and the regional land use shows a spatial gradient of decreasing built-up development from the bay area to the hinterland.
Economic and demographic data were obtained from the statistical yearbooks of each city, which contained the total GDP of primary, secondary, and tertiary industries (Table 1).The ecological control areas data of the GBA were obtained from the website of the Ministry of Natural Resources of the People's Republic of China (http://g.mnr.gov.cn/(accessed on 21 February 2022)).Land use data (Figure 1), and NDVI data (Figure 2a) were obtained from the Resource and Environmental Science and Data Center, with a resolution of 1 × 1 km (https://www.resdc.cn/(accessed on 5 December 2021)).The intensity of urban economic linkage data were obtained from Peng's paper [51], while the road network was obtained from Open Street Map (https://www.openstreetmap.orgaccessed on 2 December 2021).DEM data (Figure 2b) were obtained from the Shuttle Radar Topography Mission (SRTM, http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp (accessed on 2 December 2021)).Economic and demographic data were obtained from the statistical yearbooks of city, which contained the total GDP of primary, secondary, and tertiary industries (T 1).The ecological control areas data of the GBA were obtained from the website o Ministry of Natural Resources of the People's Republic of China (http://g.mnr.gov.cncessed on 21 February 2022)).Land use data (Figure 1), and NDVI data (Figure 2a) obtained from the Resource and Environmental Science and Data Center, with a re tion of 1 × 1 km (https://www.resdc.cn/(accessed on 5 December 2021)).The intens urban economic linkage data were obtained from Peng's paper [51], while the road work was obtained from Open Street Map (https://www.openstreetmap.orgaccesse 2 December 2021).DEM data (Figure 2b) were obtained from the Shuttle Radar Top phy Mission (SRTM, http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp (accessed December 2021)).
Based on the land use data from 1995-2015 (five years for one period), linear re sion was used to predict the total built-up land area in 2030 to be 12,371 km 2 .To en that the future agricultural land conforms to government planning requirements, we the data from the 2010-2020 master plan of Guangdong Province for the agricultural  Based on the land use data from 1995-2015 (five years for one period), linear regression was used to predict the total built-up land area in 2030 to be 12,371 km 2 .To ensure that the future agricultural land conforms to government planning requirements, we used the data from the 2010-2020 master plan of Guangdong Province for the agricultural land retention projection to estimate agricultural land retention of 8299 km 2 in 2030, excluding Hong Kong and Macao.

Method
The steps in applying the NSGA-II optimization algorithm are shown in Figure 3.The first step was to derive future land use area constraints for the optimization case study by predicting the future land use structure using a linear model of past trends.Using NSGA-II, we incorporated accessibility and ELI as modification factors and used parallel computational processing to speed up the model.

Method
The steps in applying the NSGA-II optimization algorithm are shown in Figure 3.The first step was to derive future land use area constraints for the optimization case study by predicting the future land use structure using a linear model of past trends.Using NSGA-II, we incorporated accessibility and ELI as modification factors and used parallel computational processing to speed up the model.

Method
The steps in applying the NSGA-II optimization algorithm are shown in Figure 3.The first step was to derive future land use area constraints for the optimization case study by predicting the future land use structure using a linear model of past trends.Using NSGA-II, we incorporated accessibility and ELI as modification factors and used parallel computational processing to speed up the model.We assumed that the study area was a two-dimensional regular grid (rasterized region).The land use allocation problem was to determine how to assign future (changed) land uses to specific locations, so that the optimal land use pattern balanced the conflicting objectives, based on Pareto front solutions.
To determine constrains for the overall areas of future land use, a linear model was used to extrapolate past trends of land use.In the case study area (described below), we used data of land use from every five years of the past twenty years to forecast land use constraints 15 years into the future.The modified NSGA-II algorithm was then used to determine Pareto front solutions using modified objectives and genetic algorithm procedures described in the next section.

Objectives (a) Modified GDP
GDP is widely used to measure the level of regional development, and many scholars have adopted this objective in land use optimization [13,20,21,52].Since different land uses correspond to different industries, the study uses the sector-based classification designed by Cao to calculate the GDP per unit area for different land uses [13], as shown in Table 2. Uneven development reflects complex linkages between, and drivers of, regional and local policies, trade, and market.Rapid urbanization in UA can lead to uneven development, and the strength of economic linkages between cities reflects the potential for regional economic collaborative development, which is closely related to distance, population size, and the level of economic development [1].The greater the economic linkage intensity (ELI) between cities, the greater the region's capacity for collaborative development.ELI is, therefore, a potential factor that drives change and integrated development in UA.In our study, we tested its incorporation as a spatial modification factor for GDP, in order to guide the direction of collaborative land use change towards regions of higher ELI.
Accessibility, in relation to regional commuting time and reduction in environmental impacts, such as CO 2 emissions, has significant impact on land use status and socioeconomic development [53].Within a region of several cities, accessibility is a core issue affecting connectivity, livability, and the shaping of UA [31].Accessibility was, therefore, used as a spatial modification factor for GDP, calculated according to the formula: where ELI k represents the ELI value in cell k.GDP k and Acc k represent the unit GDP value and accessibility value in cell k, respectively.K represent the total number of changed cells.
Factors included in calculating ELI were: total population, GDP, total import and export, and tourism income, based on the gravity model of urban economic linkage data obtained from Peng's paper [51].ELI data can present the economic influence degree between each city, which have high spatial autocorrelation.Considering that ELI is a function of city quality and spatial distance between cities, which decays exponentially with increasing distance.Kringing is the method use the spatial correlation between sampled point data to interpolate the value in the spatial field.Compared with other interpolation methods, Kringing can preserve spatial variability and generate estimates of uncertainty surrounding interpolated value.Thus, in this paper, Kringing was used in this paper to get the spatial distribution map of ELI.

(b) Modified ESV
ESV is a main factor determining the livability, attractiveness, and ecological health of a region.Obtaining greater economic value with minimal ESV loss is often one of the main goals of land use management.Conventionally, the same land use type is regarded as having the same ESV value, thereby ignoring differences in ESV due, for example, to vegetation densities, which could be used to temper the direction of urban expansion.Therefore, in this paper, the NDVI was used as modified index to identify the vegetation density of the selecting cells.With higher NDVI, there will be a lower chance of change.The modified ESV objective change was calculated by multiplying the unit ESV by the NDVI, using the formula: where k represents the changed cell.K represent the total number of changed cells.NDVI k is the value of NDVI in cell k.ESV k is the unit ESV in cell k.

(c) Compactness
Compactness is a measure of urban intensification and land use efficiency [14,54].For any cell, compactness is the number of cells with the same land use type in its eight neighborhoods, i.e., the more cells with the same value in its neighborhood, the higher the compactness.The formula is as follows: where N i and N i are the number of cells with the same land use type in i cell's eight neighborhoods in the solution and initial pattern, respectively.I represents the total number of cells in the region.

Constraints
Constraints on the regional optimization, detailed later in the case study, are set by the ecological control zone, slope constraint, total agricultural land area constraint, and built-up land area constraint.Additionally, each cell's land use type is unique and, therefore, cannot be converted to another land use type.

Parallel NSGA-II
The NSGA-II algorithm, widely used in land use optimization, was used to obtain the Pareto solution set from multiple non-linear objectives through non-dominated sorting.As shown in Figure 1, improved NSGA-II involved five steps: initialization, non-dominated sorting, selection, crossover, accessibility mutation, and gravity flip.These steps are described below, bearing in mind that a land use pattern represents a chromosome, and each land use cell is called a gene.To reflect economic attraction between cities shaping UA, we modified the algorithm with a two-step mutation operator and improved the overall computational efficiency through parallel processing.
Initialization: Based on the projected area of each land use in future (2030 in the case study), the land use that needs to be increased is randomly assigned to other land uses, and 100 initial chromosomes are set to form the initial population.
Non-dominated sorting and crowing distance: A selected solution in the search space is a non-dominated solution [55] only if there are no other solutions that are equal in all objectives or better than that solution in at least one objective.Crowding distance is used to choose the better solution on a front shown in Figure 4, which represents the non-dominated sorting of NSGA-II, based on the following process.Assume that the parent is Pi and population after crossover and mutation is T.There are M solutions within the population.First, P and T are combined and then sorted non-dominantly from shortest to longest crowding distance, and the top M chromosomes are selected as the new offspring P i+1 .
The NSGA-II algorithm, widely used in land use optimization, was used to obtain the Pareto solution set from multiple non-linear objectives through non-dominated sorting.As shown in Figure 1, improved NSGA-II involved five steps: initialization, non-dominated sorting, selection, crossover, accessibility mutation, and gravity flip.These steps are described below, bearing in mind that a land use pattern represents a chromosome, and each land use cell is called a gene.To reflect economic attraction between cities shaping UA, we modified the algorithm with a two-step mutation operator and improved the overall computational efficiency through parallel processing.
Initialization: Based on the projected area of each land use in future (2030 in the case study), the land use that needs to be increased is randomly assigned to other land uses, and 100 initial chromosomes are set to form the initial population.
Non-dominated sorting and crowing distance: A selected solution in the search space is a non-dominated solution [55] only if there are no other solutions that are equal in all objectives or better than that solution in at least one objective.Crowding distance is used to choose the better solution on a front shown in Figure 4, which represents the non-dominated sorting of NSGA-II, based on the following process.Assume that the parent is Pi and population after crossover and mutation is T.There are M solutions within the population.First, P and T are combined and then sorted non-dominantly from shortest to longest crowding distance, and the top M chromosomes are selected as the new offspring Pi+1.Selection: After non-dominated sorting of the parent population, roulette is performed, based on the crowding distance.The smaller the distance, the higher the probability of being selected.
Crossover: Crossover operator is used to preserve relatively better land-use change patterns under multi-objective constraints and reduce the generation of duplicate individuals.In original NSGA-II, the selected cell in the two parents will be exchanged to form the offspring.Different to the original NSGA-II, the selected cells in our model use maximum count of neighbors' land use to choose whether to change.The procedure for crossover is as follows: after determining the crossover probability, if the number of cells amongst the eight neighbors that are the same as the central cell in Parent1 is less than that in Parent2, crossover is performed by replacing the center cell in Parent1 with the central cell land use in Parent2, as illustrated in Figure 5. Selection: After non-dominated sorting of the parent population, roulette is performed, based on the crowding distance.The smaller the distance, the higher the probability of being selected.
Crossover: Crossover operator is used to preserve relatively better land-use change patterns under multi-objective constraints and reduce the generation of duplicate individuals.In original NSGA-II, the selected cell in the two parents will be exchanged to form the offspring.Different to the original NSGA-II, the selected cells in our model use maximum count of neighbors' land use to choose whether to change.The procedure for crossover is as follows: after determining the crossover probability, if the number of cells amongst the eight neighbors that are the same as the central cell in Parent1 is less than that in Parent2, crossover is performed by replacing the center cell in Parent1 with the central cell land use in Parent2, as illustrated in Figure 5.
Mutation operator is used to enhance population diversity and prevent the optimization results from falling into a local optimum trap.Original NSGA-II only has one mutation operator.For the selected cell, the land use type of the selected cell is replaced by calculating the highest number of land use types in its eight neighborhoods.Unlike the original NSGA-II, two mutations with implicit economic relationships were set for enabling agglomeration and limiting fragmentation in our model, the procedure for two-step mutation operator is as follows.
Accessibility Mutation: We developed a new mutation operator that used accessibility as a modification factor.Referring to the example in Figure 6, after determining the mutation probability, cell AA (with land use type of built-up land) was selected, and its accessibility value was 0.2, while cell BB was random selected, and its accessibility value was 0.4.Since the accessibility value of BB was higher than AA, the land use type shifted between AA and BB.Mutation operator is used to enhance population diversity and prevent the optimization results from falling into a local optimum trap.Original NSGA-II only has one mutation operator.For the selected cell, the land use type of the selected cell is replaced by calculating the highest number of land use types in its eight neighborhoods.Unlike the original NSGA-II, two mutations with implicit economic relationships were set for enabling agglomeration and limiting fragmentation in our model, the procedure for two-step mutation operator is as follows.
Accessibility Mutation: We developed a new mutation operator that used accessibility as a modification factor.Referring to the example in Figure 6, after determining the mutation probability, cell AA (with land use type of built-up land) was selected, and its accessibility value was 0.2, while cell BB was random selected, and its accessibility value was 0.4.Since the accessibility value of BB was higher than AA, the land use type shifted between AA and BB.Gravity flip: ELI is calculated from the gravity model, and its spatial distribution has a structure similar to a gravity field.Areas with strong gravity attract areas from weak gravity.The spatial linkage of cities can be optimized by gravity flip, which is used to determine if spatial flipping results in an overall stronger impact of the gravity field on modified GDP.As a simple example, where weight = mass × gravity, moving a mass closer to a stronger gravity field results in greater weight.Likewise, moving an area producing GDP closer to a stronger ELI field results in potentially higher GDP generated = original GDP × ELI.Using 3 × 3 window as unit, the total modified GDP in the original window and solution after gravity flip, which, in this case, corresponds to a transposition, are compared.If the total modified GPD after flipping is higher, then the window will be flipped (Figure 7).Note that vertical and horizontal flipping are possible as options.Mutation operator is used to enhance population diversity and prevent the optimization results from falling into a local optimum trap.Original NSGA-II only has one mutation operator.For the selected cell, the land use type of the selected cell is replaced by calculating the highest number of land use types in its eight neighborhoods.Unlike the original NSGA-II, two mutations with implicit economic relationships were set for enabling agglomeration and limiting fragmentation in our model, the procedure for two-step mutation operator is as follows.
Accessibility Mutation: We developed a new mutation operator that used accessibility as a modification factor.Referring to the example in Figure 6, after determining the mutation probability, cell AA (with land use type of built-up land) was selected, and its accessibility value was 0.2, while cell BB was random selected, and its accessibility value was 0.4.Since the accessibility value of BB was higher than AA, the land use type shifted between AA and BB.Gravity flip: ELI is calculated from the gravity model, and its spatial distribution has a structure similar to a gravity field.Areas with strong gravity attract areas from weak gravity.The spatial linkage of cities can be optimized by gravity flip, which is used to determine if spatial flipping results in an overall stronger impact of the gravity field on modified GDP.As a simple example, where weight = mass × gravity, moving a mass closer to a stronger gravity field results in greater weight.Likewise, moving an area producing GDP closer to a stronger ELI field results in potentially higher GDP generated = original GDP × ELI.Using 3 × 3 window as unit, the total modified GDP in the original window and solution after gravity flip, which, in this case, corresponds to a transposition, are compared.If the total modified GPD after flipping is higher, then the window will be flipped (Figure 7).Note that vertical and horizontal flipping are possible as options.Gravity flip: ELI is calculated from the gravity model, and its spatial distribution has a structure similar to a gravity field.Areas with strong gravity attract areas from weak gravity.The spatial linkage of cities can be optimized by gravity flip, which is used to determine if spatial flipping results in an overall stronger impact of the gravity field on modified GDP.As a simple example, where weight = mass × gravity, moving a mass closer to a stronger gravity field results in greater weight.Likewise, moving an area producing GDP closer to a stronger ELI field results in potentially higher GDP generated = original GDP × ELI.Using 3 × 3 window as unit, the total modified GDP in the original window and solution after gravity flip, which, in this case, corresponds to a transposition, are compared.If the total modified GPD after flipping is higher, then the window will be flipped (Figure 7).Note that vertical and horizontal flipping are possible as options.To summarize, accessibility mutation determines whether the location of a random built-up land cell and another randomly selected non-built-up land cell are exchanged, based on the accessibility value, while gravity flip is implemented if the total modified GDP is higher than the original value when the cells in the 3 × 3 window are rotated 180 degrees clockwise.To summarize, accessibility mutation determines whether the location of a random builtup land cell and another randomly selected non-built-up land cell are exchanged, based on the accessibility value, while gravity flip is implemented if the total modified GDP is higher than the original value when the cells in the 3 × 3 window are rotated 180 degrees clockwise.

Objective Quantification and Constraints
The spatial constraints were set, based on following restrictions: (1) the total arable land area to be no less than 8299 km 2 and according to the overall development plan of each city; (2) all changeable cells could not be changed to unused land; (3) slopes less than 15 degrees can be built-up land, less than 25 degrees can be cropland, and more than 25 degrees can only be forest; (4) all ecological control area and water could not be changed (Figure 8).All the restrictions are not only applied in the initialization process, but also play an important role in constraining the global situation.

Objective Quantification and Constraints
The spatial constraints were set, based on following restrictions: (1) the total arable land area to be no less than 8299 km 2 and according to the overall development plan of each city; (2) all changeable cells could not be changed to unused land; (3) slopes less than 15 degrees can be built-up land, less than 25 degrees can be cropland, and more than 25 degrees can only be forest; (4) all ecological control area and water could not be changed (Figure 8).All the restrictions are not only applied in the initialization process, but also an important role in constraining the global situation.Parameters were set as follows: the unit ESV of each land use is derived from the results of Xie's assessment, as shown in Table 3 [56].GDP is calculated, with reference to Cao's study, and the unit GDP values of each land use are shown in Table 3 [13].The initial spatial distributions of GDP and ESV (Figure 9) were obtained using ArcGIS.This paper uses Peng's results of ELI for the Greater Bay Area [51].Accessibility is expressed by road network density, and its spatial distribution was obtained in arcGIS (Figure 10) after applying the normalization method in Pan's model [21].Parameters were set as follows: the unit ESV of each land use is derived from the results of Xie's assessment, as shown in Table 3 [56].GDP is calculated, with reference to Cao's study, and the unit GDP values of each land use are shown in Table 3 [13].The initial spatial distributions of GDP and ESV (Figure 9) were obtained using ArcGIS.This paper uses Peng's results of ELI for the Greater Bay Area [51].Accessibility is expressed by road network density, and its spatial distribution was obtained in arcGIS (Figure 10) after applying the normalization method in Pan's model [21].

Implementation and Evaluation
A total of 100 random patterns were used as the initial population, based on probabilities of crossover, accessibility mutation, and gravity flip of 0.06, 0.1, and 0.1, respectively.A total of 600 iterations simulated, in order to optimize the three conflicting objectives under constraints on the agricultural land total and ecological control areas.We found that parallel NSGA-II outperformed classical NSGA-II in execution, taking 16.71 s in each iteration, while classical NSGA-II took 151.63 s.

Implementation and Evaluation
A total of 100 random patterns were used as the initial population, based on probabilities of crossover, accessibility mutation, and gravity flip of 0.06, 0.1, and 0.1, respectively.A total of 600 iterations were simulated, in order to optimize the three conflicting objectives under constraints on the agricultural land total and ecological control areas.We found that parallel NSGA-II outperformed classical NSGA-II in execution, taking 16.71 s in each iteration, while classical NSGA-II took 151.63 s.

Implementation and Evaluation
A total of 100 random patterns were used as the initial population, based on probabilities of crossover, accessibility mutation, and gravity flip of 0.06, 0.1, and 0.1, respectively.A total of 600 iterations were simulated, in order to optimize the three conflicting objectives under constraints on the agricultural land total and ecological control areas.We found that parallel NSGA-II outperformed classical NSGA-II in execution, taking 16.71 s in each iteration, while classical NSGA-II took 151.63 s.
The progress of iterations (Figure 11a) shows that the solutions achieved a better spread across the three objectives, with the respect to the overall optimization, but improvements become smaller and smaller.To illustrate the convergence trend in the model, we use minimum modified ESV, shown in Figure 11b, which showed a rapid decrease up to 90 iterations and then an almost static phase through successive iterations.The curve appears stable after 212 iterations, implying the solution process efficiently approaches near-optimality relatively quickly.
spread across the three objectives, with the respect to the overall optimization, but improvements become smaller and smaller.To illustrate the convergence trend in the model, we use minimum modified ESV, shown in Figure 11b, which showed a rapid decrease up to 90 iterations and then an almost static phase through successive iterations.The curve appears stable after 212 iterations, implying the solution process efficiently approaches near-optimality relatively quickly.In Figure 11a, modified ESV and compactness both change negatively to modified GDP, and the pattern with smallest compactness loss also has smallest modified ESV loss.Taking each preferred solution as an example of the optimal solution set, the following specific solutions are possible: Solution I (SI)-the pattern with the highest modified GDP; Solution II (SII)-the pattern with minimal modified ESV loss; Solution III (SIII)-the pattern with minimal compactness loss.To demonstrate the superiority of our model, solutions from the classic NSGA-II approach were also simulated.In this approach, the GDP and ESV were not affected by their spatial location, and cells selected in mutation changed their land use values randomly.Solution IV (SIV)-was obtained from the classic NSGA-II for comparison with SI, due to the almost identical increments of build-up land, ESV loss, and compactness loss.The objective values in SIV were recalculated according to our proposed method to facilitate the comparison.The comparative results from the different solutions are shown below (Figure 12, Tables 4 and 5).We did not count the percentage change in unused land because the total area of unused land is almost zero.For comparison purposes, the spatial distribution of land use and built-up land increment for each scenario are shown in Figure 12.In Figure 11a, modified ESV and compactness both change negatively to modified GDP, and the pattern with smallest compactness loss also has smallest modified ESV loss.Taking each preferred solution as an example of the optimal solution set, the following specific solutions are possible: Solution I (SI)-the pattern with the highest modified GDP; Solution II (SII)-the pattern with minimal modified ESV loss; Solution III (SIII)-the pattern with minimal compactness loss.To demonstrate the superiority of our model, solutions from the classic NSGA-II approach were also simulated.In this approach, the GDP and ESV were not affected by their spatial location, and cells selected in mutation changed their land use values randomly.Solution IV (SIV)-was obtained from the classic NSGA-II for comparison with SI, due to the almost identical increments of build-up land, ESV loss, and compactness loss.The objective values in SIV were recalculated according to our proposed method to facilitate the comparison.The comparative results from the different solutions are shown below (Figure 12, Tables 4 and 5).We did not count the percentage change in unused land because the total area of unused land is almost zero.For comparison purposes, the spatial distribution of land use and built-up land increment for each scenario are shown in Figure 12.

Discussion
In the past decades, many scholars have conducted regional land use optimization studies and focused on multi-objective trade-offs, as well as model improvements [20,[57][58][59]; however, few studies take the effect of inner-city economic linkages and ecological gradient on land use changes in UA into account [13,60].Cao have proposed land use  SI (GDP), SII (ESV), and SIII (compactness) achieved better scores, with respect to their preferred objectives.In SI (GDP), built-up land increased by 32.88%, forest and cropland decreased 8.84% and 1.36%, respectively, while grassland increased by 0.34%.Compared to SII (ESV), more built-up land increment in SI (GDP) led to a greater decline in modified ESV and compactness.In SII (ESV), although most of the built-up land increment comes from forest, grassland is used more for land use transfer thereby reducing the amount of forest transfer, in a way that minimizes the modified ESV reduction while built-up land increased.Almost identical ESV and compactness value in SI and SIV with the same built-up land increment, the modified GDP value in SI was significantly higher than SIV (7.41%), which illustrates the superiority of our model in getting better pattern with higher economic output.
Comparing the spatial distribution of built-up land increments across scenarios, the inclusion of ELI, accessibility and NDVI as modifying factors appears to positively affect the trend of urban expansion.The cluster of built-up land increments in SII (ESV) is worse than SI (GDP), reflecting that modified ESV increased for more decentralized patterns of built-up land increments.Although SIII (Compactness) have the least loss of compactness, the cluster of built-up land increments is worse than SI (GDP), reflecting the limited role of compactness objective on the distribution of built-up land clustering.In addition, the built-up land increment in SIV is spatially evenly distributed without clustering, this phenomenon arises because spatial agglomeration can only be controlled through the compactness objective, without influences from gravitation attraction of ELI, accessibility, and ESV.
The ELI increments were 1072 and 993 for SI and SIV and the accessibility increments were 773 and 565 for SI and SIV, respectively; implying that SI has higher ELI and accessibility with the same built-up land increment.Compared with SIV, the spatial distribution of built-up land increment is more clustered in SI, reflecting the role of modified GDP in spatial agglomeration proposed by our model.To explore the advantages of the modified ESV, we divided NDVI into three level by natural break method, and analyzed the amount of vegetation transfer at different levels in SI and SIV (Table 6).The total changed percent of vegetation (cropland, grassland, and forest) of high level in SIV was 7.92% higher than SI, indicating the effectiveness of our model in optimizing selection of low value vegetation for changes.

Discussion
In the past decades, many scholars have conducted regional land use optimization studies and focused on multi-objective trade-offs, as well as model improvements [20,[57][58][59]; however, few studies take the effect of inner-city economic linkages and ecological gradient on land use changes in UA into account [13,60].Cao have proposed land use optimization for urban agglomerations, but he focused on the intermediate process of change between the initial and future patterns and ignored the intrinsic characteristics of urban agglomerations [13].Some scholars have discussed the role of economic linkages for regional collaborative development [31], these studies cannot identify the spatial location preferences of land use.Meanwhile, although objectives, such as soil erosion and climate risk are still important factors for high quality development of UA [49,50], this model argues that using too many objectives would confuse the focus of the study, so only the three objectives with the highest universality are chosen (GDP, ESV, and compactness) and modified to meet the regional collaborative development needs.
In this model, ELI and accessibility were introduced into modified GDP, while NDVI was used in modified ESV; the modification not only give economic and ecological objective a quantitative meaning, but also the ability of spatial location preference for enabling agglomeration, whereas the above objectives had only a quantitative function in the previous study [13,20].Additionally, original NSGA-II used a neighbor-change method to reduce the fragmentation [40,42].In our model, spatial allocation with implicit economic relationships was enhanced through a two-step mutation operator to limit fragmentation.The results show that our model is able to promote agglomeration and optimize ecological objectives around higher value vegetation.Our framework appears to provide more relevant future land use planning solutions that can comprehensively and efficiently simulate UA around different objective priorities.
Although this paper establishes a land use optimization framework applicable to UA, it still has some limitations.Firstly, the raster-based optimization could be more practical if we adopt land use zones or patches as optimization units.Additionally, UA are complex, and the variety of urban development across regions is more pronounced than inland UA-for example, from regional linkages.Computationally, vertical and horizontal gravity flips should be tested, in addition to the rotation used, and more continuous measures of vegetation density could be used.Additionally, risk objectives, such as heat risk and soil erosion risk, may be required in the optimization process, in order to determine the long-term livability of regions.Despite shortcomings, our study provides a good reference and strategy for future planning of UA.

Conclusions
This study established a new multi-objective land use optimization of UA that explored the intensity of economic ties and ecological gradients using the multi-objective NSGA-II algorithm.Three conflicting objectives (modified GDP, modified ESV, and compactness) were proposed.ELI and accessibility factors were proposed as modifying factors for GDP.NDVI were used to modify ESV and limit development within high-value vegetation.Accessibility mutation and gravity flip mutation were constructed for enabling agglomerations.The parallel process was applied in NSGA-II by making eight chromosomes run simultaneously in each iteration, which improved the model efficiency in execution.On this basis, the model was used to optimize the land-use layout in GBA in 2030.A Pareto set was obtained, and its optimization performance was evaluated.
The optimization results showed that the performance of our model was better than original NSGA-II.The value of modified GDP of our model was 7.41% higher than the original NSGA-II for almost the same built-up land increment.Anyway, the comparison showed that our model reduced the conversion rate of high-density vegetation (7.92%) and select relatively low-and medium-density vegetation for built-up land development.The reason for these improvements is that ELI and accessibility caused the clustering of cities, which enhanced modified GDP performance, and the modified ESV preferentially optimized ecological objectives around higher value vegetation.However, prioritizing planning around the modified ESV favored decentralized built-up development.
Our model provided a set of relative optimum land-use arrangements that meet the requirements of the ecological red line, and it has shown high robustness and flexibility in multiple-objective land-use optimization problems.It could simultaneously optimize the land-use spatial layout and improve the economic benefits, without negative effects on

Figure 1 .
Figure 1.Land use map of the GBA in 2015.For modeling purposes, the land-use map is di into 509 × 417 cells, with each cell being 1 × 1 km in size.

Figure 1 .
Figure 1.Land use map of the GBA in 2015.For modeling purposes, the land-use map is divided into 509 × 417 cells, with each cell being 1 × 1 km in size.

Figure 2 .
Figure 2. The distribution of NDVI (a) and DEM (b) in the GBA.

Figure 2 .
Figure 2. The distribution of NDVI (a) and DEM (b) in the GBA.

Figure 2 .
Figure 2. The distribution of NDVI (a) and DEM (b) in the GBA.

Figure 3 .
Figure 3. Workflow procedure of the multi-objective NSGA-II land use optimization framework to derive Pareto-optimal solutions.A linear prediction of past land use was used to provide future land use initialization constraints.The NSGA-II algorithm was then applied using non-dominated sorting, cross-over, accessibility mutation, and gravity flip operators, which are described in the main text.

Figure 4 .
Figure 4. Non-dominated sorting in NSGA-II, where Pi is the parent population while T is the result population after crossover and mutation.Pi+1 represents the result after non-dominated sorting.

Figure 4 .
Figure 4. Non-dominated sorting in NSGA-II, where Pi is the parent population while T is the result population after crossover and mutation.P i+1 represents the result after non-dominated sorting.

Land 2022 , 19 Figure 5 .
Figure 5. Schematic illustration of the procedure for crossover: replacing the center cell in Parent1 if the new center cell type in Parent2 results in neighbors of the same class that outnumbers those in Parent1.

Figure 6 .
Figure 6.Procedure of mutation, where the selected cell in the parent is changed according to the accessibility value.

Figure 5 .
Figure 5. Schematic illustration of the procedure for crossover: replacing the center cell in Parent1 if the new center cell type in Parent2 results in neighbors of the same class that outnumbers those in Parent1.

Figure 5 .
Figure 5. Schematic illustration of the procedure for crossover: replacing the center cell in Parent1 if the new center cell type in Parent2 results in neighbors of the same class that outnumbers those in Parent1.

Figure 6 .
Figure 6.Procedure of mutation, where the selected cell in the parent is changed according to the accessibility value.

Figure 6 .
Figure 6.Procedure of mutation, where the selected cell in the parent is changed according to the accessibility value.

Figure 8 .
Figure 8. Constraints in GBA, in relation to (a) the slope constraint and (b) ecological control area, that cannot be changed.

Figure 8 .
Figure 8. Constraints in GBA, in relation to (a) the slope constraint and (b) ecological control area, cannot be changed.

Figure 9 .
Figure 9.The initial spatial distribution of ESV (a) and GDP (b).

Figure 10 .
Figure 10.The spatial distribution of normalized accessibility (a) and ELI (b).

Figure 9 .
Figure 9.The initial spatial distribution of ESV (a) and GDP (b).

Figure 9 .
Figure 9.The initial spatial distribution of ESV (a) and GDP (b).

Figure 10 .
Figure 10.The spatial distribution of normalized accessibility (a) and ELI (b).

Figure 10 .
Figure 10.The spatial distribution of normalized accessibility (a) and ELI (b).

Figure 11 .
Figure 11.The joint distribution of three objectives (a), and the convergence curve of modified ESV change (b) with iteration.In (a), the objective values are shown at different iteration stages: the blue points represent the initial random population; the green and orange points represent the 20-th and 50-th generation solutions, respectively; the red points are the final result.

Figure 11 .
Figure 11.The joint distribution of three objectives (a), and the convergence curve of modified ESV change (b) with iteration.In (a), the objective values are shown at different iteration stages: the blue points represent the initial random population; the green and orange points represent the 20-th and 50-th generation solutions, respectively; the red points are the final result.

Figure 12 .
Figure 12.Distribution of land use for simulation scenarios: SI (GDP)-optimal GDP objective; SII (ESV)-optimal ESV objective; SIII (compactness)-optimal compactness objective; SIV (NSGA-II)-baseline NSGA-II algorithm, without ELI, accessibility, or NDVI.Left panels show final simulated land use patterns, while right side panels show the corresponding built-up land increment changes.

Figure 12 .
Figure 12.Distribution of land use for simulation scenarios: SI (GDP)-optimal GDP objective; SII (ESV)-optimal ESV objective; SIII (compactness)-optimal compactness objective; SIV (NSGA-II)baseline NSGA-II algorithm, without ELI, accessibility, or NDVI.Left panels show final simulated land use patterns, while right side panels show the corresponding built-up land increment changes.

Table 1 .
The GDP of different cities in the Great Bay Area.

Table 3 .
The unit ESV and GDP in each land use type.

Table 3 .
The unit ESV and GDP in each land use type.

Table 3 .
The unit ESV and GDP in each land use type.

Table 4 .
The objective values, under different scenarios, described in the main text.

Table 4 .
The objective values, under different scenarios, described in the main text.

Table 5 .
The Area change of land use under different scenarios.

Table 6 .
The change of vegetation transfer at different levels in SI and SIV.