Method and System for Defining Elastic Development Boundaries of Urban Agglomerations

Abstract

The present invention discloses a method for delimiting the elastic development boundary of an urban agglomeration, including: S1, establishing a spatial database of evaluation indicators for the importance of ecological protection; S2, performing principal component analysis based on the spatial database of evaluation indicators; S3, simulating the evaluation of the importance of ecological protection in multiple scenarios based on the OWA operator; S4, extracting and selecting the ecological protection priority areas in each scenario; S5, establishing a spatial database of driving factors for land use change; S6, predicting the quantity structure of land use in the target year under the natural development scenario and the ecological protection scenario; S7, simulating the spatial structure of land use in the base year and verifying the simulation accuracy; S8, simulating the spatial structure of land use in the target year under the natural development scenario and the ecological protection scenario, and delimiting the elastic development boundary. The present invention simulates multiple scenarios for evaluating the importance of ecological protection, selects the ecological protection priority areas in the optimal scenario based on the protection efficiency and trade-off degree, simulates the spatial structure of land use in the target year, and delimits the elastic development boundary of the urban agglomeration.

Description

Technical Field

[0001] The present invention relates to the field of land space planning and management technology, and in particular to a method and system for delineating the flexible development boundaries of urban agglomerations based on ecological protection importance evaluation. Background Art

[0002] Urban sprawl is a form of excessive urban spatial growth, which can easily lead to cities deviating from the reasonable development track and causing problems such as large-scale erosion of cultivated land, low land utilization rate and ecological environment damage. The urban flexible development boundary refers to the spatial boundary between the areas where urban development and construction are allowed and prohibited within a certain period of time. The essence of this boundary is to determine the scope of urban development space that balances the contradiction between ecological protection and economic development. While guarding the natural ecological security boundary, it guides the orderly development of urban construction land, thereby effectively governing urbanization-derived problems such as urban sprawl, reserving more ecological space and agricultural space for sustainable development, and ultimately achieving a larger scope and higher level of land space protection. Urban agglomerations are a spatial organizational form that appears when the urbanization process develops to an advanced stage, and are also the main spatial form that carries regional development elements. Scientifically delineating the urban agglomeration's flexible development boundary is of great significance for ensuring regional ecological security and guiding the coordinated development of land space.

[0003] In the context of increasing natural resource constraints, it has become a general consensus to define elastic development boundaries based on anti-planning thinking based on ecological constraints. The key to defining elastic development boundaries lies in two points: one is the perspective of ecological constraints, and the other is the land use change simulation model used in the boundary delineation process. On the one hand, previous studies, whether judging ecological constraints from a single ecological perspective such as ecosystem services or from multiple ecological perspectives such as comprehensive ecological patterns obtained from multiple dimensions, are essentially to clearly assign ecological constraints to various regions through indirect evaluation of the ecological protection level of each region. At present, there is still a lack of technical practices for directly evaluating the ecological protection level of each region, especially from the external functionality and internal stability of the ecosystem, that is, the importance of ecological protection based on ecosystem services and ecological fragility, to determine the ecological constraint intensity of each region.

[0004] On the other hand, many mathematical models, such as the CA-Markov model, the constrained CA model, the neural network model, the GeoSOS model, and the FLUS model, etc., have been widely used in the research on land use change simulation and then delimiting the elastic development boundary. However, the above models are usually only applicable to simulating the spatial changes of urban land, and have weak capabilities in dynamically simulating natural land and synchronously evolving multiple land use patches. The Patch-generating Land Use Simulation (PLUS) model can be used to explore land use transfer rules and predict the patch-level evolution of multiple land use spatial structures. Currently, there is still a lack of corresponding exploration of this model in the research on delimiting the elastic development boundary. In addition, the previous domestic and foreign research on the elastic development boundary mainly focuses on the meso-scale cities, districts, counties, and autonomous prefectures in terms of research scale and research object. There are relatively few studies on delimiting the development boundary at the macro-scale of urban agglomerations. Summary of the Invention

[0005] The purpose of the present invention is to solve the problem that the prior art ignores delimiting the elastic development boundary based on the evaluation of ecological protection importance and the PLUS model at the scale of urban agglomerations, and to provide a method for delimiting the elastic development boundary of urban agglomerations. The steps are scientific, general, and replicable, and are suitable for delimiting the elastic development boundary of urban agglomerations under a specific time frame in any region.

[0006] In order to achieve the above purpose, the technical solutions adopted by the present invention are as follows:

[0007] A method for delimiting the elastic development boundary of urban agglomerations, comprising the following steps:

[0008] S1: Select the evaluation indicators of ecological protection importance, construct an evaluation index system of ecological protection importance with the evaluation of ecological protection importance as the target layer and ecosystem services and ecological vulnerability as the criterion layer, carry out the specific evaluation of the evaluation indicators of ecological protection importance, and establish a spatial database of the evaluation indicators of ecological protection importance;

[0009] S2: Conduct principal component analysis according to the spatial database of the evaluation indicators of ecological protection importance, generate principal component index data and a principal component analysis report, conduct KMO and Bartlett's tests and the cumulative value test of the principal component variance contribution rate, and determine the principal component index data and the criterion weights of each principal component index;

[0010] S3: Set multiple decision risk coefficient scenarios, calculate the rank weights and trade-off degrees of the principal component indicators according to the determined principal component index data, aggregate the rank weights and criterion weights of the principal component indicators under different scenarios by using the OWA operator, and conduct multi-scenario simulation of the evaluation of ecological protection importance;

[0011] S4: According to the simulation results of the ecological protection importance evaluation, extract the ecological protection priority areas under various scenarios, compare the protection efficiency and trade-off degree of the ecological protection priority areas of different scenarios for the ecological protection importance evaluation results of each scenario, and screen out the ecological protection priority areas with the jointly optimal protection efficiency and trade-off degree;

[0012] S5: Select the driving factors of land use change, construct a driving factor system of land use change with natural factors and human factors as the criterion layer, conduct a specific evaluation of the driving factors of land use change, and establish a spatial database of driving factors of land use change;

[0013] S6: According to the historical land use data, obtain the land use quantity structure and land use transfer probability matrix of the target year in the natural development scenario. According to the historical land use data, extract the expansion patch data of each land use type, calculate the land use transfer probability matrix of the ecological protection scenario, and obtain the land use quantity structure of the target year in the ecological protection scenario;

[0014] S7: According to the spatial database of driving factors of land use change and the expansion patch data of each land use type, output the development probability data of each land use type. According to the historical land use data and the development probability data of each land use type, simulate the land use spatial structure data of the base year and conduct a simulation accuracy verification;

[0015] S8: Based on the land use data of the base year and the development probability data of each land use type, simulate the land use spatial structure data of the target year in the natural development scenario; based on the land use data of the base year and the development probability data of each land use type, simulate the land use spatial structure data of the target year in the ecological protection scenario and delimit the elastic development boundary.

[0016] Furthermore, step S1 is specifically as follows: Select the evaluation indicators of ecological protection importance, construct an evaluation index system of ecological protection importance with ecological protection importance evaluation as the target layer and ecosystem services and ecological vulnerability as the criterion layer; combine multi-source data, ArcGIS software and InVEST model to conduct a specific evaluation of the evaluation indicators of ecological protection importance; uniformly normalize the evaluation results to eliminate the influence of dimensions between different evaluation indicators; integrate all normalized evaluation indicators with the help of the ArcGIS platform to establish a spatial database of evaluation indicators of ecological protection importance in the study area;

[0017] Among them, the supply service, regulation service, support service belonging to ecosystem services and the exposure, sensitivity and adaptability belonging to ecological vulnerability are the first-level index layers, and the selected specific indicators are the second-level index layers of the evaluation index system of ecological protection importance. According to the influence on the target layer, the indicators belonging to adaptability among the selected specific indicators are negative indicators, and other indicators are positive indicators;

[0018] The calculation method for normalizing positive indicators is as follows:

[0019]

[0020] The calculation method for normalizing negative indicators is as follows:

[0021]

[0022] In the formula, X norm is the normalized index pixel value, X x is the index pixel value, X max is the maximum value of all pixel values of the index, X min is the minimum value of all pixel values of the index.

[0023] Furthermore, step S2 is specifically as follows: Use the raster sampling tool of QGIS software to sample the raster layers of all evaluation indicators to a point layer converted from a raster with the resolution required by the research, and then import the attribute table data of the point layer into SPSS software for principal component analysis to generate principal component index data and a principal component analysis report; conduct KMO and Bartlett's tests, use the eigenvalue greater than 1 as the criterion for selecting principal components, and conduct a cumulative value test of the principal component variance contribution rate; determine the principal component index data based on the ecological protection importance evaluation indicators, and determine the criterion weights of each principal component index whose sum is 1 according to the principal component variance contribution rate;

[0024] Among them, conducting KMO and Bartlett's tests is specifically as follows:

[0025] If the KMO and Bartlett's test results in the principal component analysis report show that: the KMO sampling adequacy measure is greater than 0.7, and the significance of the Bartlett sphericity test is less than 0.001, then the research data is suitable for principal component analysis; otherwise, the research data is not suitable for principal component analysis;

[0026] Conducting a cumulative value test of the principal component variance contribution rate is specifically as follows:

[0027] If the total variance explanation result in the principal component analysis report shows that: the cumulative value of the principal component variance contribution rate with an eigenvalue greater than 1 exceeds 60%, then the generated principal component index data is relatively reasonable and effective; otherwise, the generated principal component index data is unreasonable and ineffective.

[0028] Furthermore, step S3 is specifically as follows: Set z (z≥5) decision risk coefficient scenarios, and according to the determined principal component index data, use the monotonic rule increasing method to calculate the rank weights of the principal component indices. The specific calculation method is as follows:

[0029]

[0030] In the formula, j is the sequence number; v j is the sequence weight, v j ∈[0,1]; n is the number of principal component indicators; α is the decision-making risk coefficient, α∈(0,∞); w k is the importance level of the principal component indicator; r k is the assignment of the principal component indicator. The indicators are assigned according to the size of the pixel values of the indicators after superimposing the criterion weights. The maximum value is assigned 1, the second largest value is assigned 2, and the minimum value is assigned n;

[0031] The trade-off degree is the compensation degree between the principal component indicators under different decision-making risk coefficients. The specific calculation method is as follows:

[0032]

[0033] Among them, tradeoff is the trade-off degree, and 0≤tradeoff≤1; n is the number of principal component indicators; w k is the importance level of the k-th principal component indicator;

[0034] Using the OWA operator mechanism of the MCE module of TerrSet software, aggregate the sequence weights and criterion weights of the principal component indicators under different scenarios, so as to obtain the spatial evaluation results of the ecological protection importance in z scenarios.

[0035] Furthermore, step S4 is specifically as follows: Based on the simulation results of the ecological protection importance evaluation in z scenarios, calculate the ratio of the cumulative ecological protection importance pixel value to the total ecological protection importance pixel value in descending order, and use the ecological protection importance pixel value corresponding to the ratio of 50% as the ecological protection importance threshold for each scenario;

[0036] In ArcGIS software, use the aggregation tool to aggregate the relatively aggregated or adjacent patches in the area greater than or equal to the ecological protection importance threshold into relatively complete and contiguous patches. The aggregation distance is the required resolution for the study, and the aggregation result is used as the preliminary ecological protection priority area for each scenario;

[0037] To reduce the fragmentation degree of the ecological protection priority areas, the area threshold of the patches to be removed is determined by using the ecological protection priority area patches in scenario z1. Among them, scenario z1 is the scenario where the ordinal weights of the principal component indicators in z scenarios are equal. The ratio of the number of patches with an area smaller than the threshold in the ecological protection priority area of scenario z1 to the total number of all patches is statistically analyzed, as well as the change of its total area with the area threshold. The minimum area threshold corresponding to the inflection point is used as the final patch area threshold, and the corresponding independent and scattered small patches are removed according to the threshold. On this basis, existing various ecological protection areas such as national and provincial nature reserves, forest parks, and scenic spots are merged, and finally the final ecological protection priority areas of each scenario are obtained.

[0038] Compare the protection efficiency and trade-off degree of the final ecological protection priority areas of different scenarios with the evaluation results of the ecological protection importance of scenario z1, so as to screen out the final ecological protection priority areas of the scenario with the jointly optimal protection efficiency and trade-off degree.

[0039] Among them, the calculation method of the protection efficiency is as follows:

[0040]

[0041] In the formula, z is the scenario serial number, z = 1, 2,..., z, P z is the protection efficiency of scenario z, is the average value of the ecological protection importance of scenario z1 within the scope of the ecological protection priority area of scenario z, is the average value of the ecological protection importance of scenario z1 within the entire study area.

[0042] Furthermore, step S5 is specifically as follows: Select the driving factors of land use change, and construct a land use change driving factor system with natural factors and human factors as the criterion layers; Combine multi-source data and ArcGIS software to carry out the specific evaluation of the driving factors of land use change; The obtained evaluation results are uniformly normalized to eliminate the dimensional influence between different driving factors; Integrate all normalized driving factors with the help of the ArcGIS platform to establish a spatial database of the driving factors of land use change in the study area.

[0043] Among them, the terrain and water area belonging to natural factors and the population, economy, location, and infrastructure belonging to human factors are the first-level indicator layers, and the infrastructure includes transportation, medical care, education, and commercial services. The specific indicators selected are the land use change driving factor system of the second-level indicator layer.

[0044] Furthermore, step S6 is specifically as follows: Based on two-phase historical land use data, use the Markov Chain module of the PLUS model to obtain the target-year land use quantity structure and land use transfer probability matrix of the natural development scenario.

[0045] Based on two - period historical land - use data, the land - use expansion module of the PLUS model is used to extract the expansion patch data of each land - use type. Based on this, the proportion of each non - construction land in the expansion patch area of urban construction land is calculated, and the obtained proportion data is used as the negative growth rate of the transfer probability from non - construction land to urban construction land. Thus, the land - use transfer probability matrix under the ecological protection scenario is calculated based on the land - use transfer probability matrix under the natural development scenario;

[0046] Based on the Markov Chain principle and the land - use transfer probability matrix under the ecological protection scenario, the land - use quantity structure in the target year under the ecological protection scenario is calculated;

[0047] Among them, the historical land - use data includes seven land - use types: cultivated land, forest land, grassland, water area, urban construction land, rural residential areas, and unused land. Among them, urban construction land and rural residential areas belong to construction land, and the rest of the land - use types belong to non - construction land.

[0048] Furthermore, step S7 is specifically as follows: Based on the land - use change driving factor spatial database and the expansion patch data of each land - use type, the LEAS module of the PLUS model is used to output the development probability data of each land - use type; calculate and set the parameters required for the CARS module of the PLUS model; use the land - use quantity structure of the second - period historical land - use data as the land - use demand data, and based on the first - period historical land - use data and the development probability data of each land - use type, the CARS module of the PLUS model is used to simulate the land - use spatial structure in the base year; based on the two - period historical land - use data and the base - year land - use spatial structure data, the Confusion Matrix and Fom module of the PLUS model is used to verify the simulation accuracy;

[0049] Among them, calculating and setting the parameters required for the CARS module of the PLUS model is specifically as follows:

[0050] Based on the expansion patch data of each land - use type, calculate the proportion of the expansion patch area of each land - use type in the total expansion area, and use the obtained proportion data as the neighborhood weight parameter; set the transfer - allowed value of each land - use type in the land - transfer matrix to 1, without restricting the transfer between any land - use types; for parameters such as the attenuation coefficient of the decreasing threshold, the probability of random patch seeds, the maximum proportion of random patch seeds, the neighborhood range, and the number of parallel threads, default values are adopted;

[0051] Verifying the simulation accuracy is specifically as follows:

[0052] If the Kappa coefficient is greater than 0.8 and the Fom coefficient is greater than 0.2, the simulation results of the PLUS model are basically consistent with the actual results, and the simulation results are credible; otherwise, the input data and relevant parameters of the PLUS model need to be continuously adjusted until the above Kappa coefficient and Fom coefficient test requirements are met before proceeding with the subsequent steps.

[0053] Further, step S8 is specifically as follows: Using the land use quantity structure of the natural development scenario target year as the land use demand data, and the existing various ecological protection areas as the restricted development areas, based on the base year land use data and the development probability data of each land use type, use the CARS module of the PLUS model to simulate the land use spatial structure data of the natural development scenario target year;

[0054] Using the land use quantity structure of the ecological protection scenario target year as the land use demand data, and the ecological protection priority areas as the restricted development areas, based on the base year land use data and the development probability data of each land use type, use the CARS module of the PLUS model to simulate the land use spatial structure data of the ecological protection scenario target year;

[0055] Designate the urban construction land boundaries in the land use spatial structure data of the natural development scenario and the ecological protection scenario target years as the elastic development boundaries of the research area in the natural development scenario and the ecological protection scenario target years, respectively.

[0056] An urban agglomeration elastic development boundary delineation system includes:

[0057] An evaluation index module for establishing a spatial database of ecological protection importance evaluation indexes;

[0058] A principal component analysis module for performing principal component analysis based on the spatial database of ecological protection importance evaluation indexes;

[0059] An OWA operator module for simulating the ecological protection importance evaluation of multiple scenarios based on the OWA operator;

[0060] A scenario screening module for extracting the ecological protection priority areas of each scenario and screening out the ecological protection priority areas of the optimal scenario;

[0061] A driving factor module for establishing a spatial database of land use change driving factors;

[0062] A demand prediction module for predicting the land use quantity structure of the natural development scenario and the ecological protection scenario target years;

[0063] A simulation and verification module for simulating the land use spatial structure of the base year and verifying the simulation accuracy;

[0064] The boundary demarcation module is used to simulate the land use spatial structure in the target year of the natural development scenario and the ecological protection scenario and demarcate the elastic development boundary.

[0065] Compared with the prior art, the present invention has the following advantages and beneficial effects:

[0066] (1) Through the ecological protection importance evaluation that couples ecosystem service indicators and ecological vulnerability indicators, the present invention directly determines the ecological constraint intensity and ecological protection level of each region, then extracts the ecological protection priority areas based on the ecological protection importance evaluation results, and then uses the restricted development areas obtained from the ecological protection priority areas to conduct land use simulation and further demarcate the development boundary, effectively avoiding the problem that the prior art ignores the demarcation of the elastic development boundary based on the ecological protection importance evaluation. The development boundary demarcated based on the technology of the present invention helps to protect high-quality ecological space, maintain regional ecological security and ensure human well-being. At the same time, an innovative method for screening the optimal ecological protection priority areas from multiple scenarios based on protection efficiency and trade-off degree is proposed, providing valuable technical guidance for the decision-making of natural resource management and ecological protection.

[0067] (2) The present invention uses the PLUS model for land use simulation and further demarcates the urban agglomeration development boundary, providing a mechanism and method for using the PLUS model in the demarcation of the urban agglomeration development boundary, effectively avoiding the problem that the prior art ignores the demarcation of the elastic development boundary of the urban agglomeration based on the PLUS model, helping to make up for the deficiencies of the existing development boundary demarcation technology, and improving the scientificity of the development boundary demarcation. At the same time, the present invention is applicable to the demarcation of the development boundary at the urban agglomeration scale, making up for the deficiencies in the research on the development boundary demarcation technology at the macro scale such as urban agglomerations in the prior art.

[0068] (3) The present invention provides a method for demarcating the development boundary of the natural development scenario and the ecological protection scenario, meeting the requirements for the boundary demarcation methods of different scenarios, facilitating the comparative analysis of the boundary demarcation results of different scenarios, thereby improving the scientificity of the boundary demarcation, and helping to construct a more scientific and reasonable pattern of territorial space development and protection. At the same time, the steps of the present invention are scientific, general and reproducible. The data processing and spatial analysis involved in the process can be realized through software such as ArcGIS widely used in multiple fields such as ecology, geography, and planning at present, and it is applicable to the demarcation of the urban agglomeration development boundary in the study area under a specific time frame in any region. Description of the Drawings

[0069] Figure 1 It is a flow schematic diagram of the method for demarcating the elastic development boundary of the urban agglomeration in the embodiment.

[0070] Figure 2 It is a statistical analysis diagram of the ratio of the number of patches with an area less than the threshold to all patches and its total area in the ecological protection priority area in Scenario 5 in the embodiment changing with the area threshold. Specific implementation manners

[0071] The following further explains the method and system for delimiting the elastic development boundary of the urban agglomeration of the present invention in conjunction with the accompanying drawings and specific embodiments.

[0072] Please refer to Figure 1 , the present invention discloses a method for delimiting the elastic development boundary of an urban agglomeration, including the following steps:

[0073] S1: Select evaluation indicators of the importance of ecological protection, construct an evaluation index system of the importance of ecological protection with the evaluation of the importance of ecological protection as the target layer and ecosystem services and ecological vulnerability as the criterion layer, conduct specific evaluations of the evaluation indicators of the importance of ecological protection, and establish a spatial database of evaluation indicators of the importance of ecological protection.

[0074] S2: Conduct principal component analysis based on the spatial database of evaluation indicators of the importance of ecological protection, generate principal component index data and a principal component analysis report, conduct KMO and Bartlett's tests and tests on the cumulative value of the principal component variance contribution rate, and determine the principal component index data and the criterion weights of each principal component index.

[0075] S3: Set multiple decision risk coefficient scenarios, calculate the rank weights and trade-off degrees of the principal component indicators according to the determined principal component index data, aggregate the rank weights and criterion weights of the principal component indicators in different scenarios by using the OWA operator, and conduct multi-scenario simulation evaluations of the importance of ecological protection.

[0076] S4: According to the simulation results of the evaluation of the importance of ecological protection, extract the ecological protection priority areas in various scenarios, compare the protection efficiency and trade-off degrees of the ecological protection priority areas in different scenarios for the evaluation results of the importance of ecological protection in each scenario, and screen out the ecological protection priority areas with the jointly optimal protection efficiency and trade-off degree.

[0077] S5: Select the driving factors of land use change, construct a driving factor system of land use change with natural factors and human factors as the criterion layer, conduct specific evaluations of the driving factors of land use change, and establish a spatial database of driving factors of land use change.

[0078] S6: According to the historical land use data, obtain the target-year land use quantity structure and land use transfer probability matrix of the natural development scenario, extract the expansion patch data of each land use type according to the historical land use data, calculate the land use transfer probability matrix of the ecological protection scenario, and obtain the target-year land use quantity structure of the ecological protection scenario.

[0079] S7: Based on the spatial database of land use change driving factors and the expansion patch data of each land use type, output the development probability data of each land use type. According to the historical land use data and the development probability data of each land use type, simulate the land use spatial structure data in the base year and conduct simulation accuracy verification.

[0080] S8: Based on the land use data in the base year and the development probability data of each land use type, simulate the land use spatial structure data in the target year of the natural development scenario; based on the land use data in the base year and the development probability data of each land use type, simulate the land use spatial structure data in the target year of the ecological protection scenario and delimit the elastic development boundary.

[0081] In step S1, select the evaluation indicators of ecological protection importance, and construct an evaluation index system of ecological protection importance with the evaluation of ecological protection importance as the target layer, and ecosystem services and ecological vulnerability as the criterion layer. The supply service, regulation service, support service belonging to ecosystem services, and exposure, sensitivity and adaptability belonging to ecological vulnerability are the first-level index layers, and the selected specific indicators are the second-level index layers of the evaluation index system of ecological protection importance. According to the influence on the target layer, the indicators belonging to adaptability among the selected specific indicators are divided into negative indicators, and the rest of the indicators are divided into positive indicators. Then, combined with multi-source data, ArcGIS software and InVEST model, conduct the specific evaluation of the evaluation indicators of ecological protection importance. The obtained evaluation results are uniformly normalized to eliminate the influence of dimensions between different evaluation indicators. Integrate all normalized evaluation indicators with the help of the ArcGIS platform to establish a spatial database of evaluation indicators of ecological protection importance in the study area.

[0082] For positive indicators and negative indicators, the range method normalization is carried out using formula (1) and formula (2) respectively. Among them, the calculation method of range method normalization for positive indicators is:

[0083]

[0084] The calculation method of range method normalization for negative indicators is:

[0085]

[0086] In the formula, X norm is the normalized index pixel value, X x is the index pixel value, X max is the maximum value of all index pixels, X min is the minimum value of all index pixels.

[0087] In step S2, based on the spatial database of ecological protection importance evaluation indicators, the raster sampling tool of QGIS software is used to sample the raster layers of all evaluation indicators into a point layer converted from a raster with the resolution required by the research. Then, the attribute table data of the point layer is imported into SPSS software for principal component analysis to generate principal component index data and a principal component analysis report. Then, KMO and Bartlett's tests are performed, and the criterion for selecting principal components is that the eigenvalue is greater than 1, and the cumulative value of the principal component variance contribution rate is tested. Finally, the principal component index data based on the ecological protection importance evaluation indicators is determined, and the criterion weights of each principal component index with a sum of 1 are determined according to the principal component variance contribution rate.

[0088] If the results of the KMO and Bartlett's tests in the principal component analysis report show that the KMO sampling adequacy measure is greater than 0.7 and the significance of the Bartlett sphericity test is less than 0.001, then the research data is suitable for principal component analysis; otherwise, the research data is not suitable for principal component analysis.

[0089] If the total variance explained result in the principal component analysis report shows that the cumulative value of the principal component variance contribution rate with an eigenvalue greater than 1 exceeds 60%, then the generated principal component index data is relatively reasonable and effective; otherwise, the generated principal component index data is unreasonable and ineffective.

[0090] In step S3, z (z≥5) decision risk coefficient scenarios are set. Based on the principal component index data, the rank weights of the principal component indexes are calculated using the monotonically increasing rule method. The specific calculation methods are shown in formulas (3) and (4):

[0091]

[0092] In the formula, j is the rank; v j is the rank weight, v j ∈[0,1]; n is the number of principal component indexes; α is the decision risk coefficient, α∈(0,∞); w k is the importance level of the principal component index; r k is the assignment of the principal component index. The indexes are assigned according to the size of the pixel value of the index after superimposing the criterion weight. The maximum value is assigned 1, the second largest value is assigned 2, and the minimum value is assigned n.

[0093] The trade-off represents the compensation degree between principal component indexes under different decision risk coefficients. The specific calculation method is shown in formula (5):

[0094]

[0095] Among them, tradeoff is the trade-off, and 0≤tradeoff≤1; n is the number of principal component indexes; wk is the importance level of the k-th principal component index.

[0096] Using the OWA operator mechanism of the MCE module in TerrSet software, aggregate the rank weights and criterion weights of the principal component indices under different scenarios, so as to obtain the spatial evaluation results of ecological protection importance under z scenarios.

[0097] In step S4, based on the simulation results of ecological protection importance evaluation under z scenarios, calculate the ratio of the cumulative ecological protection importance pixel value to the total ecological protection importance pixel value in descending order, and take the ecological protection importance pixel value corresponding to the ratio of 50% as the ecological protection importance threshold for each scenario. Then, use the aggregation tool in ArcGIS software to aggregate the relatively aggregated or adjacent patches in the area greater than or equal to the ecological protection importance threshold into relatively complete and contiguous patches. The aggregation distance is the required resolution of the study, and the aggregation result is the preliminary ecological protection priority area for each scenario.

[0098] To reduce the fragmentation degree of the ecological protection priority area, use the patch of the ecological protection priority area in scenario z1 to determine the threshold of the patch area to be removed. Among them, scenario z1 is the scenario where the rank weights of each index are equal among z scenarios. Taking 2 km 2 as the initial value and 2 km 2 as the step length, statistically analyze the ratio of the number of patches with an area less than the threshold to the total number of patches in the ecological protection priority area of scenario z1, and the change of its total area with the area threshold. Take the minimum area threshold corresponding to the inflection point (the turning point where the growth rate first becomes 0) as the final patch area threshold, and remove the corresponding independent and scattered small patches according to the threshold. On this basis, merge existing various ecological protection areas such as national and provincial nature reserves, forest parks, and scenic spots, and finally obtain the final ecological protection priority area for each scenario. Compare the protection efficiency and trade-off degree of the final ecological protection priority areas of different scenarios with the ecological protection importance evaluation results of scenario z1, so as to screen out the final ecological protection priority area of the scenario with the jointly optimal protection efficiency and trade-off degree (the higher the protection efficiency, the lower the trade-off degree, and the better the corresponding scenario);

[0099] The specific calculation method of the protection efficiency is shown in formula (6):

[0100]

[0101] In the formula, z is the scenario number, z = 1, 2,... z, P z is the protection efficiency of scenario z, is the average value of the ecological protection importance of scenario z1 within the scope of the ecological protection priority area of scenario z, is the average value of the ecological protection importance of scenario z1 within the entire study area.

[0102] In step S5, the driving factors of land use change are selected to construct a driving factor system of land use change with natural factors and human factors as the criterion layer. The terrain and water area belonging to natural factors, and population, economy, location, and infrastructure (including transportation, medical care, education, and commercial services) belonging to human factors are the first-level index layer, and the specific selected indicators are the second-level index layer of the driving factor system of land use change. Then, combined with multi-source data and ArcGIS software, a specific evaluation of the driving factors of land use change is carried out. The obtained evaluation results are uniformly normalized using formula (1) to eliminate the dimensionality influence between different driving factors. With the help of the ArcGIS platform, all normalized driving factors are integrated to establish a spatial database of the driving factors of land use change in the study area.

[0103] In step S6, based on two-phase historical land use data (including seven land use types: cultivated land, forest land, grassland, water area, urban construction land, rural residential areas, and unused land. Among them, urban construction land and rural residential areas belong to construction land, and the remaining land use types belong to non-construction land), the Markov Chain module of the PLUS model is used to obtain the target-year land use quantity structure and land use transfer probability matrix under the natural development scenario.

[0104] Then, based on the two-phase historical land use data, the land use expansion module of the PLUS model is used to extract the expansion patch data of each land use type. Based on this, the proportion of each non-construction land in the urban construction land expansion patch area is calculated, and the obtained proportion data is used as the negative growth rate of the transfer probability from non-construction land to urban construction land. Thus, based on the land use transfer probability matrix under the natural development scenario, the land use transfer probability matrix under the ecological protection scenario is calculated. Finally, based on the Markov Chain principle and the land use transfer probability matrix under the ecological protection scenario, the target-year land use quantity structure under the ecological protection scenario is calculated.

[0105] In step S7, based on the spatial database of the driving factors of land use change and the expansion patch data of each land use type, the LEAS module of the PLUS model is used to output the development probability data of each land use type. Then, the parameters required for the CARS module of the PLUS model are calculated and set. Next, taking the land use quantity structure of the second-phase historical land use data as the land use demand data, based on the first-phase historical land use data and the development probability data of each land use type, the CARS module of the PLUS model is used to simulate the land use spatial structure in the year to which the second-phase historical land use data belongs, that is, the base year. Finally, based on the two-phase historical land use data and the simulated land use spatial structure data of the base year, the simulation accuracy is verified using the Confusion Matrix and Fom module of the PLUS model.

[0106] Based on the expansion patch data of each land use type, calculate the proportion of the expansion patch area of each land use type in the total expansion area, and use the obtained proportion data as the neighborhood weight parameter. Then, set the allowable transfer value of each land use type in the land transfer matrix to 1, without restricting the transfer between any land use types. Parameters such as the attenuation coefficient of the decreasing threshold, the probability of random patch seeds, the maximum proportion of random patch seeds, the neighborhood range, and the number of parallel threads adopt default values.

[0107] If the Kappa coefficient is greater than 0.8 and the Fom coefficient is greater than 0.2, the simulation results of the PLUS model are basically consistent with the real results, and the simulation results are credible; otherwise, it is necessary to continuously adjust the input data and related parameters of the PLUS model until the above Kappa coefficient and Fom coefficient test requirements are met before proceeding with the subsequent steps.

[0108] In step S8, using the land use quantity structure of the natural development scenario target year as the land use demand data, and the existing various ecological protection areas as the restricted development areas, based on the land use data of the base year and the development probability data of each land use type, use the CARS module of the PLUS model to simulate the land use spatial structure data of the natural development scenario target year.

[0109] Then, using the land use quantity structure of the ecological protection scenario target year as the land use demand data, and the ecological protection priority areas as the restricted development areas, based on the land use data of the base year and the development probability data of each land use type, use the CARS module of the PLUS model to simulate the land use spatial structure data of the ecological protection scenario target year.

[0110] Finally, demarcate the urban construction land boundaries in the land use spatial structure data of the natural development scenario and the ecological protection scenario target years as the elastic development boundaries of the research area in the natural development scenario and the ecological protection scenario target years respectively.

[0111] As can be seen from the above, a method for delimiting the elastic development boundary of an urban agglomeration provided by the present invention first selects indicators to construct an evaluation index system for the importance of ecological protection, and conducts a specific evaluation of the evaluation indicators for the importance of ecological protection to establish a spatial database of the evaluation indicators for the importance of ecological protection. Then, principal component analysis is carried out based on the spatial database of the evaluation indicators, and KMO and Bartlett's tests as well as the cumulative value test of the principal component variance contribution rate are carried out. On this basis, the principal component index data and its criterion weights are determined. Set z (z≥5) decision risk coefficient scenarios, calculate the rank weights and trade-off degrees of the principal component indicators, and on this basis, conduct a multi-scenario simulation of the importance evaluation of ecological protection based on the OWA operator. Extract the ecological protection priority areas under various scenarios according to the simulation results of the importance evaluation of ecological protection, and calculate the protection efficiency of the ecological protection priority areas in each scenario, so as to screen out the ecological protection priority areas with the jointly optimal protection efficiency and trade-off degree.

[0112] Then, select indicators to construct a system of driving factors for land use change, and establish a spatial database of the driving factors for land use change. Then, use the PLUS model to predict the quantity structure of land use in the target year under the natural development scenario and the ecological protection scenario, and use the PLUS model to simulate the spatial structure of land use in the base year, and conduct a simulation accuracy verification. Finally, use the PLUS model to simulate the spatial structure of land use in the target year under the natural development scenario and the ecological protection scenario, and use the urban construction land boundaries in the two scenarios as the elastic development boundaries of the study area in the target year. The method steps for delimiting the elastic development boundary of the urban agglomeration of the present invention are scientific, general and replicable, and are suitable for being applied to the delimitation of the elastic development boundary of the urban agglomeration under a specific time frame in any region, providing a theoretical basis and technical support for the development boundary decision-making and spatial control of the territorial space planning.

[0113] Next, select the Guangdong-Hong Kong-Macao Greater Bay Area as the research area to elaborate in detail the method for delimiting the elastic development boundary of the urban agglomeration of the present invention.

[0114] Step 1: Establish a spatial database of the evaluation indicators for the importance of ecological protection.

[0115] Based on the ecological environment characteristics of the research area, as well as relevant literature induction and policy interpretation, 17 indicators are selected to construct an evaluation index system for the importance of ecological protection with the importance evaluation of ecological protection as the target layer and ecosystem services and ecological vulnerability as the criterion layer. The supply service, regulation service, and support service belonging to ecosystem services, as well as the exposure, sensitivity, and adaptability belonging to ecological vulnerability are the first-level index layers, and the specific indicators selected are the second-level index layers of the evaluation index system for the importance of ecological protection. According to the influence on the target layer, the indicators belonging to adaptability among the selected specific indicators are classified as negative indicators, and the rest of the indicators are classified as positive indicators. The evaluation index system for the importance of ecological protection is shown in Table 1.

[0116] Table 1 Evaluation Index System for the Importance of Ecological Protection

[0117]

[0118]

[0119] Then, combined with multi-source data (the specific data types and sources are shown in Table 2), ArcGIS software, and the InVEST model, the specific evaluation of the ecological protection importance evaluation indicators is carried out.

[0120] Table 2 Data Description

[0121]

[0122]

[0123] Among them, the specific evaluation methods for each evaluation indicator are shown in Table 3.

[0124] Table 3 Explanation of the Evaluation Methods for the Ecological Protection Importance Evaluation Indicators in the Guangdong-Hong Kong-Macao Greater Bay Area

[0125]

[0126]

[0127] All the above data are resampled or spatially interpolated to a 100m spatial resolution through ArcGIS software.

[0128] The obtained evaluation results are uniformly normalized to eliminate the dimensional influence between different evaluation indicators. With the help of the ArcGIS platform, all the normalized evaluation indicators are integrated to establish a spatial database of the ecological protection importance evaluation indicators in the study area.

[0129] The positive and negative indicators are normalized by the range method using Formula (1) and Formula (2) respectively. Among them, the calculation method for normalizing the positive indicator by the range method is:

[0130]

[0131] The calculation method for normalizing the negative indicator by the range method is:

[0132]

[0133] Among them, X norm is the pixel value of the normalized indicator, X x is the pixel value of the indicator, X max is the maximum value of all the pixel values of the indicator, and X min is the minimum value of all the pixel values of the indicator.

[0134] Step 2: Conduct principal component analysis based on the evaluation index spatial database.

[0135] Based on the ecological protection importance evaluation index spatial database, using the raster sampling tool of QGIS software, sample all the evaluation index raster layers to the point layer converted from the 100m resolution raster, and then import the attribute table data of the point layer into SPSS software for principal component analysis to generate principal component index data and a principal component analysis report. Then conduct KMO and Bartlett's tests, and then use the criterion that the eigenvalue is greater than 1 as the standard for selecting principal components to conduct the cumulative value test of the principal component variance contribution rate. Finally, 5 principal component index data based on the ecological protection importance evaluation index are determined, and the criterion weights of each principal component index with a sum of 1 are determined according to the principal component variance contribution rate, as shown in Table 4.

[0136] Table 4 Eigenvalues, variance contribution rates, and criterion weights of each principal component index

[0137]

[0138]

[0139] The principal component analysis report of this embodiment shows that the KMO sampling adequacy measure is 0.813 (>0.7), and the significance of the Bartlett sphericity test is 0.000 (<0.001). Therefore, the data of this embodiment passes the KMO test and Bartlett's test and is suitable for principal component extraction. At the same time, the cumulative value of the variance contribution rates of the 5 principal component indicators with eigenvalues greater than 1 exceeds 60%. Therefore, the data of this embodiment passes the cumulative value test of the principal component variance contribution rate, and the generated principal component index data is relatively reasonable and effective.

[0140] Step 3: Conduct multi-scenario ecological protection importance evaluation simulations based on the OWA operator.

[0141] Set 9 decision risk coefficient scenarios (α = 0.00001, 0.1, 0.2, 0.5, 1, 2, 5, 10, 100000), and then use the monotonic rule increasing method to calculate the rank weights based on the 5 principal component indicators. The specific calculation methods are shown in Formulas 3 and 4:

[0142]

[0143] In the formula, j is the rank; v j is the rank weight, v j ∈[0,1]; n is the number of principal component indicators; α is the decision risk coefficient, α ∈ (0, ∞); w k is the importance level of the principal component indicator; r kThe assignment of the principal component indicators is to assign values to the indicators according to the size after weighting by the pixel value superposition criterion of the indicators. The maximum value is assigned 1, the second largest value is assigned 2, and the minimum value is assigned n.

[0144] The trade-off represents the compensation degree between the principal component indicators under different decision risk coefficients. The specific calculation method is shown in Formula 5:

[0145]

[0146] where tradeoff is the trade-off, and 0 ≤ tradeoff ≤ 1; n is the number of principal component indicators; w k is the importance level of the k-th principal component indicator. The decision risk coefficients, trade-offs, and ordinal weights in 9 scenarios are shown in Table 5.

[0147] Table 5 Decision risk coefficients, trade-offs, and ordinal weights in different scenarios

[0148]

[0149]

[0150] Using the OWA operator mechanism of the MCE module in TerrSet software, aggregate the ordinal weights and criterion weights of the principal component indicators in different scenarios, so as to obtain the spatial evaluation distribution map of ecological protection importance in 9 scenarios.

[0151] Step 4: Extract and select the ecological protection priority areas in each scenario.

[0152] Based on the simulation results of the ecological protection importance evaluation in 9 scenarios, calculate the ratio of the cumulative ecological protection importance pixel value to the total ecological protection importance pixel value in descending order. Take the ecological protection importance pixel value corresponding to the ratio of 50% as the ecological protection importance threshold in each scenario. Then, use the aggregation tool in ArcGIS software to aggregate the relatively aggregated or adjacent patches in the area greater than or equal to the ecological protection importance threshold into relatively complete contiguous patches. The aggregation distance is 100m, and the aggregation result is used as the preliminary ecological protection priority area in each scenario.

[0153] To reduce the fragmentation degree of the ecological protection priority area, use the patch map of the ecological protection priority area in Scenario 5 to determine the threshold of the patch area to be removed. Among them, Scenario 5 is the scenario where the ordinal weights of each indicator are equal in 9 scenarios. Taking 2km 2 as the initial value and 2km 2 as the step size, statistically analyze the ratio of the number of patches with an area less than the threshold to all patches in the ecological protection priority area in Scenario 5, and the change of its total area with the area threshold. See Figure 2 . This embodiment willFigure 2 The minimum area threshold corresponding to the inflection point (the first turning point where the growth rate is 0) is 26 km 2 Taking the inflection point as the final patch area threshold, and removing the corresponding independent and scattered small patches according to the threshold. On this basis, existing various ecological protection areas such as national and provincial nature reserves, forest parks, and scenic spots are merged, and finally the spatial distribution map of the final ecological protection priority areas in each scenario is obtained. Compare the protection efficiency and trade-off degree of the final ecological protection priority areas in different scenarios for the ecological protection importance evaluation results of Scenario 5, so as to screen out the final ecological protection priority areas of the scenario with the jointly optimal protection efficiency and trade-off degree (the higher the protection efficiency and the lower the trade-off degree, the better the corresponding scenario). The specific calculation method of the protection efficiency is shown in formula (6):

[0154]

[0155] In the formula, z is the scenario number, z = 1, 2,..., z, P z is the protection efficiency of scenario z, is the average value of the ecological protection importance of scenario z1 within the ecological protection priority area of scenario z, is the average value of the ecological protection importance of scenario z1 within the entire study area, which is approximately 0.4069. The average value of the ecological protection importance of Scenario 5 and the protection efficiency within the ecological protection priority areas of each scenario are shown in Table 6.

[0156] Table 6 Average value of the ecological protection importance of Scenario 5 and protection efficiency within the ecological protection priority areas of different scenarios

[0157]

[0158] Based on the above analysis, from the perspective of protection efficiency, the protection efficiencies of Scenario 4 and Scenario 5 are the highest, both being 1.0553. Considering the trade-off degree again, the trade-off degree of Scenario 4 is lower than that of Scenario 5, that is, Scenario 4 compensates for the trade-off impact between the principal component indicators to a certain extent compared with Scenario 5. Therefore, in this embodiment, the ecological protection priority area of Scenario 4 is finally selected as the ecological protection priority area of the optimal scenario.

[0159] Step Five: Establish a spatial database of land use change driving factors.

[0160] Based on the natural background conditions and social and economic conditions of the study area, 19 driving factors are selected to construct a land use change driving factor system with natural factors and human factors as the criterion layers. The terrain and water area belonging to natural factors and population, economy, location, and infrastructure (including transportation, medical care, education, and commercial services) belonging to human factors are the first-level index layers, and the specific selected indicators are the second-level index layers of the land use change driving factor system. The land use change driving factor system in the Guangdong-Hong Kong-Macao Greater Bay Area is shown in Table 7.

[0161] Table 7 Driving factor system for land use change in the Guangdong-Hong Kong-Macao Greater Bay Area

[0162]

[0163]

[0164] Then, combined with multi-source data (see Table 2 for the types and sources of relevant data) and ArcGIS software, the specific evaluation of the driving factors of land use change is carried out. The specific evaluation methods for each driving factor are shown in Table 8.

[0165] Table 8 Explanation of the evaluation methods for the driving factors of land use change in the Guangdong-Hong Kong-Macao Greater Bay Area

[0166]

[0167]

[0168] The evaluation results are uniformly normalized using a formula to eliminate the dimensionality effects between different driving factors. With the help of the ArcGIS platform, all normalized driving factors are integrated to establish a spatial database of the driving factors of land use change in the study area.

[0169] Step 6: Predict the quantity structure of land use in the target year under the natural development scenario and the ecological protection scenario.

[0170] Based on the historical land use data in 2010 and 2020 (including seven land use types: cultivated land, forest land, grassland, water area, urban construction land, rural settlements, and unused land. Among them, urban construction land and rural settlements belong to construction land, and the rest of the land use types belong to non-construction land), the Markov Chain module of the PLUS model is used to obtain the quantity structure of land use in 2050 under the natural development scenario (as shown in Table 12) and the land use transfer probability matrix (as shown in Table 9).

[0171] Table 9 Land use transfer probability matrix under the natural development scenario

[0172]

[0173] Next, based on the historical land use spatial structure data in 2010 and 2020, the expansion patch data of each land use type was extracted using the land use expansion module of the PLUS model. Based on this, the proportion of each non-construction land in the expansion patch area of urban construction land was calculated, as shown in Table 10. And the obtained proportion data was used as the negative growth rate of the transfer probability from non-construction land to urban construction land, that is, the transfer probabilities of cultivated land, forest land, grassland, water area, and unused land to urban construction land decreased by 45.64%, 27.87%, 2.21%, 20.19%, and 0.19% respectively.

[0174] Table 10 Proportion of each non-construction land type in the expansion area of urban construction land from 2010 to 2020

[0175]

[0176]

[0177] Based on the land use transfer probability matrix under the natural development scenario, the land use transfer probability matrix under the ecological protection scenario was calculated, as shown in Table 11.

[0178] Table 11 Land use transfer probability matrix under the ecological protection scenario

[0179]

[0180] Finally, based on the Markov Chain principle and the land use transfer probability matrix under the ecological protection scenario, the land use quantity structure in 2050 under the ecological protection scenario was calculated. As shown in Table 12.

[0181] Table 12 Land use quantity structure in 2050 predicted based on Markov Chain

[0182]

[0183] Step 7: Simulate the land use spatial structure in the base year and verify the simulation accuracy.

[0184] Based on the land use change driving factor spatial database and the expansion patch data of each land use type, the development probability spatial distribution map of each land use type was output using the LEAS module of the PLUS model.

[0185] Then, calculate and set the parameters required for the CARS module of the PLUS model. Calculate the proportion of the expansion patch area of each land use type in the total expansion area based on the expansion patch data of each land use type, and use the obtained proportion data as the neighborhood weight parameter (as shown in Table 13). Then, set the transfer permission value of each land use type in the land transfer matrix to 1, without restricting the transfer between any land use types. Parameters such as the decay coefficient (0.9) of the decreasing threshold, the probability (0.1) of the random patch seed, the maximum proportion (0.001) of the random patch seed, the neighborhood range (3), and the number of parallel threads (1) adopt the default values.

[0186] Table 13 Neighborhood weight parameters of each land use type

[0187]

[0188] Next, use the historical land use quantity structure in 2020 as the land use demand data. Based on the historical land use data in 2010 and the development probability data of each land use type, use the CARS module of the PLUS model to simulate the land use spatial structure data in 2020. Finally, based on the historical land use data in 2010 and 2020 and the simulated land use spatial structure data in 2020, use the Confusion Matrix and Fom module of the PLUS model to verify the simulation accuracy. The output result of the Confusion Matrix and Fom module shows that the Kappa coefficient is 0.9389 (>0.8). Therefore, the simulation result of the PLUS model is basically consistent with the real result, and the simulation result is credible.

[0189] Step 8: Simulate the land use spatial structure of the target year in the natural development scenario and the ecological protection scenario and delimit the elastic development boundary.

[0190] Taking the land use quantity structure in the natural development scenario in 2050 as the land use demand data, and the existing various ecological protection areas as the restricted development areas, based on the land use spatial structure data in 2020 and the development probability data of each land use type, the spatial distribution map of the land use spatial structure in 2050 in the natural development scenario is simulated by using the CARS module of the PLUS model. Then, taking the land use quantity structure in the ecological protection scenario in 2050 as the land use demand data, and the ecological protection priority areas as the restricted development areas, based on the land use spatial structure data in 2020 and the development probability data of each land use type, the spatial distribution map of the land use spatial structure in 2050 in the ecological protection scenario is simulated by using the CARS module of the PLUS model. Finally, the urban construction land boundaries in the spatial distribution maps of the land use spatial structures in 2050 in the natural development scenario and the ecological protection scenario are respectively designated as the elastic development boundaries in 2050 in the natural development scenario and the ecological protection scenario of the study area.

[0191] In summary, the method for delimiting the elastic development boundary of the urban agglomeration of the present invention has the following advantages and beneficial effects:

[0192] (1) Through the ecological protection importance evaluation that couples the ecosystem service index and the ecological vulnerability index, the present invention directly determines the ecological constraint intensity and the ecological protection level of each region, then extracts the ecological protection priority areas based on the ecological protection importance evaluation results, and then uses the restricted development areas obtained from the ecological protection priority areas to conduct land use simulation and further delimit the development boundary, effectively avoiding the problem that the prior art ignores the delimitation of the elastic development boundary based on the ecological protection importance evaluation. The development boundary delimited based on the technology of the present invention helps to protect high-quality ecological space, maintain regional ecological security and ensure human well-being. At the same time, a method for screening the ecological protection priority areas of the optimal scenario from multiple scenarios based on the protection efficiency and trade-off degree is innovatively proposed, providing valuable technical guidance for the decision-making of natural resource management and ecological protection.

[0193] (2) The present invention uses the PLUS model to conduct land use simulation and further delimit the development boundary of the urban agglomeration, providing a mechanism and method for using the PLUS model to delimit the development boundary of the urban agglomeration, effectively avoiding the problem that the prior art ignores the delimitation of the elastic development boundary of the urban agglomeration based on the PLUS model, helping to make up for the deficiencies of the existing development boundary delimitation technology, and improving the scientificity of the development boundary delimitation. At the same time, the present invention is applicable to the delimitation of the development boundary at the urban agglomeration scale, making up for the deficiencies in the research on the development boundary delimitation technology at the macro scale such as urban agglomerations in the prior art.

[0194] (3) The present invention provides a method for delimiting the development boundary in natural development scenarios and ecological protection scenarios, meeting the requirements for delimiting methods in different scenarios, facilitating the comparative analysis of the delimitation results in different scenarios, thereby improving the scientificity of boundary delimitation and contributing to the construction of a more scientific and reasonable pattern for the development and protection of territorial space. Meanwhile, the steps of the present invention are scientific, general, and replicable. The data processing and spatial analysis involved in the process can be achieved through software such as ArcGIS widely used in multiple fields such as ecology, geography, and planning at present, and it is applicable to the delimitation of the development boundary of the urban agglomeration in the study area under a specific time frame in any region.

[0195] The present invention also discloses a system for delimiting the flexible development boundary of an urban agglomeration, including: an evaluation index module for establishing a spatial database of ecological protection importance evaluation indexes; a principal component analysis module for performing principal component analysis based on the spatial database of ecological protection importance evaluation indexes; an OWA operator module for simulating the ecological protection importance evaluation in multiple scenarios based on the OWA operator; a scenario screening module for extracting the ecological protection priority areas in each scenario and screening out the ecological protection priority areas of the optimal scenario; a driving factor module for establishing a spatial database of land use change driving factors; a demand prediction module for predicting the quantity structure of land use in the target year of natural development scenarios and ecological protection scenarios; a simulation and verification module for simulating the land use spatial structure in the base year and verifying the simulation accuracy; and a boundary delimitation module for simulating the land use spatial structure in the target year of natural development scenarios and ecological protection scenarios and delimiting the flexible development boundary. The system for delimiting the flexible development boundary of an urban agglomeration according to the present invention can execute the method for delimiting the flexible development boundary of an urban agglomeration of the present invention, can execute any combination of the implementation steps of the method embodiments, and has the corresponding functions and beneficial effects of the method.

[0196] Although the present invention has been described in the context of functional modules, it should be understood that, unless otherwise stated, one or more of the functions and / or features can be integrated in a single physical device and / or software module, or one or more functions and / or features can be implemented in separate physical devices or software modules. It can also be understood that a detailed discussion of the actual implementation of each module is not necessary for understanding the present invention. More precisely, considering the attributes, functions, and internal relationships of the various functional modules in the system disclosed herein, the actual implementation of the module will be understood within the routine skills of an engineer. Therefore, those skilled in the art can implement the present invention as set forth in the claims without undue experimentation using ordinary skills. It can also be understood that the specific concepts disclosed are merely illustrative and are not intended to limit the scope of the present invention, which is determined by the full scope of the appended claims and their equivalents.

[0197] If a function is implemented in the form of a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, in essence, or the part that contributes to the prior art, or a part of this technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions for causing a computer device (which may be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the methods of various embodiments of the present invention. The foregoing storage medium includes: various media such as USB flash drives, mobile hard disks, read-only memories (ROMs), random access memories (RAMs), magnetic disks, or optical discs that can store program codes.

[0198] More specific examples (nonexhaustive list) of computer-readable media include the following: electrical connection parts (electronic devices) having one or more wirings, portable computer disk cartridges (magnetic devices), random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), optical fiber devices, and portable compact disc read-only memory (CDROM). Additionally, the computer-readable medium can even be paper or other suitable media on which a program can be printed, because the program can be obtained electronically, for example, by optically scanning the paper or other media, then editing, interpreting, or processing it in other suitable ways as necessary, and then storing it in a computer memory.

[0199] Each part of the present invention can be implemented by hardware, software, firmware, or a combination thereof. In the above-described embodiments, multiple steps or methods can be implemented by software or firmware stored in a memory and executed by a suitable instruction execution system. For example, if implemented by hardware, as in another embodiment, any one or a combination of the following techniques well known in the art can be used: discrete logic circuits having logic gate circuits for implementing logical functions on data signals, application-specific integrated circuits having suitable combinational logic gate circuits, programmable gate arrays (PGAs), field-programmable gate arrays (FPGAs), etc.

[0200] The above description is a detailed description of the preferred and feasible embodiments of the present invention, but the embodiments are not intended to limit the scope of the patent application of the present invention. Any equivalent changes or modifications made under the technical spirit disclosed by the present invention shall fall within the scope of the patent covered by the present invention.

Claims

1. A method for delimiting the elastic development boundary of an urban agglomeration, characterized in that, It includes the following steps: S1: Select the evaluation indicators of ecological protection importance, construct an evaluation index system of ecological protection importance with the evaluation of ecological protection importance as the target layer, and ecosystem services and ecological vulnerability as the criterion layers, conduct the specific evaluation of the evaluation indicators of ecological protection importance, and establish a spatial database of evaluation indicators of ecological protection importance; S2: Conduct principal component analysis based on the spatial database of evaluation indicators of ecological protection importance, generate principal component index data and a principal component analysis report, conduct KMO and Bartlett's tests and the cumulative value test of the principal component variance contribution rate, and determine the principal component index data and the criterion weights of each principal component index; S3: Set multiple decision risk coefficient scenarios, calculate the rank weights and trade-off degrees of the principal component indicators according to the determined principal component index data, aggregate the rank weights and criterion weights of the principal component indicators under different scenarios using the OWA operator, and conduct multi-scenario simulation of the evaluation of ecological protection importance; S4: According to the simulation results of the evaluation of ecological protection importance, extract the ecological protection priority areas under various scenarios, compare the protection efficiency and trade-off degrees of the ecological protection priority areas in different scenarios for the evaluation results of ecological protection importance in each scenario, and screen out the ecological protection priority areas with the jointly optimal protection efficiency and trade-off degree; S5: Select the driving factors of land use change, construct a driving factor system of land use change with natural factors and human factors as the criterion layers, conduct the specific evaluation of the driving factors of land use change, and establish a spatial database of driving factors of land use change; S6: According to the historical land use data, obtain the target-year land use quantity structure and land use transfer probability matrix of the natural development scenario, extract the expansion patch data of each land use type according to the historical land use data, calculate the land use transfer probability matrix of the ecological protection scenario, and obtain the target-year land use quantity structure of the ecological protection scenario; S7: According to the spatial database of driving factors of land use change and the expansion patch data of each land use type, output the development probability data of each land use type, simulate the land use spatial structure data of the base year according to the historical land use data and the development probability data of each land use type, and conduct simulation accuracy verification; S8: Based on the land use data of the base year and the development probability data of each land use type, simulate the land use spatial structure data of the target year of the natural development scenario; based on the land use data of the base year and the development probability data of each land use type, simulate the land use spatial structure data of the target year of the ecological protection scenario, and delimit the flexible development boundary.

2. The method for delimiting the elastic development boundary of the urban agglomeration according to claim 1, characterized in that Specifically, step S1 is as follows: Select the evaluation indicators of ecological protection importance, and construct an evaluation index system of ecological protection importance with the evaluation of ecological protection importance as the target layer, and ecosystem services and ecological vulnerability as the criterion layer; combine multi-source data, ArcGIS software and InVEST model to conduct specific evaluations of the evaluation indicators of ecological protection importance; the evaluation results are uniformly normalized to eliminate the dimensional influence between different evaluation indicators; integrate all normalized evaluation indicators with the help of the ArcGIS platform to establish a spatial database of evaluation indicators of ecological protection importance in the study area. Among them, the supply service, regulation service, and support service belonging to ecosystem services, as well as the exposure, sensitivity, and adaptability belonging to ecological vulnerability are the first-level index layers, and the specific selected indicators are the second-level index layers of the evaluation index system of ecological protection importance. According to the influence on the target layer, the indicators belonging to adaptability among the selected specific indicators are negative indicators, and other indicators are positive indicators. The calculation method for normalizing positive indicators is: The calculation method for normalizing negative indicators is: where X norm is the normalized index pixel value, X x is the index pixel value, X max is the maximum value of all the index pixel values, and X min is the minimum value of all the index pixel values.

3. The method for delimiting the elastic development boundary of an urban agglomeration according to claim 1, characterized in that Step S2 is specifically as follows: Use the raster sampling tool of QGIS software to sample all the raster layers of evaluation indicators to a point layer converted from a raster with the resolution required by the research, and then import the attribute table data of the point layer into SPSS software for principal component analysis to generate principal component index data and a principal component analysis report; conduct KMO and Bartlett's tests, and use the eigenvalue greater than 1 as the standard for selecting principal components, and conduct a cumulative value test of the principal component variance contribution rate. Determine the principal component index data based on the evaluation indicators of ecological protection importance, and determine the criterion weights of each principal component index with a sum of 1 according to the principal component variance contribution rate. Among them, conducting KMO and Bartlett's tests is specifically as follows: If the KMO and Bartlett's test results in the principal component analysis report show that the KMO sampling adequacy measure is greater than 0.7 and the significance of the Bartlett sphericity test is less than 0.001, then the research data is suitable for principal component analysis; otherwise, the research data is not suitable for principal component analysis. Conducting a cumulative value test of the principal component variance contribution rate is specifically as follows: If the total variance explanation result in the principal component analysis report shows that the cumulative value of the principal component variance contribution rate with an eigenvalue greater than 1 exceeds 60%, then the generated principal component index data is relatively reasonable and effective; otherwise, the generated principal component index data is unreasonable and ineffective.

4. The method for delimiting the elastic development boundary of an urban agglomeration according to claim 1, wherein Step S3 is specifically as follows: Set z (z≥5) decision risk coefficient scenarios, and according to the determined principal component index data, use the monotonically increasing rule method to calculate the rank weights of the principal component indicators. The specific calculation method is as follows: where j is the rank; v j is the rank weight, v j ∈[0,1]; n is the number of principal component indicators; α is the decision-making risk coefficient, α∈(0,∞); w k is the importance level of the principal component indicator; r k is the assignment of the principal component indicator. The indicators are assigned according to the size after weighting by the pixel value superposition criterion of the indicators. The maximum value is assigned 1, the second largest value is assigned 2, and the minimum value is assigned n; The trade-off degree is the compensation degree between principal component indicators under different decision risk coefficients. The specific calculation method is as follows: where tradeoff is the tradeoff degree, and 0 ≤ tradeoff ≤ 1; n is the number of principal component indicators; w k is the importance level of the k-th principal component indicator; Use the OWA operator mechanism of the MCE module of TerrSet software to aggregate the rank weights and criterion weights of the principal component indicators under different scenarios, so as to obtain the spatial evaluation results of ecological protection importance in z scenarios.

5. The method for delimiting the elastic development boundary of an urban agglomeration according to claim 1, characterized in that Step S4 is specifically as follows: Based on the simulation results of the ecological protection importance evaluation under z scenarios, calculate the ratio of the cumulative ecological protection importance pixel value to the total ecological protection importance pixel value in descending order, and take the ecological protection importance pixel value corresponding to a ratio of 50% as the ecological protection importance threshold for each scenario; In the ArcGIS software, use the aggregation tool to aggregate relatively aggregated or adjacent patches in the area greater than or equal to the ecological protection importance threshold into relatively complete and contiguous patches. The aggregation distance is the resolution required for the study, and the aggregation result is used as the preliminary ecological protection priority area for each scenario; To reduce the fragmentation degree of the ecological protection priority area, use the patch of the ecological protection priority area in scenario z1 to determine the threshold of the patch area to be removed; among them, scenario z1 is the scenario where the ordinal weights of the main component indicators in the z scenarios are equal. Statistically analyze the ratio of the number of patches with an area smaller than the threshold to all patches in the ecological protection priority area of scenario z1, and the change of its total area with the area threshold. Take the minimum area threshold corresponding to the inflection point as the final patch area threshold, and remove the corresponding independent and scattered small patches according to the threshold. On this basis, merge the existing various ecological protection areas such as national and provincial nature reserves, forest parks, and scenic spots, and finally obtain the final ecological protection priority area for each scenario; Compare the protection efficiency and trade-off degree of the final ecological protection priority areas of different scenarios for the ecological protection importance evaluation results of scenario z1, so as to screen out the final ecological protection priority area of the scenario with the jointly optimal protection efficiency and trade-off degree; Among them, the calculation method of the protection efficiency is as follows: where z is the scenario number, z = 1, 2,..., z, P z is the protection efficiency of scenario z, is the average ecological protection importance of scenario z1 within the ecological protection priority area of scenario z, is the average ecological protection importance of scenario z1 within the entire study area.

6. The method for delimiting the elastic development boundary of an urban agglomeration according to claim 1, wherein Step S5 is specifically as follows: Select the driving factors of land use change, and construct a land use change driving factor system with natural factors and human factors as the criterion layer; combine multi-source data and the ArcGIS software to carry out the specific evaluation of the land use change driving factors; the obtained evaluation results are uniformly normalized to eliminate the dimensional influence between different driving factors; integrate all the normalized driving factors with the help of the ArcGIS platform to establish a spatial database of land use change driving factors in the study area; Among them, the terrain and water area belonging to natural factors and the population, economy, location, and infrastructure belonging to human factors are the first-level index layers, and the infrastructure includes transportation, medical care, education, and commercial services. The specific selected indicators are the land use change driving factor system of the second-level index layer.

7. The method for delimiting the elastic development boundary of the urban agglomeration according to claim 1, wherein, Step S6 is specifically as follows: Based on two-phase historical land use data, use the Markov Chain module of the PLUS model to obtain the land use quantity structure and land use transfer probability matrix of the target year in the natural development scenario; Based on two-phase historical land use data, use the land use expansion module of the PLUS model to extract the expansion patch data of each land use type. Based on this, calculate the proportion of each non-construction land in the expansion patch area of urban construction land, and use the obtained proportion data as the negative growth rate of the transfer probability of non-construction land to urban construction land. Thus, calculate the land use transfer probability matrix of the ecological protection scenario based on the land use transfer probability matrix of the natural development scenario; Based on the Markov Chain principle and ecological protection scenarios, calculate the land use quantity structure in the target year of the ecological protection scenario using the land use transfer probability matrix; Among them, the historical land use data includes seven land use types: cultivated land, forest land, grassland, water area, urban construction land, rural residential areas, and unused land. Among them, urban construction land and rural residential areas belong to construction land, and the rest of the land use types belong to non-construction land.

8. The method for delimiting the elastic development boundary of an urban agglomeration according to claim 1, wherein Step S7 is specifically as follows: Based on the land use change driving factor spatial database and the expansion patch data of each land use type, use the LEAS module of the PLUS model to output the development probability data of each land use type; calculate and set the parameters required for the CARS module of the PLUS model; Use the land use quantity structure of the second-phase historical land use data as the land use demand data. Based on the first-phase historical land use data and the development probability data of each land use type, use the CARS module of the PLUS model to simulate the land use spatial structure in the base year; based on the two-phase historical land use data and the base-year land use spatial structure data, use the Confusion Matrix and Fom module of the PLUS model to conduct simulation accuracy verification; Among them, calculating and setting the parameters required for the CARS module of the PLUS model is specifically as follows: Calculate the proportion of the expansion patch area of each land use type in the total expansion area based on the expansion patch data of each land use type, and use the obtained proportion data as the neighborhood weight parameter; set the allowable transfer value of each land use type in the land transfer matrix to 1, without restricting the transfer between any land use types; the attenuation coefficient of the decreasing threshold, the probability of the random patch seed, the maximum proportion of the random patch seed, the neighborhood range, and the number of parallel thread parameters use the default values; Conduct simulation accuracy verification, specifically as follows: If the Kappa coefficient is greater than 0.8 and the Fom coefficient is greater than 0.2, the simulation results of the PLUS model are basically consistent with the real results, and the simulation results are credible; otherwise, it is necessary to continuously adjust the input data and relevant parameters of the PLUS model until the above Kappa coefficient and Fom coefficient test requirements are met before proceeding with the subsequent steps.

9. The method for delimiting the elastic development boundary of an urban agglomeration according to claim 1, wherein Step S8 is specifically as follows: Use the land use quantity structure in the target year of the natural development scenario as the land use demand data, and the existing various ecological protection areas as restricted development areas. Based on the base-year land use data and the development probability data of each land use type, use the CARS module of the PLUS model to simulate the land use spatial structure data in the target year of the natural development scenario; Use the land use quantity structure in the target year of the ecological protection scenario as the land use demand data, and the ecological protection priority areas as restricted development areas. Based on the base-year land use data and the development probability data of each land use type, use the CARS module of the PLUS model to simulate the land use spatial structure data in the target year of the ecological protection scenario; The urban construction land boundaries in the land use spatial structure data of the natural development scenario and the ecological protection scenario target years are respectively designated as the elastic development boundaries of the natural development scenario and the ecological protection scenario target years in the study area.

10. A system for delimiting the elastic development boundary of an urban agglomeration, characterized in that, Applying the method for delimiting the elastic development boundary of the urban agglomeration according to any one of claims 1 to 9, comprising: An evaluation index module for establishing a spatial database of ecological protection importance evaluation indexes; A principal component analysis module for performing principal component analysis based on the spatial database of ecological protection importance evaluation indexes; An OWA operator module for simulating the ecological protection importance evaluation of multiple scenarios based on the OWA operator; A scenario screening module for extracting the ecological protection priority areas of each scenario and screening out the ecological protection priority areas of the optimal scenario; A driving factor module for establishing a spatial database of land use change driving factors; A demand prediction module for predicting the land use quantity structure of the natural development scenario and the ecological protection scenario target years; A simulation and verification module for simulating the land use spatial structure of the base year and verifying the simulation accuracy; A boundary delimitation module for simulating the land use spatial structure of the natural development scenario and the ecological protection scenario target years and delimiting the elastic development boundary.

Images