Land utilization pattern evaluation model construction method based on dynamic climate and application

Abstract

The invention relates to the technical field of land space pattern optimization, in particular to a land utilization pattern evaluation model construction method based on dynamic climate and application. The invention provides a dynamic climate-based land utilization pattern evaluation model construction method and application, and the method achieves the effective evaluation of the influence of the land utilization change on the near-surface climate through the coupling of an SD-CA model and a near-surface climate dynamic correlation model LCDR. According to the method, land utilization changes in different scenes can be simulated, the land utilization changes can be fed back to the SD model through the climate correlation module, land utilization pattern configuration is optimized, and influences brought by climate changes are relieved. In addition, the method also considers factors such as environmental factors, neighborhood influence, inertia transition probability and random disturbance, so that the model is more comprehensive and accurate. By implementing the method, sustainable utilization of land resources can be effectively guided, and a scientific basis and technical support are provided for coping with climate changes.

Description

Technical Field

[0001] The present invention relates to the technical field of land space pattern optimization, and specifically to a method for constructing a land use pattern assessment model based on dynamic climate and its application. Background Art

[0002] Land use and its pattern changes caused by human activities, such as the expansion of built-up areas, the contiguous aggregation of impervious surfaces, and the fragmentation of ecological land, not only put tremendous pressure on the ecological environment and food security on the surface, but also have an impact on climate change that has been confirmed in many studies. Therefore, how to build a land use pattern that is friendly to climate conditions and mitigate climate change caused by land use patterns has become an important issue that needs to be solved urgently.

[0003] Among the existing land use models, the cellular automation model (CA) is the most widely used. This model regards the minimum land grid as a cell, and each cell evolves according to the rules. At the same time, the change in land use demand is the driving force of the spatiotemporal change of land use. The overall land use classification is used for quantitative control. The system dynamics model (SD) can abstract the impact of the complex relationship between socio-economic elements on land use and form a connection between socio-economic development and land use change. Although many of these models have simulated the impact of socio-economic and geographical factors on land use change, there is no effective method to construct a feedback loop between the change in land use and the change in near-surface climate and the land use pattern, and to establish a method to mitigate the impact of land use pattern change on surface climate. Summary of the invention

[0004] The technical problem solved by the present invention is to provide a method for constructing a land use pattern assessment model based on dynamic climate and its application so as to feed back the changes in near-surface climate caused by land use changes to the land use pattern and mitigate the impact of land use pattern changes on the surface climate.

[0005] The basic scheme provided by the present invention is a method for constructing a land use pattern assessment model based on dynamic climate, comprising the steps of:

[0006] S100: Establish the system boundary and construct the land use quantity demand system dynamic model based on the population subsystem, economic industry subsystem and land use subsystem;

[0007] S200: Construct the environmental impact probability, domain impact probability, inertial conversion probability and random disturbance probability of the cellular automaton model to simulate the spatial change of land use;

[0008] S300: The SD-CA model is formed by coupling the system dynamics model and the cellular automaton model, and the quantitative simulation results of each land on the corresponding time scale are obtained according to the quantitative structural changes of land use types under different scenarios;

[0009] S400: Construct a dynamic correlation model between land use pattern and near-surface climate (LCDR) to quantify the correlation between the dynamic impact of land use change on near-surface climate;

[0010] S500: A climate correlation constraint module is constructed based on the dynamic correlation model between land use pattern and near-surface climate to quantify the correlation between land use pattern changes and climate indicators. The climate correlation constraint module couples this correlation to the SD-CA model through the constraints of dynamic climate factors on land use pattern, thus realizing the land use spatiotemporal simulation optimization model SD-CA-CCOM for mitigating climate change.

[0011] Further, the S100 includes the steps of:

[0012] S110: Establish system boundaries including spatial boundaries, temporal boundaries, and content boundaries:

[0013] S120: Construct the causal relationship among the population subsystem, economic industry subsystem and land use subsystem respectively, establish a system dynamic model according to the causal relationship among the subsystems, and generate a fitting function between the parameters that change over time.

[0014] Further, the S110 includes the steps of:

[0015] S111: The external matrix of the study area is used as the spatial boundary;

[0016] S112: The simulation start and end times are used as time boundaries, the time when data is input is used as the start time, and the time required for simulation is used as the end time;

[0017] S113: The population subsystem, economic industry subsystem and land use subsystem are used as the content boundaries of the entire system.

[0018] Further, the S120 includes the steps of:

[0019] S121: The causal relationship of population change in the population subsystem is:

[0020]

[0021] In the formula, Pop t is the population base; is the natural increase of population, Natural increase of population For birth population and deaths The difference between is the mechanical growth rate of population, Mechanical growth of population For immigrants With the out-migration population The difference between

[0022] S122: The economic industry subsystem includes the industrial labor force coefficient divided by the three industries, the labor force divided by the three industries, the per capita output value coefficient divided by the three industries, and the industrial output value divided by the three industries. The gross national product can be obtained by summing up the output values ​​of the three industries. The causal relationship of the gross national product is:

[0023] GDP t =FGDP t +IGDP t +SGDP t

[0024] In the formula, GDP at time t is composed of the output value of the primary industry in the current period FGDP t , Secondary Industry Output Value IGDP t and the output value of the tertiary industry SGDP t The sum of FGDP t =FWPop×FWGDPC, FWPop is the labor force in the primary industry, FWGDPC is the output value coefficient of the primary industry; IGDP t =IWPop×IWGDPC, where IWPop is the labor force in the secondary industry, and IWGDPC is the output value coefficient of the secondary industry; SGDP t =SWPop×SWGDPC, where SWPop is the labor force in the tertiary industry and SWGDPC is the output value coefficient of the tertiary industry;

[0025] S123: In the land use subsystem, land use is divided into six types: cultivated land, forest land, grassland, water, construction land and unused land. The causal relationship between their areas is as follows:

[0026] TotalAREA=ACPL+AFST+AGSL+AWTR+AUBL+ABAL

[0027] In the formula, the total area TotalAREA determined by the spatial boundary is a constant value, the cultivated land area ACPL is the regression value of RURPop, agricultural investment coefficient CPLINVESTR, and grain yield per unit area GRAINYP, where RURPop is the rural population obtained by the urbanization rate of Enjia; the forest area AFST is the regression value of RURPop, forestry output value LYGDP and forestry investment coefficient FSTINVESTR; the grassland area AGSL is the regression value of RURPop, unit animal husbandry output value HUSYP and animal husbandry investment coefficient GSLINEVSTR; the water area AWTR is replaced by a table function that changes over time; the construction land area AUBL is the regression value of GDP and fixed asset investment coefficient UBLINVESTR; the unused land area ABAL is set as the difference between TotalAREA and various other land areas;

[0028] S124: Establish a system dynamics model based on the cause-effect relationship of each subsystem, and generate a fitting function of the time-varying parameters.

[0029] Further, the S200 includes the steps of:

[0030] S210: Constructing the environmental impact probability of cellular automaton model:

[0031] S211: All social-natural factor data are converted into Lambert projection coordinate system, and based on Arcpy environment, the Euclidean distance from each raster to each social-natural factor is calculated through Euclidean Distance function; the social-natural factors include transportation connection factors such as roads (vector line data) and railways (vector line data), economic connection factors such as town centers (vector point data) and banks (vector point data), social connection factors such as hospitals (vector point data) and schools (vector point data), natural connection factors such as rivers (vector line data) and natural scenic spots (vector point data), and terrain factors such as DEM (raster data) and slope (raster data);

[0032] S212: Logistic regression method is used to establish the transition probability PE of each central cell's surrounding environment affecting the conversion of the cell's land use type into various land use types i,j,m :

[0033] PE i,j,m =[p1,p2,…p k …,p n ]

[0034] Where, PE i,j,m The probability p that the land use type of cell (i, j) is m and is transformed into the kth land use type due to environmental factors is k ;

[0035] The probability p for each transition type n Calculated by Logistic regression model:

[0036]

[0037] This formula is the Sigmoid function expression in the Logistic model;

[0038] y=β1X1+β2X2+…+β n X n +γ

[0039] In the formula, β is the weight corresponding to the environmental factor X, where X is the social-natural factor, and γ is the regression constant. The weight of each cell converted into other land use types is calculated based on this, and then the transfer probability of the cell type to other land use types is calculated based on the weight. When m=k, it means that the current land use type remains unchanged;

[0040] Before performing logistic regression, the values ​​of each environmental factor were normalized to remove dimensional differences. The normalization algorithm used the following formula:

[0041]

[0042] In the formula, It represents the normalized value of the factor X corresponding to the i-th cell, Min(X) represents the minimum value of all X factors, and Max(X) represents the maximum value of all X factors.

[0043] S220: Constructing domain impact probability of cellular automaton model:

[0044] S221: Define the probability of the central cell transforming to this type of land use affected by the proportion of the neighboring land use type:

[0045] PN i,j,m =[η1,η2,…η k …,η n ]

[0046] Where PN i,j,m The probability p that the land use type of cell (i, j) is m and is transformed into the kth land use type due to the influence of neighboring cells is k , because it is determined by the type of neighborhood cells defined by the central cell (i, j), so n is uncertain, but the land use type is classified into 6 categories, so n≤6;

[0047] According to the neighborhood window size N (N ≥ 1), the probability of Moore's neighborhood influence can be defined as:

[0048]

[0049] Where η k It means that the central cell is attracted by the transformation of the kth neighboring cell, and W(k) represents the number of k-type land use types under the neighborhood window;

[0050] S230: Constructing the inertial transition probability of the cellular automaton model:

[0051] S231: Based on the Markov chain method, the land use transition matrix is ​​used to show the probability of inertial change in land use change.

[0052]

[0053] The probability of a certain type transforming into another specific type;

[0054]

[0055] Where N mk represents the area converted from land type m to land type k between time t-1 and t, It represents the sum of all land types transformed from land type m during this time. This value is the total area of ​​the region. The probability value of each land type change

[0056]

[0057] In the formula, It represents the total number of land types converted to land type k during this period. This value is also the total area of ​​the region. The probability value of each land type change

[0058] S240: Constructing random perturbation probability of cellular automaton model, random perturbation probability PR m,t :

[0059] PR m,t =1+(-lnθ) α

[0060] In the formula, is a random number between 0 and 1, and α is a positive integer between 1 and 10, which is used as the intensity coefficient to control the random disturbance. Thus, the transition probability P of each cell is obtained. m,t :

[0061] P m,t =PE i,j,m ×PN i,j,m ×PS mk ×PR m,t

[0062] Where Pm,t It represents the probability of converting the cell state m to other land use types at time t, according to the conversion probability P of each cell land use type m,t In the CA simulation, all cells are iterated simultaneously without any order, and the result of the iteration is the spatial distribution of the probability of land use transformation.

[0063] Further, the S300 includes the steps of:

[0064] S301: Starting from SD, driven by the macro subsystem variables in SD, the quantitative simulation results of each land use in the land use subsystem on a time scale under different scenarios are obtained according to different population changes, economic development trends, land use changes, etc.;

[0065] S302: taking the quantity ratio of each land use as a macro constraint to act on the land use transformation rule function of the CA model, controlling the number of spatial iterations of the CA model to affect the neighborhood state, environmental factors, random disturbances and cellular inertia probability, and thus changing the land use type of each cell;

[0066] S303: Count the land use quantity ratio on the spatial scale of CA simulation, and adjust the starting value of the dynamic equation of the SD model at the next moment to drive the simulation of the SD system at the next moment;

[0067] S304: Through the bidirectional interactive coupling method of S301-S303, the CA simulation calculation formula is rewritten as follows:

[0068] P m,t =F[f(PE i,j,m,t ,PN i,j,m,t ,PS mk,t ,PR m,t ),SD t+1 ]

[0069] In the formula, the transition probability P of each cell at time t is m,t , not only affected by the corresponding environmental factor PE at time t i,j,m,t , Neighborhood state PN i,j,m,t , Cellular inertia PS mk,t and random perturbation PR m,t The influence of t+1 is also affected by the parameter variable SD of the SD model simulation corresponding to the next time t+1. t+1 impact.

[0070] Further, the S400 includes the steps of:

[0071] S401: Calculate the centrality index of the land landscape pattern index through the principle of graph network betweenness centrality, and select the index items with the highest centrality index of several index categories; then use the sensitivity experiment method to test the range in which the index of the index item weakens with the landscape grain size, and determine the grain size that is insensitive to the grain size change as the appropriate landscape grain size;

[0072] S402: Setting a simulation scheme for near-surface climate response under land use pattern change, and simulating and obtaining near-surface climate data under the influence of the corresponding land use pattern;

[0073] S403: According to the index of the indicator item at the appropriate landscape granularity and the near-surface climate data under the influence of the land use pattern, the least square fitting coefficient of the landscape pattern index and the corresponding climate data is obtained:

[0074]

[0075] In the formula, represents the near-surface climate factor corresponding to grid ij; is the land landscape index corresponding to grid ij; β is the land landscape index corresponding to grid ij and The regression coefficient is the correlation between the land landscape index and the climate factor on the ij grid; ε is the residual as the intercept; the above fitting is evolved into the corresponding relationship between the change of the land landscape index and the change of the climate factor:

[0076]

[0077] In the formula, and The WRF simulated near-surface climate factors corresponding to the grid ij are the next time results and the starting time results of the study. and is the landscape pattern index result at the next time corresponding to grid ij and the result at the start time of the study, and β′ is the dynamic correlation coefficient between the two.

[0078] Furthermore, the climate research and forecasting model is adopted in S402, and the land use raster data corresponding to the simulation time is replaced into the climate research and forecasting model to form the underlying surface data formed by the land use spatial pattern, and then the seven parameters in the climate research and forecasting model including the microphysical equation, the longwave radiation equation, the shortwave radiation equation, the underlying surface parameter scheme, the surface parameter scheme, the atmospheric boundary layer scheme and the cumulus scheme are set to simulate and obtain the near-surface climate data under the influence of the corresponding land use pattern.

[0079] Furthermore, in S500, the CCOM module optimizes the land use pattern configuration and makes modifications during the CA calculation process, adding the following modifications:

[0080]

[0081] In the formula, It is the normalized calculation result of the dynamic correlation coefficient DR of climate factors corresponding to the land use type m and landscape pattern index r calculated by the LCDR model. The smaller the DR is, the smaller the impact on the climate is, the smaller the constraint is, and the greater the probability of conversion is allowed, so the DR is negative; is the dynamic change of land use pattern, is the difference between the land use pattern at time t and the previous time t-1, which can be expressed as follows:

[0082]

[0083] In the formula, LM m,r,t represents the land landscape index r, w for land use type m at time t F is the calculation weight of the original SD-CA, w l is the calculation weight of the CCOM coupling module. If w F >w l This means that the model is more inclined to the original natural resource and socio-economic driven model. F <w l It is more inclined to optimize the model.

[0084] The application of the dynamic climate-based land use pattern assessment model, wherein any of the above-mentioned dynamic climate-based land use pattern assessment models is applied, comprises the steps of:

[0085] S1: The SDT1 variable value at time T1 is obtained by simulating the starting data of the SD model, and is passed to the CA model according to the SD-CA connection variable;

[0086] S2: Based on the SDT1 variable value, the initial land use grid data T0 is guided to generate the original SD-CA simulated land use grid data T1 through the CA model;

[0087] S3: Based on the land use grid data at time T1 and T0, the land landscape index of the two times is calculated, and then the difference of the land use patterns of the two phases is calculated to obtain the Δ spatial pattern;

[0088] S4: Based on the calculation of the dynamic correlation coefficient between land use pattern and climate factors and the Δ spatial pattern, the environmental constraints of climate factors are integrated to obtain the optimized land use grid data T1';

[0089] S5: Based on the land use grid data of T1', the number of each land use type is counted and passed into the SD model as the starting data after the SD model is fine-tuned at time T2. The process in S1 is iterated to calculate the land use spatial pattern at the next time.

[0090] The principle and advantage of the present invention are as follows: the present invention provides a method and application for constructing a land use pattern assessment model based on dynamic climate, by establishing the system boundary and constructing a land use quantity demand system dynamic model according to the population subsystem, the economic industry subsystem and the land use subsystem, then constructing the environmental impact probability, the domain impact probability, the inertia conversion probability and the random disturbance probability of the cellular automaton model, simulating the spatial change of land use, and by coupling the SD-CA model with the near-surface climate dynamic correlation model LCDR, an effective assessment of the impact of land use change on the near-surface climate is achieved. The method can not only simulate the land use changes under different scenarios, but also feed back the land use changes to the SD model through the climate correlation module, optimize the land use pattern configuration, and mitigate the impact of climate change. In addition, the invention also takes into account factors such as environmental factors, neighborhood influences, inertia conversion probability and random disturbances, making the model more comprehensive and accurate. By implementing the invention, the sustainable use of land resources can be effectively guided, and a scientific basis and technical support can be provided for responding to climate change. BRIEF DESCRIPTION OF THE DRAWINGS

[0091] Figure 1 The present invention is a flowchart of the steps of an embodiment of a method for constructing a land use pattern assessment model based on dynamic climate.

[0092] Figure 2 A flow diagram of a socio-economic-land use change system dynamics model according to an embodiment of a method for constructing a land use pattern assessment model based on dynamic climate of the present invention.

[0093] Figure 3 A map of environmental factors affecting land use changes in the Beijing-Tianjin-Hebei urban agglomeration according to an embodiment of a method for constructing a land use pattern assessment model based on dynamic climate in the present invention.

[0094] Figure 4 This is a probability map of the impact of natural environmental factors on different land use types in the Beijing-Tianjin-Hebei urban agglomeration according to an embodiment of a method for constructing a land use pattern assessment model based on dynamic climate of the present invention.

[0095] Figure 5 A coupling interaction mechanism diagram of the SD-CA model in an embodiment of a method for constructing a land use pattern assessment model based on dynamic climate in the present invention.

[0096] Figure 6 This is a land shape landscape index LIBC correlation matrix diagram of an embodiment of a method for constructing a land use pattern assessment model based on dynamic climate in the present invention.

[0097] Figure 7This is a multi-granularity change and multi-scale change slope diagram of the land landscape index of the Beijing-Tianjin-Hebei urban agglomeration in an embodiment of a method for constructing a land use pattern assessment model based on dynamic climate in the present invention.

[0098] Figure 8 The present invention discloses a dynamic correlation table of land landscape pattern indexes according to an embodiment of a method for constructing a land use pattern assessment model based on dynamic climate.

[0099] Fig. 9 It is a schematic diagram of the operation iteration of the SD-CA-CCOM model of an embodiment of a method for constructing a land use pattern assessment model based on dynamic climate in the present invention.

[0100] Fig.10 This is a comparison chart between the SD-CA-CCOM simulation results of land use and the actual land use results of an embodiment of a method for constructing a land use pattern assessment model based on dynamic climate in the present invention.

[0101] Fig.11 This is a NST difference comparison diagram of an embodiment of a method for constructing a land use pattern assessment model based on dynamic climate in the present invention.

[0102] Fig.12 This is a logical diagram of the operation of the SD-CA-CCOM model of an embodiment of a method for constructing a land use pattern assessment model based on dynamic climate in the present invention. DETAILED DESCRIPTION

[0103] The following is further described in detail through specific implementation methods:

[0104] The specific implementation process is as follows:

[0105] Embodiment 1

[0106] Embodiment 1 is as attached Figure 1 As shown, a method for constructing a land use pattern assessment model based on dynamic climate includes the following steps:

[0107] S100: Establish the system boundary and construct the land use quantity demand system dynamic model based on the population subsystem, economic industry subsystem and land use subsystem;

[0108] S200: Construct the environmental impact probability, domain impact probability, inertial conversion probability and random disturbance probability of the cellular automaton model to simulate the spatial change of land use;

[0109] S300: The SD-CA model is formed by coupling the system dynamics model and the cellular automaton model, and the quantitative simulation results of each land on the corresponding time scale are obtained according to the quantitative structural changes of land use types under different scenarios;

[0110] S400: Construct a dynamic correlation model between land use pattern and near-surface climate (LCDR) to quantify the correlation between the dynamic impact of land use change on near-surface climate;

[0111] S500: A climate correlation constraint module is constructed based on the dynamic correlation model between land use pattern and near-surface climate to quantify the correlation between land use pattern changes and climate indicators. The climate correlation constraint module couples this correlation to the SD-CA model through the constraints of dynamic climate factors on land use pattern, thus realizing the land use spatiotemporal simulation optimization model SD-CA-CCOM for mitigating climate change.

[0112] Wherein, S100 comprises the steps of:

[0113] S110: Establish system boundaries including spatial boundaries, temporal boundaries, and content boundaries.

[0114] S120: Construct the causal relationship among the population subsystem, economic industry subsystem and land use subsystem respectively, establish a system dynamic model according to the causal relationship among the subsystems, and generate a fitting function between the parameters that change over time.

[0115] The SD model is used to abstractly and simplifies the impact of socio-economic development on various types of land use, as well as the feedback of changes in the quantitative demand for various types of land use on socio-economic development, and to simulate the changes in the quantitative structure of various land use types under different future development scenarios.

[0116] Specifically, the establishment of the system boundary in S110 includes:

[0117] S111: The external matrix of the study area is used as the spatial boundary.

[0118] S112: The simulation start and end times are used as time boundaries, the time when data is input is used as the start time, and the time required for simulation is used as the end time; the data is iterated every year.

[0119] S113: Population subsystem, economic industry subsystem and land use subsystem are used as the content boundaries of the whole system. Because the social population and economic industry are selected as the main driving factors of urban land use change through literature reference and Delphi method, and the change of land use quantity ratio is the result, the system content boundaries are the above three subsystems.

[0120] The population subsystem, economic industry subsystem and land use subsystem in this embodiment have causal relationships between their parameters. Specifically, the derivation process of the causal relationship in S120 is:

[0121] Population subsystem: Population changes are the basis for influencing economic growth and the quantitative structure of land use. They promote economic development through the urban-rural ratio and the flow of labor in various industries, and change the quantitative structure of land use. Construct parameters including total population (RKZL), birth rate (BIRRATE), death rate (DIERATE), in-migration population (IMGINPop) and out-migration population (IMGOUTPop). The birth population (BIRPop) and IMGINPop form the population increment (INCPop), while the death population (DIEPop) and IMGOUTPop form the population decrement (DECPop). The difference between the two and the sum of the total population in the base period is the total population each year. Therefore, the causal relationship of the total population change in the population subsystem is:

[0122]

[0123] In the formula, Pop t is the population base; is the natural increase of population, Natural increase of population For birth population and deaths The difference between is the mechanical growth rate of population, Mechanical growth of population For immigrants With the out-migration population The output of the population subsystem includes the labor force (LABPop) obtained according to the labor force coefficient, which will affect the labor force and gross domestic product of each industry; the rural population (RURPop) obtained according to the urbanization rate will affect the cultivated land area, forest area and grassland area in the land use quantity structure.

[0124] Economic industry subsystem: Changes in investment in economic industries will indirectly or directly drive the transformation of various land use structures. It is an economic activity system formed by the development of different industries through labor in human society. It is the reflection and basis of human activities on the utilization and transformation of different land types. The economic industry subsystem in this embodiment includes the industrial labor force coefficient divided by the three industries, the labor force divided by the three industries, the per capita output value coefficient divided by the three industries, and the industrial output value divided by the three industries. The gross national product can be obtained by summing the output values ​​of the three industries. The causal relationship of the gross national product is:

[0125] GDP t =FGDP t +IGDP t +SGDP t

[0126] In the formula, GDP at time t is composed of the output value of the primary industry in the current period FGDP t , Secondary Industry Output Value IGDP t and the output value of the tertiary industry SGDP t The sum of FGDP t =FWPop×FWGDPC, FWPop is the labor force in the primary industry, FWGDPC is the output value coefficient of the primary industry; IGDP t =IWIPop×IWGDPC, where IWPop is the labor force in the secondary industry, and IWGDPC is the output value coefficient of the secondary industry; SGDP t =SWPop×SWGDPC, where SWPop is the labor force in the tertiary industry and SWGDPC is the output value coefficient of the tertiary industry. In the output of the economic industry subsystem, FGDP t According to the proportional coefficient, it is divided into agricultural output value (NYGDP), forestry output value (LYGDP) and animal husbandry output value (XMYGDP), which will further affect the cultivated land area, forest land area and grassland area. At the same time, the growth of GDP will drive the construction land area through the fixed asset investment coefficient, and then drive the transformation of the land use subsystem.

[0127] Land use subsystem: The urbanization process has caused the transformation of urban population structure and the improvement of economic and industrial structure. It has also accelerated the expansion of construction land, the reduction of cultivated land, the degradation of forest land and other changes in land use pattern. In this subsystem, the development of industrial output value, the changes in the investment coefficients of various industries, the flow of rural population and other indicators will affect the main land use changes. However, the total area of ​​each type of land use is a fixed value. Therefore, in the land use subsystem, land use is divided into 6 types of land use: cultivated land, forest land, grassland, water, construction land and unused land. The causal relationship between their areas is as follows:

[0128] TotalAREA=ACPL+AFST+AGSL+AWTR+AUBL+ABAL

[0129] In the formula, the total area TotalAREA determined by the spatial boundary is a constant value, the cultivated land area ACPL is the regression value of RURPop, agricultural investment coefficient CPLINVESTR, and grain yield per unit area GRAINYP, where RURPop is the rural population obtained by the urbanization rate of Enjia; the forest area AFST is the regression value of RURPop, forestry output value LYGDP and forestry investment coefficient FSTINVESTR; the grassland area AGSL is the regression value of RURPop, unit animal husbandry output value HUSYP and animal husbandry investment coefficient GSLINEVSTR; the water area AWTR is replaced by a table function that changes over time; the construction land area AUBL is the regression value of GDP and fixed asset investment coefficient UBLINVESTR; the unused land area ABAL is set as the difference between TotalAREA and various other land areas;

[0130] Coupling construction of each subsystem: Based on the cause-effect relationship, the system dynamics model is established using the system dynamics principle to link the three subsystems. The mathematical relationship between parameters is fitted using existing statistical data as prior data to form a fitting function that changes over time between parameters. The main forms of the function include linear, exponential, logarithmic, sine and cosine functions, etc. The highest fitting accuracy is used as the selection criterion. Over time, all parameters in this system will change with the change of time and other parameters.

[0131] In this example, the Beijing-Tianjin-Hebei urban agglomeration is used as the research area. The SD-CA-CCOM model is used to simulate the land use data of the Beijing-Tianjin-Hebei urban agglomeration from 2000 to 2020 to obtain the simulated data for 2020, which is verified with the real land use data for 2020 to test the accuracy of the model. Based on the data from the Social and Economic Statistical Yearbook, the required data items in the population subsystem, economic industry subsystem, and land use subsystem are obtained. According to the positive and negative correlations between the relevant data items, the following is established: Figure 2 The SD model shown has a mathematical relationship determined based on the fit coefficients of the existing data.

[0132] Since there are 38 data items involved, based on the data of the Social and Economic Statistical Yearbook, the data from 2000 to 2016 are the most complete, so this study takes the data from 2000 to 2016 as historical reference data (Table 1), as the initial value of the SD model, and makes parameter corrections to obtain trends close to the actual social and economic development and verify the model. The data obtained include the following variables in the model:

[0133] Table 1 Variable data and related descriptions in the system dynamics model

[0134]

[0135]

[0136] In S200, the environmental impact probability, domain impact probability, inertial conversion probability and random disturbance probability of the cellular automaton model are constructed to simulate the spatial change of land use.

[0137] Among them, the CA model abstracts the rules of spatial land use change, constructs environmental impact rules, neighborhood impact rules, inertial transformation effects and random disturbance factors, simulates the spatial change of land use, and defines the transformation rules as:

[0138] Where P m,t It indicates the probability of converting the cell state m to other land use types at time t, which is determined by the environmental impact probability PE m,t , Neighborhood influence probability PN m,t , inertial conversion probability PS m,t and the random perturbation probability PR m,t The result of a joint impact.

[0139] Specifically, S200 includes the following steps:

[0140] S210: Constructing the environmental impact probability of cellular automaton model:

[0141] S211: All social-natural factor data are converted into Lambert projection coordinate system to ensure the standardized measurement of distance, and based on Arcpy environment, the Euclidean distance from each raster to each social-natural factor is calculated through Euclidean Distance function; the social-natural factors include transportation connection factors such as roads (vector line data) and railways (vector line data), economic connection factors such as town centers (vector point data) and banks (vector point data), social connection factors such as hospitals (vector point data) and schools (vector point data), natural connection factors such as rivers (vector line data) and natural scenic spots (vector point data), and terrain factors such as DEM (raster data) and slope (raster data);

[0142] S212: Logistic regression method is used to establish the weight of the surrounding environment of each central cell affecting the conversion of the land use type of the cell into various land use types. By determining the distance of each factor and the binary determination of each land use classification, the probability of the environmental factors of each cell being converted into various land use types is obtained. The transition probability PE of each cell affected by environmental factors i,j,m :

[0143] PE i,j,m =[p1,p2,…p k …,p n ]

[0144] Where, PE i,j,m The probability p that the land use type of cell (i, j) is m and is transformed into the kth land use type due to environmental factors is k ; There are 6 land use classifications, so n=6 in this embodiment.

[0145] The probability p for each transition type n Calculated by Logistic regression model:

[0146]

[0147] This formula is the Sigmoid function expression in the Logistic model;

[0148] y=β1X1+β2X2+…+β n X n +γ

[0149] In the formula, β is the weight corresponding to the environmental factor X, where X is the social-natural factor, and γ is the regression constant. The weight of each cell converted into other land use types is calculated based on this, and then the transfer probability of the cell type to other land use types is calculated based on the weight. When m=k, it means that the current land use type remains unchanged. In this embodiment, binary logistic regression is used, so when calculating for each type of land use, this type of land use is set to 1, and the other land use values ​​are set to 0. For example, if the cell is construction land, it is set to 1, and the rest are set to 0 to indicate non-construction land.

[0150] Before performing logistic regression, the values ​​of each environmental factor were normalized to remove dimensional differences. The normalization algorithm used the following formula:

[0151]

[0152] In the formula, It represents the normalized value of the factor X corresponding to the i-th cell, Min(X) represents the minimum value of all X factors, and Max(X) represents the maximum value of all X factors.

[0153] S220: Constructing domain impact probability of cellular automaton model:

[0154] S221: In this embodiment, in a window, if the proportion of the surrounding land use types is higher for the central cell, the attraction to the central cell is greater, resulting in a greater probability that the central cell will be converted to this land type. Therefore, the probability of the central cell being converted to this type of land use under the influence of the proportion of the neighboring land use types is defined as:

[0155] PN i,j,m=[η1,η2,…η k …,η n ]

[0156] Where PN i,j,m The probability p that the land use type of cell (i, j) is m and is transformed into the kth land use type due to the influence of neighboring cells is j , because it is determined by the type of neighborhood cells defined by the central cell (i, j), so n is uncertain, but the land use type is classified into 6 categories, so n≤6;

[0157] According to the neighborhood window size N (N ≥ 1), the probability of Moore's neighborhood influence can be defined as:

[0158]

[0159] Where η j Indicates that the central cell is attracted by the transformation of the kth neighboring cell, W(k) indicates the number of k-type land use types under the neighborhood window; in a matrix such as (2N+1)×(2N+1), all cell types except the central cell are summed. Because of the principle of spatial autocorrelation, this study uses the tightest 3×3 window, that is, N=1.

[0160] S230: Constructing the inertial transition probability of the cellular automaton model:

[0161] S231: Because there is also the influence of inertial change probability in land use change, in order to adjust the inheritance probability of each cellular land use type, a land use transfer matrix is ​​used in this embodiment. It is based on the Markov chain method and not only quantitatively measures the number of transfers between land use types, but also can clearly and intuitively display the change rate between land use types.

[0162]

[0163] Where, PS m,k The set matrix represents the probability of land type m at time t-1 transforming into land type k at time t, and ω represents the probability of a specific type transforming into another specific type; 1q represents the probability of the first land type being transformed into the qth land type, ω qq Represents the probability that the land type remains unchanged during the changing time period. ω mk The calculation of is obtained by the following formula:

[0164]

[0165] Where N mk represents the area converted from land type m to land type k between time t-1 and t, It represents the sum of all land types transformed from land type m during this time. This value is the total area of ​​the region. The probability value of each land type change

[0166]

[0167] In the formula, It represents the sum of the land types converted from other land types to land type k during this time. This value is also the total area of ​​the region. The probability value of each land type change

[0168] S240: Random disturbance can just cover the land use changes caused by uncertain factors; on the other hand, in the CA model, some land use types will be overfitted quickly and the simulation iteration will end. Through random disturbance probability, the CA model can achieve the best effect for multiple land use types in the simulation. Therefore, the random disturbance probability of the cellular automaton model is constructed, and the random disturbance probability PR m,t :

[0169] PR m,t =1+(-lnθ) α

[0170] In the formula, is a random number between 0 and 1, and α is a positive integer between 1 and 10, which is used as the intensity coefficient to control the random disturbance. Therefore, because the influence of each probability is superimposed rather than accumulated, the conversion rule function is defined as the product of each probability to obtain the conversion probability P of each cell. m,t :

[0171] P m,t =PE i,j,m ×PN i,j,m ×PS mk ×PR m,t

[0172] Where P m,t It represents the probability of converting the cell state m to other land use types at time t, according to the conversion probability P of each cell land use type m,t In the CA simulation, all cells are iterated simultaneously without any order, and the result of the iteration is the spatial distribution of the probability of land use transformation.

[0173] In this embodiment, the spatial data required by the CA model include land use data, POI dataset, road network data and elevation data. According to the corresponding data requirements in the second part S1(1) of the method, the land use data adopts the image data of 2000 and 2020 of the Resource and Environmental Science and Data Center of the Chinese Academy of Sciences (RESDC, https: / / www.resdc.cn / ); the POI dataset adopts the Amap POI dataset, and the latest data is the 2020 dataset; the road network data adopts the 2018 OSM road network data; the elevation data adopts the DEM data of STRM_90m.

[0174] The 15 land use environmental factors required were extracted from the above data sets to obtain the distribution map of environmental factors affecting land use change in the Beijing-Tianjin-Hebei urban agglomeration, as shown in the figure. Figure 3 shown.

[0175] According to the land use environmental change factors, the regression values ​​of the impact of environmental factors on different land use types are obtained through the Logistic regression method, such as Figure 4 shown.

[0176] The probability of land use conversion is based on the RESDC land use data of 2000 and 2020. The probability of land use conversion between time nodes is established based on the land use matrix as the inertial conversion probability of land use conversion. The inertial conversion probability table is shown in Table 2 below.

[0177] Table 2 Probability of land use type transfer in the Beijing-Tianjin-Hebei urban agglomeration from 2000 to 2020 (%)

[0178] arable land woodland grassland Water Construction land Unused land arable land 95.424 0.094 0.093 0.252 4.065 0.073 woodland 0.094 99.369 0.054 0.019 0.465 0 grassland 0.059 0.125 98.717 0.134 0.95 0.016 Water 1.159 0.153 1.892 86.789 8.802 1.205 Construction land 0.519 0.038 0.132 0.033 99.264 0.014 Unused land 0.064 0.128 0.385 7.628 1.731 90.064 Unused land 0.064 0.128 0.385 7.628 1.731 90.064

[0179] In S300, the system dynamics model and the cellular automation model are coupled to form the SD-CA model, and the quantitative simulation results of each land on the corresponding time scale are obtained according to the quantitative structure changes of land use types under different scenarios.

[0180] Specifically, this method couples the interaction mechanism between SD and CA models as follows Figure 1 As shown, in order to enable the SD and CA models to feedback each other, the parameter variables of the SD and CA connection transmit the simulation results at each moment to each other as cause and effect, so as to realize the coupling of SD-CA, so S300 includes the steps:

[0181] S301: Starting from SD, driven by the macro subsystem variables in SD, the quantitative simulation results of each land use in the land use subsystem on a time scale under different scenarios are obtained according to different population changes, economic development trends, land use changes, etc.;

[0182] S302: taking the quantity ratio of each land use as a macro constraint to act on the land use transformation rule function of the CA model, controlling the number of spatial iterations of the CA model to affect the neighborhood state, environmental factors, random disturbances and cellular inertia probability, and thus changing the land use type of each cell;

[0183] S303: Count the land use quantity ratio on the spatial scale of CA simulation, and adjust the starting value of the dynamic equation of the SD model at the next moment to drive the simulation of the SD system at the next moment;

[0184] S304: Through the bidirectional interactive coupling method of S301-S303, the CA simulation calculation formula is rewritten as follows:

[0185] P m,t =F[f(PE i,j,m,t ,PN i,j,m,t ,PS mk,t ,PR m,t ),SD t+1 ]

[0186] In the formula, the transition probability P of each cell at time t is m,t , not only affected by the corresponding environmental factor PE at time t i,j,m,t , Neighborhood state PN i,j,m,t , Cellular inertia PS mk,t and random perturbation PR m,t The influence of t+1 is also affected by the parameter variable SD of the SD model simulation corresponding to the next time t+1. t+1 impact.

[0187] In this embodiment, the interaction mechanism is as follows: Figure 5 As shown. The construction process of the coupled SD-CA model can be divided into the following 6 steps: 1) clarify the model boundaries and objectives; 2) clarify the interactive variables of the SD-CA model; 3) construct the SD model; 4) construct the CA model; 5) interactive communication of SD-CA variables; 6) adjust the model parameters. This study uses Python as the programming language for model implementation. The SD model part uses the Vensim platform to sort out and construct the variable relationship, and calls PySD to construct the mathematical equations between variables; the data processing of the CA model part is mainly processed by Arcpy extended by the ArcGIS platform, and the calculation and simulation of the probability of land use cell conversion also call GDAL, NumPy, Pandas, Math, Statsmodels and other Python libraries.

[0188] In S400, based on the research on static correlation, a land use pattern and near-surface climate dynamic correlation quantitative model (Landscape and Climate Dynamic Relevance model, LCDR) was proposed, which consists of three parts: land use change landscape calculation, near-surface climate environment dynamic simulation, and dynamic correlation fitting, in order to quantify the correlation degree of the dynamic impact of land use change on the near-surface climate.

[0189] Specifically, S400 includes the following steps:

[0190] S401: Calculate the centrality index of the land landscape pattern index through the principle of graph network betweenness centrality, and select the index items with the highest centrality index of several index categories; then use the sensitivity experiment method to test the range in which the index of the index item weakens with the landscape grain size, and determine the grain size that is insensitive to grain size changes as the appropriate landscape grain size.

[0191] Since there are multiple similar land landscape pattern indices, among which different landscape pattern indices have similarities, the most representative land landscape pattern index should be selected to reduce repeated calculations of similarity indices. The principle of graph network betweenness centrality is introduced into this method to calculate the Pearson correlation coefficient of each land landscape pattern index. The correlation coefficients of a land landscape pattern index and other indices are accumulated, and the sum is the centrality index (LIBC) of the index. The higher the centrality index, the more representative it is. This method selects the landscape indexes in three types of indicators, namely, land landscape area size, landscape shape index, and landscape aggregation, and finally selects the indicator item with the highest centrality index of the three types of indicators as the subsequent calculation item.

[0192] In addition, because there is a significant scale effect in the land use pattern, the land landscape index will also change with the change of grain size, and the scale effect will affect the index's description of the landscape pattern in the region, so it is necessary to screen out specific research scales for different research areas. This method selects the sensitivity experiment method, that is, by testing the range in which each landscape index weakens with the landscape grain size, so as to determine the grain size that is insensitive to grain size changes as the appropriate model grain size.

[0193] S402: Setting a simulation scheme for near-surface climate response under changes in land use pattern, and simulating near-surface climate data under the influence of the corresponding land use pattern. In this embodiment, a climate research and forecasting mode is adopted, and the land use grid data corresponding to the simulation time is replaced into the climate research and forecasting mode to form underlying surface data formed by the land use spatial pattern, and then seven parameters including microphysics equation (Microphysics), longwave radiation equation (Longwave radiation), shortwave radiation equation (Shortwave radiation), underlying surface parameter scheme (Surface-layer option), surface parameter scheme (Land-surface option), atmospheric boundary layer scheme (Boundary-layer option) and cumulus option (Cumulus option) in the climate research and forecasting mode are set to simulate and obtain near-surface climate data under the influence of the corresponding land use pattern.

[0194] S403: According to the index of the indicator item at the appropriate landscape granularity and the near-surface climate data under the influence of the land use pattern, the least square fitting coefficient of the landscape pattern index and the corresponding climate data is obtained:

[0195]

[0196] In the formula, represents the near-surface climate factor corresponding to grid ij; is the land landscape index corresponding to grid ij; β is the land landscape index corresponding to grid ij and The regression coefficient is the correlation between the land landscape index and the climate factor on the ij grid; ε is the residual as the intercept; the above fitting is evolved into the corresponding relationship between the change of the land landscape index and the change of the climate factor:

[0197]

[0198] In the formula, and The WRF simulated near-surface climate factors corresponding to the grid ij are the next time results and the starting time results of the study. and is the landscape pattern index result at the next time corresponding to grid ij and the result at the start time of the study, and β′ is the dynamic correlation coefficient between the two.

[0199] In this embodiment, the remote sensing image classification data corresponding to the Beijing-Tianjin-Hebei urban agglomeration in 2000 is calculated using Fragstats4.2 software, and the remote sensing image classification data corresponding to 2000 is preferably determined according to the calculation process of each sub-index in the spatial pattern index to determine the concentration landscape pattern sub-index, shape landscape pattern sub-index, and area / edge landscape pattern sub-index. Among them, these three sub-indices characterize the composition pattern and spatial pattern of the land use pattern, and generate the following Table 3 (Land Area / Edge Landscape Index LIBC Correlation Matrix), Table 4 (Land Shape Landscape Index LIBC Correlation Matrix) and Table 5 (Land Shape Landscape Index LIBC Correlation Matrix) Landscape index correlation matrix, and the representative index is selected by the landscape index betweenness centrality (the betweenness centrality value LIBC of the sub-index).

[0200] Table 3 Land area / edge landscape index LIBC correlation matrix

[0201] Area_MN LPI ED GYRATE_MN LIBC AREA_MN 1.000 0.980 0.655 0.918 3.553 LPI 0.980 1.000 0.536 0.944 3.460 ED 0.655 0.536 1.000 0.487 2.678 GYRATE_MN 0.918 0.944 0.487 1.000 3.349

[0202] Table 4 LIBC correlation matrix of land shape landscape index

[0203] FRAC_MN SHAPE_MN PARA_MN CIRCLE_MN CONTIG_MN PAFRAC LIBC FRAC_MN 1.000 0.975 -0.783 0.629 0.803 0.594 4.785 SHAPE_MN 0.975 1.000 -0.824 0.760 0.844 0.547 4.950 PARA_MN -0.783 -0.824 1.000 -0.691 -0.996 0.372 4.666 CIRCLE_MN 0.629 0.760 -0.691 1.000 0.682 0.176 3.938 CONTIG_MN 0.803 0.844 -0.996 0.682 1.000 -0.210 4.534 PAFRAC 0.594 0.547 0.372 0.176 -0.210 1.000 2.898

[0204] From the above Table 3, Table 4 and Figure 6 It can be seen that among the calculated area / edge landscape pattern sub-indexes, the calculated index AREA_MN has the largest LIBC (3.553); among the shape landscape pattern indexes, the calculated index SHAPE_MN has a large LIBC (4.950); among the aggregation landscape pattern sub-indexes, the calculated index AI has the largest LIBC (7.437). In the spatial pattern index, the average patch area (AREA_MN) is used to characterize the area / edge landscape pattern index, the average shape index (SHAPE_MN) is used to characterize the shape landscape pattern index index, and the aggregation (AI) is used to characterize the aggregation landscape pattern index. In the composition pattern index, the proportion of the class area in the landscape unit (PLAND) index is directly used to characterize the composition pattern index, which represents the proportion of the area occupied by a landscape class in the patch. PLAND does not belong to the spatial pattern index. In order to characterize the land use pattern from different dimensions, specifically, the above four landscape pattern indices are selected. The detailed explanation of the four indices is shown in Table 6 (Description of the representative index of land landscape in the Beijing-Tianjin-Hebei urban agglomeration).

[0205] Table 6 Description of the representative index of land landscape in the Beijing-Tianjin-Hebei urban agglomeration

[0206]

[0207]

[0208] To determine the suitability granularity of the land use pattern in the study area, a sensitivity experiment method is generally used. Specifically, the method for determining the preset granularity in the above embodiment is used to calculate the various granularities of the four selected landscape pattern indices, PLAND, Area_MN, SHAPE_MN and AI, for example, the 11 granularities of n = 3, 5, 7, 9, 11, 15, 21, 25, 31, 41, 51 (km), to determine the granularity interval where the landscape pattern index tends to be stable with the change of granularity. Since the intervals of granularity are not completely equal, the rate of change (slope) of the index is further calculated, that is, the derivative of the granularity change and the corresponding index. Please refer to the fitting curves (a, b, c, d) and the derivative curves (e, f, g, h). Figure 6 , determine the particle size n = 11km as the appropriate particle size.

[0209] After simulation and calibration, the physical parameterization scheme of the WRF model for the Beijing-Tianjin-Hebei urban agglomeration is shown in Table 7.

[0210] Table 7 Physical parameterization scheme setting table of WRF model of Beijing-Tianjin-Hebei urban agglomeration

[0211] parameter plan Solution parameter value Microphysics equations WSM6 Solution 6 Longwave Radiation Equation RRTM Solution 1 Shortwave Radiation Equation Dudhia Solution 1 Surface layer options Revised MM5 Monin-Obukhov scheme 1 Land surface options Noah land surface plan 2 Boundary Layer Options YSU Program 1 Cumulus Options Kain-Fritsch scheme 1

[0212] The improved least squares model was used to fit the correlation between the change in land use pattern and the change in near-surface climate factors, and the dynamic correlation coefficient between the land use landscape pattern and near-surface climate factors was obtained. Taking the grid as the statistical unit, the index change value of the grid where the land use pattern changed and the NST change value of different months were extracted. The least squares model was used to calculate the dynamic correlation coefficient (DR) of each land landscape pattern index in January and July. The results are as follows: Figure 6 As shown, PR is the Pearson correlation coefficient.

[0213] In S500, a climate correlation constraint module is constructed based on the dynamic correlation model between land use pattern and near-surface climate to quantify the correlation between land use pattern changes and climate indicators. The climate correlation constraint module couples this correlation to the SD-CA model through the constraints of dynamic climate factors on land use pattern, thus realizing the land use spatiotemporal simulation optimization model SD-CA-CCOM for mitigating climate change.

[0214] Specifically, in S500, a climate correlation constraint model CCOM (ClimateConstraint Optimization Model) with a dynamic correlation model of land use pattern and near-surface climate as the core is constructed on the basis of LCDR. The correlation between land use pattern changes and climate indicators is quantified according to the LCDR results. CCOM feeds this correlation back to the land use pattern simulation, and couples the constraints of land use pattern by dynamic climate factors to the land use model integrated with SD-CA, thus realizing the land use spatiotemporal simulation optimization model SD-CA-CCOM for climate change mitigation. The CCOM module optimizes the land use pattern configuration and makes modifications during the CA calculation process, adding the following modifications:

[0215]

[0216] In the formula, It is the normalized calculation result of the dynamic correlation coefficient DR of climate factors corresponding to the land use type m and landscape pattern index r calculated by the LCDR model. The smaller the DR is, the smaller the impact on the climate is, the smaller the constraint is, and the greater the probability of conversion is allowed, so the DR is negative; is the dynamic change of land use pattern, is the difference between the land use pattern at time t and the previous time t-1, which can be expressed as follows:

[0217]

[0218] In the formula, LM m,e,r represents the land landscape index r, w for land use type m at time t F is the calculation weight of the original SD-CA, w l is the calculation weight of the CCOM coupling module. If w F >w l This means that the model is more inclined to the original natural resource and socio-economic driven model. F <w l It is more inclined to optimize the model.

[0219] In this embodiment, the model continuously iterates the S100-S500 process, and finally stops the model simulation when the simulation time of SD exceeds the research time. The specific iterative operation process is as follows: Figure 7 As shown, each iteration goes through 4 cycles (marked in red in the figure).

[0220] Loop 1: It is a mutually driven loop between SD and CA. The SD model simulates the socio-economic data from the starting time t, and simulates the number of land use types at t′ through the socio-demographic subsystem-economic-industrial subsystem-land use subsystem, guiding the CA model to iterate and enter loop 2. When the number of land use spatial statistics simulated by the CA model reaches the simulated value of the SD model at t′, the value is returned to the SD model, t′=t+1, and then returns to the starting point of loop 1 for the next loop.

[0221] Loop 2: This is the internal loop of the CA model. The t′ value calculated by the SD model is used as the control threshold. The CA model continuously iterates the spatial pattern of land use to meet the land use driving factor value of the SD model, forming a land use spatial pattern corresponding to t′. Through the statistics of the number of classifications, the SD land use driving factor is passed to the SD model for the next cycle.

[0222] Loop 3: It is the LCDR coupling loop of the CCOM climate constraint optimization module. The land use spatial patterns at t′ and t generated by the CA model drive the WRF model to solve the climate change caused by the change of land use spatial pattern. At the same time, the land use spatial pattern landscape index is solved. The DR coefficient at t′ is fitted by the LCDR model to form a spatial constraint, which is integrated into the CA model for spatial optimization and promotes the next cycle of CA.

[0223] Loop 4: It is a sub-loop of the CCOM climate constraint optimization module. When the land use spatial patterns generated by the CA model at time t′ and time t are slightly different, the impact on climate factors is small. At the same time, the WRF solution process takes a relatively long time. According to the experiment, this study sets the time interval to 5, that is, it enters loop 4 only when t′-t=5 to generate a new DR coefficient. When t′-t<5, the spatial constraint is generated by calculating the landscape pattern index and the DR coefficient at time t to promote CA to the next cycle.

[0224] In summary, the above four cycles are operated, and the model simulation ends when t+1=tEND. The generated land use prediction data will be combined with the GIS analysis model for spatiotemporal evaluation and analysis.

[0225] This study conducted a historical test on the simulation of the SD model. If the simulation results differ greatly from the actual values, the parameters need to be adjusted. The model will be improved and the simulation results close to the actual results and consistent with the change trend will be obtained through debugging. First, the three key socio-economic factor variables of the Beijing-Tianjin-Hebei urban agglomeration, the total population, the gross regional product (GDP), and the urbanization rate are selected for comparison. The real statistical values, simulation values, and their error percentages from 2000 to 2016 are shown in Table 9 below. From the results, it can be seen that the error between the simulation value and the historical real value is within a controllable range, all within 7%; only the population and GDP around 2012 have large errors, about 6%; the urbanization rate, as a rate variable that changes over time, is not affected by other variable factors, and is only obtained by mathematical abstraction based on the historical real value, so it is close to the real value and has a small error.

[0226] Table 9 Comparison of historical real values ​​and system dynamics simulation values ​​of socio-economic variables

[0227]

[0228]

[0229] Secondly, this study conducted a simulation test on the number of six land use types. Because there is no continuous data for land use classification data, the classification data of RESDC every 5 years was used as standard data for comparison, and the SD simulation of land use classification was tested using the data from 2000 / 2005 / 2010 / 2015. The comparison results are shown in Table 10. The results show that the error between the SD model and the historical true value is within 5%. It is generally believed that an error <5% indicates that the model can reflect the changes in the true value, has good consistency with the historical trend, and can effectively predict future changes.

[0230] Table 10 Comparison between historical real values ​​and system dynamics simulation values ​​of land use variables

[0231]

[0232] The SD-CA-CCOM model with the CCOM module optimizes the spatial pattern of land use, and uses near-surface temperature as a reference climate factor to test the effect of mitigating climate warming. This section verifies the historical development of land use change and its climate effect to ensure that SD-CA-CCOM can meet the real process of land use change on the one hand, and optimize the spatial pattern of land use to achieve the effect of mitigating climate warming on the other hand.

[0233] The land use data in 2000 were used as the first phase of data, and the SD-CA-CCOM model was used to simulate the land use in 2020. The simulation results were compared with the actual land use data in 2020, as shown below. Figure 8 As shown in the figure, the comparison shows that the SD-CA-CCOM model can better reflect the spatial evolution of land use. The Kappa coefficient of the simulation result is 0.784 and the overall accuracy is 0.823.

[0234] By extracting the 2m near-surface air temperature of the simulation results and making the difference, that is, SD-CA simulation 2020-real 2020 and SD-CA-CCOM simulation 2020-real 2020, the difference between the land use of the model before and after optimization and the climate conditions driven by the real land use is obtained. The spatial difference is as follows Fig. 9 As shown in the figure, it can be found that the temperature difference changes are mainly concentrated in the plain areas in the east and south of the Beijing-Tianjin-Hebei urban agglomeration. These areas are also the areas where land use changes mainly occurred from 2000 to 2020, especially the areas where construction land increased and concentrated. The difference is that the SD-CA simulation formed dense warming areas in the south of Beijing and the northeast of Shijiazhuang, and there were concentrated cooling areas in the southeast coastal areas, which is quite different from the actual simulation results. Compared with the difference of SD-CA-CCOM, except for a few warming areas concentrated in the north of Beijing and a few cooling areas concentrated in the southeast, the rest are alternating between warming and cooling, indicating that SD-CA-CCOM can reduce the contiguous warming and concentrated areas through the spatial pattern of land use in land use optimization through the constraints of spatial pattern.

[0235] The WRF simulation results of the minimum, maximum and average values ​​of NST corresponding to the simulation results before and after land use optimization in the Beijing-Tianjin-Hebei urban agglomeration in 2020 are shown in Table 11. From the NST changes corresponding to the real land use in the table, it can be found that the land use changes in 2020 compared with 2000 caused the minimum and maximum values ​​of NST to decrease in January (-0.027℃ and -0.121℃), while the average value increased (0.05℃); the corresponding minimum, maximum and average values ​​in July all showed an upward trend (0.051℃, 0.286℃ and 0.177℃). The simulation results of SD-CA and SD-CA-CCOM in this study are consistent with the NST changes corresponding to the real land use. SD-CA increased the average temperature in January by 0.111℃, while the optimized SD-CA-CCOM only increased by 0.083℃, which is 0.05℃ higher than the real land use in 2020. In July of summer, the unoptimized SD-CA caused the average temperature to rise by 0.179℃, which was higher than the 0.177℃ caused by the actual 2020 land use, while the optimized SD-CA-CCOM caused the average temperature to rise by 0.142℃, which was lower than the warming caused by the actual 2020 land use, and the minimum change (0.041℃) and maximum change (0.281℃) were also lower than the warming caused by the actual 2020 (0.051℃ and 0.39℃).

[0236] Table 11 Maximum and average NST values ​​corresponding to land use in SD-CA and SD-CA-CCOM Unit: °C

[0237]

[0238] In summary, although the NST corresponding to the land use simulated by SD-CA-CCOM after landscape pattern optimization is slightly higher than the NST corresponding to the real land use in January, it is lower than the NST corresponding to the land use simulated by SD-CA without optimization; and for the optimization results in July, the NST corresponding to the land use simulated by SD-CA-CCOM is lower than the NST corresponding to the real and SD-CA simulated land use. This shows that the SD-CA-CCOM model has the effect of alleviating climate warming through the optimization of land use pattern.

[0239] Embodiment 2

[0240] Embodiment 2 is basically as shown in the attached Fig.12 As shown, the difference between the second embodiment and the first embodiment is that the second embodiment is an application of a land use pattern assessment model based on dynamic climate, and any of the above land use pattern assessment models based on dynamic climate is applied. The SD-CA model generates the land use spatial pattern at the next time by CA based on the land use spatial structure and SD results at the initial time. According to the difference in the land use spatial pattern at two times, combined with the dynamic correlation between the land landscape index calculated in the CCOM module and climate change, the climate and environmental constraints of the CA model are generated, land use is optimized, and then a land use pattern that takes into account climate change is generated. It specifically includes the following steps:

[0241] S1: The SDT1 variable value at time T1 is obtained by simulating the starting data of the SD model, and is passed to the CA model according to the SD-CA connection variable;

[0242] S2: Based on the SDT1 variable value, the initial land use grid data T0 is guided to generate the original SD-CA simulated land use grid data T1 through the CA model;

[0243] S3: Based on the land use grid data at time T1 and T0, the land landscape index of the two times is calculated, and then the difference of the land use patterns of the two phases is calculated to obtain the Δ spatial pattern;

[0244] S4: Based on the calculation of the dynamic correlation coefficient between land use pattern and climate factors and the Δ spatial pattern, the environmental constraints of climate factors are integrated to obtain the optimized land use grid data T1';

[0245] S5: Based on the land use grid data of T1', the number of each land use type is counted and passed into the SD model as the starting data after the SD model is fine-tuned at time T2. The process in S1 is iterated to calculate the land use spatial pattern at the next time.

[0246] The above are only embodiments of the present invention. Common knowledge such as the known specific structures and characteristics in the scheme is not described in detail here. Ordinary technicians in the relevant field know all the common technical knowledge in the technical field to which the invention belongs before the application date or priority date, can obtain all the existing technologies in the field, and have the ability to apply conventional experimental means before that date. Ordinary technicians in the relevant field can improve and implement this scheme in combination with their own abilities under the enlightenment given by this application. Some typical known structures or known methods should not become obstacles for ordinary technicians in the relevant field to implement this application. It should be pointed out that for those skilled in the art, without departing from the structure of the present invention, several deformations and improvements can be made, which should also be regarded as the scope of protection of the present invention, and these will not affect the effect of the implementation of the present invention and the practicality of the patent. The scope of protection required by this application shall be based on the content of its claims, and the specific implementation methods and other records in the specification can be used to interpret the content of the claims.

Claims

1. A method for constructing a land use pattern assessment model based on dynamic climate, characterized in that: Includes steps: S100: Establish the system boundary and construct the land use quantity demand system dynamic model based on the population subsystem, economic industry subsystem and land use subsystem; S200: Construct the environmental impact probability, domain impact probability, inertial conversion probability and random disturbance probability of the cellular automaton model to simulate the spatial change of land use; S300: The SD-CA model is formed by coupling the system dynamics model and the cellular automaton model, and the quantitative simulation results of each land on the corresponding time scale are obtained according to the quantitative structural changes of land use types under different scenarios; S400: Construct a dynamic correlation model between land use pattern and near-surface climate (LCDR) to quantify the correlation between the dynamic impact of land use change on near-surface climate; S5 00: Based on the dynamic correlation model between land use pattern and near-surface climate, a climate correlation constraint module is constructed to quantify the correlation between land use pattern changes and climate indicators. The climate correlation constraint module couples this correlation to the SD-CA model through the constraints of dynamic climate factors on land use pattern, thus realizing the land use spatiotemporal simulation optimization model SD-CA-CCOM for mitigating climate change.

2. The method for constructing a land use pattern assessment model based on dynamic climate according to claim 1, characterized in that: The S100 comprises the steps of: S110: Establish system boundaries including spatial boundaries, temporal boundaries, and content boundaries: S120: Construct the causal relationship among the population subsystem, economic industry subsystem and land use subsystem respectively, establish a system dynamic model according to the causal relationship among the subsystems, and generate a fitting function between the parameters that change over time.

3. The method for constructing a land use pattern assessment model based on dynamic climate according to claim 2, characterized in that: The S110 includes the following steps: S111: The external matrix of the study area is used as the spatial boundary; S112: The simulation start and end times are used as time boundaries, the time when data is input is used as the start time, and the time required for simulation is used as the end time; S113: The population subsystem, economic industry subsystem and land use subsystem are used as the content boundaries of the entire system.

4. The method for constructing a land use pattern assessment model based on dynamic climate according to claim 3, characterized in that: The S120 includes the following steps: S121: The causal relationship of population change in the population subsystem is: In the formula, Popt is the population base; is the natural increase of population, Natural increase of population For birth population and deaths The difference between is the mechanical growth rate of population, Mechanical growth of population For immigrants With the out-migration population The difference between S122: The economic industry subsystem includes the industrial labor force coefficient divided by the three industries, the labor force divided by the three industries, the per capita output value coefficient divided by the three industries, and the industrial output value divided by the three industries. The gross national product can be obtained by summing up the output values ​​of the three industries. The causal relationship of the gross national product is: GDP t =FGDP t +IGDP t +SGDP t In the formula, GDP at time t is composed of the output value of the primary industry in the current period FGDP t , Secondary Industry Output Value IGDP t and the output value of the tertiary industry SGDP t The sum of FGDP t =FWPop×FWGDPC, FWPop is the labor force in the primary industry, FWGDPC is the output value coefficient of the primary industry; IGDP t =IWPop×IWGDPC, where IWPop is the labor force in the secondary industry, and IWGDPC is the output value coefficient of the secondary industry; SGDP t =SWPop×SWGDPC, where SWPop is the labor force in the tertiary industry and SWGDPC is the output value coefficient of the tertiary industry; S123: In the land use subsystem, land use is divided into six types: cultivated land, forest land, grassland, water, construction land and unused land. The causal relationship between their areas is as follows: TotalAREA=ACPL+AFST+AGSL+AWTR+AUBL+ABAL In the formula, the total area TotalAREA determined by the spatial boundary is a constant value, the cultivated land area ACPL is the regression value of RURPop, agricultural investment coefficient CPLINVESTR, and grain yield per unit area GRAINYP, where RURPop is the rural population obtained by the urbanization rate of Enjia; the forest area AFST is the regression value of RURPop, forestry output value LYGDP and forestry investment coefficient FSTINVESTR; the grassland area AGSL is the regression value of RURPop, unit animal husbandry output value HUSYP and animal husbandry investment coefficient GSLINEVSTR; the water area AWTR is replaced by a table function that changes over time; the construction land area AUBL is the regression value of GDP and fixed asset investment coefficient UBLINVESTR; the unused land area ABAL is set as the difference between TotalAREA and various other land areas; S124: Establish a system dynamics model based on the cause-effect relationship of each subsystem, and generate a fitting function of the time-varying parameters.

5. The method for constructing a land use pattern assessment model based on dynamic climate according to claim 4, characterized in that: The S200 includes the steps of: S210: Constructing the environmental impact probability of cellular automaton model: S211: All social-natural factor data are converted into Lambert projection coordinate system, and based on Arcpy environment, the Euclidean distance from each raster to each social-natural factor is calculated through Euclidean Distance function; the social-natural factors include transportation connection factors such as roads (vector line data) and railways (vector line data), economic connection factors such as town centers (vector point data) and banks (vector point data), social connection factors such as hospitals (vector point data) and schools (vector point data), natural connection factors such as rivers (vector line data) and natural scenic spots (vector point data), and terrain factors such as DEM (raster data) and slope (raster data); S212: Logistic regression method is used to establish the transition probability PE of each central cell's surrounding environment affecting the conversion of the cell's land use type into various land use types i,j,m : INSTEAD i,j,m [p1,p2,...p k ...,p n ] Where, PE i,j,m The probability p that the land use type of cell (i, j) is m and is transformed into the kth land use type due to environmental factors is k ; The probability p for each transition type n Calculated by Logistic regression model: This formula is the Sigmoid function expression in the Logistic model; y=β1X1+β2X2+…+β n X n +g In the formula, β is the weight corresponding to the environmental factor X, where X is the social-natural factor, and γ is the regression constant. The weight of each cell converted into other land use types is calculated based on this, and then the transfer probability of the cell type to other land use types is calculated based on the weight. When m=k, it means that the current land use type remains unchanged; Before performing logistic regression, the values ​​of each environmental factor were normalized to remove dimensional differences. The normalization algorithm used the following formula: In the formula, represents the normalized value of the factor X corresponding to the i-th cell, Min(X) represents the minimum value of all X factors, and Max(X) represents the maximum value of all X factors; S220: Constructing domain impact probability of cellular automaton model: S221: Define the probability of the central cell transforming to this type of land use affected by the proportion of the neighboring land use type: PN i,j,m =[η1,η2,...η k ...,or n ] Where PN i,j,m The probability p that the land use type of cell (i, j) is m and is transformed into the kth land use type due to the influence of neighboring cells is k , because it is determined by the type of neighborhood cells defined by the central cell (i, j), so n is uncertain, but the land use type is classified into 6 categories, so n≤6; According to the neighborhood window size N (N ≥ 1), the probability of Moore's neighborhood influence can be defined as: Where η k It means that the central cell is attracted by the transformation of the kth neighboring cell, and W(k) represents the number of k-type land use types under the neighborhood window; S230: Constructing the inertial transition probability of the cellular automaton model: S231: Based on the Markov chain method, the land use transfer matrix is ​​used to display the probability of inertial change in land use change. The transfer matrix of the land use type from the quantitative status at time t-1 to the quantitative status at time t can be expressed as: Where, PS m,k The set matrix represents the probability of land type m at time t-1 transforming into land type k at time t, and ω represents the probability of a specific type transforming into another specific type; Where N mk represents the area converted from land type m to land type k between time t-1 and t, It represents the sum of all land types transformed from land type m during this time. This value is the total area of ​​the region. The probability value of each land type change In the formula, It represents the sum of the land types converted from other land types to land type k during this time. This value is also the total area of ​​the region. The probability value of each land type change S240: Constructing the random perturbation probability of the cellular automaton model, random perturbation probability PR m,t : PR m,t =1+(-lnθ) α In the formula, is a random number between 0 and 1, and α is a positive integer between 1 and 10, which is used as the intensity coefficient to control the random disturbance. Thus, the transition probability P of each cell is obtained. m，t : P m，t =PE i,j,m ×PN i，j，m ×PS mk ×PR m,t Where P m，t It represents the probability of converting the cell state m to other land use types at time t, according to the conversion probability P of each cell land use type m，t , according to the conversion probability P of each cell land use type m，t In the CA simulation, all cells are iterated simultaneously without any order, and the result of the iteration is the spatial distribution of the probability of land use transformation.

6. The method for constructing a land use pattern assessment model based on dynamic climate according to claim 5, characterized in that: The S300 includes the following steps: S301: Starting from SD, driven by the macro subsystem variables in SD, the quantitative simulation results of each land use in the land use subsystem on a time scale under different scenarios are obtained according to different population changes, economic development trends, land use changes, etc.; S302: taking the quantity ratio of each land use as a macro constraint to act on the land use transformation rule function of the CA model, controlling the number of spatial iterations of the CA model to affect the neighborhood state, environmental factors, random disturbances and cellular inertia probability, and thus changing the land use type of each cell; S303: Count the land use quantity ratio on the spatial scale of CA simulation, and adjust the starting value of the dynamic equation of the SD model at the next moment to drive the simulation of the SD system at the next moment; S304: Through the bidirectional interactive coupling method of S301-S303, the CA simulation calculation formula is rewritten as follows: P m，t =F[f(PE i,j,m,t ,PN i,j,m,t ,P.S. mk，t ,PR m,t ), SD t+1 ] In the formula, the transition probability P of each cell at time t is m，t , not only affected by the corresponding environmental factor PE at time t i,j,m,t , Neighborhood state PN i，j，m，t , Cellular inertia PS mk，t and random perturbation PR m,t The influence of t+1 is also affected by the parameter variable SD of the SD model simulation corresponding to the next time t+1. t+1 impact.

7. The method for constructing a land use pattern assessment model based on dynamic climate according to claim 6, characterized in that: The S400 step: S401: Calculate the centrality index of the land landscape pattern index through the principle of graph network betweenness centrality, and select the index items with the highest centrality index of several index categories; then use the sensitivity experiment method to test the range in which the index of the index item weakens with the landscape grain size, and determine the grain size that is insensitive to the grain size change as the appropriate landscape grain size; S402: Setting a simulation scheme for near-surface climate response under land use pattern change, and simulating and obtaining near-surface climate data under the influence of the corresponding land use pattern; S403: According to the index of the indicator item at the appropriate landscape granularity and the near-surface climate data under the influence of the land use pattern, the least square fitting coefficient of the landscape pattern index and the corresponding climate data is obtained: In the formula, represents the near-surface climate factor corresponding to grid ij; is the land landscape index corresponding to grid ij; β is the land landscape index corresponding to grid ij and The regression coefficient is the correlation between the land landscape index and the climate factor on the ij grid; ε is the residual as the intercept; the above fitting is evolved into the corresponding relationship between the change of the land landscape index and the change of the climate factor: In the formula, and The WRF simulated near-surface climate factors corresponding to the grid ij are the next time results and the starting time results of the study. and is the landscape pattern index result at the next time corresponding to grid ij and the result at the start time of the study, and β′ is the dynamic correlation coefficient between the two.

8. The method for constructing a land use pattern assessment model based on dynamic climate according to claim 7, characterized in that: In the S402, a climate research and forecasting model is adopted to replace the land use grid data corresponding to the simulation time into the climate research and forecasting model to form underlying surface data formed by the spatial pattern of land use, and then seven parameters including microphysical equations, longwave radiation equations, shortwave radiation equations, underlying surface parameter schemes, surface parameter schemes, atmospheric boundary layer schemes and cumulus schemes in the climate research and forecasting model are set to simulate and obtain near-surface climate data under the influence of the corresponding land use pattern.

9. The method for constructing a land use pattern assessment model based on dynamic climate according to claim 8, characterized in that: In S500, the CCOM module optimizes the land use pattern configuration and makes modifications during the CA calculation process, adding the following modifications: In the formula, It is the normalized calculation result of the dynamic correlation coefficient DR of climate factors corresponding to the land use type m and landscape pattern index r calculated by the LCDR model. The smaller the DR is, the smaller the impact on the climate is, the smaller the constraint is, and the greater the probability of conversion is allowed, so the DR is negative; is the dynamic change of land use pattern, is the difference between the land use pattern at time t and the previous time t-1, which can be expressed as follows: Where, LM m，r，t represents the land landscape index r, w for land use type m at time t F is the calculation weight of the original SD-CA, w l is the calculation weight of CC0M coupling module, if W F >w l This means that the model is more inclined to the original natural resource and socio-economic driven model. F <w l It is more inclined to optimize the model.

10. Application of a dynamic climate-based land use pattern assessment model, characterized by: The land use pattern assessment model based on dynamic climate described in any one of claims 1 to 9 is applied, comprising the steps of: S1: The SDT1 variable value at time T1 is obtained by simulating the starting data of the SD model, and is passed to the CA model according to the SD-CA connection variable; S2: Based on the SDT1 variable value, the initial land use grid data T0 is guided to generate the original SD-CA simulated land use grid data T1 through the CA model; S3: Based on the land use grid data at time T1 and T0, the land landscape index of the two times is calculated, and then the difference of the land use patterns of the two phases is calculated to obtain the Δ spatial pattern; S4: Based on the calculation of the dynamic correlation coefficient between land use pattern and climate factors and the Δ spatial pattern, the environmental constraints of climate factors are integrated to obtain the optimized land use grid data T1'; S5: Based on the land use grid data of T1', the number of each land use type is counted and passed into the SD model as the starting data after the SD model is fine-tuned at time T2. The process in S1 is iterated to calculate the land use spatial pattern at the next time.

Images