A method, medium and device for urban land use change simulation

Abstract

The application discloses a kind of urban land use change simulation method, medium and equipment, it is related to urban land use change simulation field, the urban land use change simulation method is combined with a variety of pre-training time sequence model and cellular automaton, constructs optimal time sequence model to calculate the conversion probability array of each plot, the overall development probability array of each plot is obtained by combining traditional vector cellular automaton, the final conversion probability array of each plot is obtained, finally, the land use class of each plot conversion is obtained by merging. The urban land use change simulation method provided by the application can take into account the time dependence of urban land use change, improve the accuracy of urban land use change simulation.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of urban land use change simulation, and more particularly, to a method, medium and device for simulating urban land use change. BACKGROUND

[0002] Exploring the mechanism of land use change can guide the healthy and sustainable development of cities, but the consideration of time sequence factors in land use change still needs to be explored. Mining the time sequence change rule of land use helps to deeply understand the internal mechanism of land use change. As a kind of spatial dynamics model, cellular automata (CA) has the characteristics of "bottom-up", which can retain more basis and details in dynamic system, and gradually occupies a dominant position in urban land change simulation. However, the traditional CA model based on grid is sensitive to the size and shape of the image element, and cannot truly express the objective geographical entity. At the same time, urban land use change has time dependence, which leads to the fact that the existing technology cannot accurately simulate urban land use change. SUMMARY

[0003] The present application aims to provide a method, medium and device for simulating urban land use change, which can take into account the time dependence of urban land use change and improve the accuracy of simulating urban land use change.

[0004] The technical solution adopted by the present application to solve its technical problems is to provide a method for simulating urban land use change, comprising the following steps: S1: obtaining a spatial driving factor dataset for model training and conversion probability calculation; S2: constructing an optimal time sequence model for calculating a conversion probability array of each plot; S3: inputting the spatial driving factor dataset into the optimal time sequence model to obtain the conversion probability array of each plot; S4: constructing a vector cellular automaton according to the spatial driving factor dataset to obtain a total development probability array of each plot; S5: merging the total development probability array of each plot with the conversion probability array of each plot to obtain a final conversion probability array of each plot, and performing a roulette algorithm operation on the final conversion probability array of each plot to obtain a converted land use class of each plot.

[0005] Further, the above-mentioned obtaining of the spatial driving factor dataset comprises the following steps: obtaining plot scale data, the plot scale data including the number of plots, the land use class of each plot, point of interest data, open road network data, and night light data, wherein the point of interest data includes dining, bus station, and life service site information; calculating the Euclidean distance between plots, such as the formula: d=sqrt[(x2-x1) 2 +(y2-y1) 2], where d is the Euclidean distance between the two plots, sqrt[] is the square root function, (x1, y1) is the centroid of the first plot, and (x2, y2) is the centroid of the second plot; the point of interest data, open road network data, and nighttime light data of the unsampled plots are estimated using spatial interpolation; the point density of each plot is calculated using the formula D = n / A, where D is the point density, n is the number of points in the calculation area, and A is the area of ​​the calculation area; the point of interest data, open road network data, and nighttime light data of each plot are normalized using the formula: X norm =(XX) min ) / (X max -X min Among them, X norm This is the normalized value, where X is the original value. min and X max These are the minimum and maximum values ​​of the calculated data, respectively, to obtain the spatial driving factor dataset. The spatial driving factor dataset includes land use data and spatial auxiliary variables. The land use data includes the number of plots, land use categories, leaf nodes of the decision tree propagation path of individual samples in the DT model, the total number of DTs trained by the RF model, the probability vector of each tree, the original output and calibration parameters of the SVM model, the number of individuals closest to the sample in the KNN model, the Euclidean distance between individuals, the classification label and all its values, and the sub-vectors of the LSTM and GRU models. The spatial auxiliary variables include digital elevation model, slope, distance to road, distance to railway, distance to highway, density of catering facilities, density of medical facilities, density of entertainment facilities, density of OSM road network, density of shopping facilities, density of factory facilities, density of public transportation facilities, density of park facilities, density of commercial facilities, and density of nighttime light.

[0006] Furthermore, the above-mentioned construction of the optimal time series model includes the following steps: S301: Obtain the spatial driving factor dataset, divide the spatial driving factor dataset into training set and test set according to time sequence, and set the training set:test set data ratio to 7:3; S302: Obtain the DT model, RF model, SVM model, KNN model, LSTM and GRU model, and adjust the training weight parameters of the DT model, RF model, SVM model, KNN model, LSTM and GRU model as shown in the formula: Where, α i Let i be the training weights for the i-th land use category. Let n be the number of samples for all land use categories, and K be the number of land use categories; jis the number of samples of the jth land use category; the land use category is: road or water body, unused land, industrial land, public service management land, commercial land, residential land; S303: input the training set and the test set into the DT model, the RF model, the SVM model, the KNN model, the LSTM and the GRU model respectively to obtain a pre-trained DT model, a pre-trained RF model, a pre-trained SVM model, a pre-trained KNN model, a pre-trained LSTM and a pre-trained GRU model; S304: input the training set into the pre-trained DT model, the pre-trained RF model, the pre-trained SVM model, the pre-trained KNN model, the pre-trained LSTM and the pre-trained GRU model respectively to obtain a conversion probability array of each plot under the DT model, the RF model, the SVM model, the KNN model, the LSTM model and the GRU model respectively; S305: construct a vector cellular automaton according to the training set to obtain a total development probability array of each plot; S306: merge the total development probability array of each plot with the conversion probability array of each plot under the DT model, the RF model, the SVM model, the KNN model, the LSTM model and the GRU model respectively and perform a roulette algorithm operation to obtain a conversion land use category of each plot under the DT model, the RF model, the SVM model, the KNN model, the LSTM model and the GRU model; S307: obtain an optimal time sequence model according to the conversion land use category of each plot under the DT model, the RF model, the SVM model, the KNN model, the LSTM model and the GRU model and in combination with the test set.

[0007] Further, the above step S304 comprises the following steps: inputting the training set into the pre-trained DT model to obtain a conversion probability array of each plot under the pre-trained DT model, as shown in the following formula: wherein p i[DT] is the conversion probability array of each plot under the pre-trained DT model, L is a leaf node satisfying a propagation path of an individual decision tree of a probability sample to be calculated, N L is the number of samples contained in the leaf node L, n i is the number of i-category samples contained in the leaf node L; inputting the training set into the pre-trained RF model, wherein the RF model is an ensemble learning method composed of multiple DTs, each DT model can output the probability of each category of sample data, and a conversion probability array of each plot under the pre-trained RF model is obtained, as shown in the following formula: wherein P [RF] is the conversion probability array of each plot under the pre-trained RF model, N is the total number of DT models of the pre-trained RF model, p n is the probability vector of the nth tree; inputting the test set into the pre-trained SVM model to obtain a conversion probability array of each plot under the pre-trained SVM model, as shown in the following formula: Where P(Y=1|Z) is the transition probability array of each land parcel under the pre-trained SVM model, Z represents the original output of the pre-trained SVM model, and A and B are calibration parameters estimated using the training dataset. The test set is input into the pre-trained KNN model, which predicts the category of the tested sample based on the categories of the K individuals closest to the land parcel, thus obtaining the transition probability array of each land parcel under the pre-trained KNN model, as shown in the formula: Among them, P i[KNN] w represents the transition probability array for each land parcel under the pre-trained KNN model. m These are the weighting coefficients. m represents the K closest individuals to the sample, and d m Let m be the Euclidean distance from individual m to the individual being predicted. Y For an individual with classification label Y, where y represents all possible values ​​for the classification label; the test set is input into a pre-trained LSTM model to obtain the transition probability array for each land parcel under the pre-trained LSTM model, as shown in the formula: Among them, P LSTM (Y=i|Z) is the transformation probability array for each plot under the pre-trained LSTM model, which is the probability that a sample belongs to class i, z i It is the i-th value in the score vector. This is the sum of the score vectors of all categories after performing an exponential operation, where C is the number of land parcel categories; the test set is input into the pre-trained GRU model to obtain the transition probability array of each land parcel under the pre-trained GRU model, as shown in the formula: Among them, P GRU (Y=i|Z) is the transition probability array for each plot under the pre-trained GRU model, which is the probability that a sample belongs to class i, z i It is the i-th value in the score vector. It is the sum of the score vectors of all categories after performing an exponential operation, where C is the number of land parcel categories.

[0008] Further, step S305 above includes the following steps: calculating the distance between the centroids of each plot before and after the conversion, determining the two plots with the smallest distance as the same plot before and after the change, and using them for plot matching operations; grouping the plots before the land use change according to land use categories, constructing a machine learning classification model for each group, using the spatial auxiliary variable data corresponding to each group of plots as attribute values, and the land category data after the land use change as label values, and calculating the development probability Pg of each plot through the machine learning classification model. i Using each plot as a unit, calculate the neighborhood effect of each plot, as shown in the formula: in, is the neighborhood effect of the jth plot on the ith plot at the tth time, d ij is the distance between the ith plot and the jth plot, S i represents the area of the ith plot, S j represents the area of the jth plot, d represents the field search radius, S max represents the maximum plot area, S min represents the minimum plot area; the conversion rule of the vector cellular automaton model is constructed as follows: wherein, is the overall development probability of the ith plot converted into the kth plot, the overall development probability of the ith plot converted into the kth plot is integrated into a vector, i.e., the overall development probability array P of each plot is obtained i[VGA] , is the development probability of the ith plot converted into the kth plot, is the neighborhood effect of the kth plot on the ith plot, is the development restriction factor of the ith plot converted into the kth plot, when the plot is a restricted development area, its development restriction factor is 0, and when the plot is a non-restricted development area, its development restriction factor is 1, RA is a random value.

[0009] Further, the above step S306 includes the following steps: merging the overall development probability array of each plot with the conversion probability array of each plot under the pre-trained DT model, the pre-trained RF model, the pre-trained SVM model, the pre-trained KNN model, the pre-trained LSTM model and the pre-trained GRU model, respectively, to sequentially obtain the first selection probability array of each plot, the second selection probability array of each plot, the third selection probability array of each plot, the fourth selection probability array of each plot, the fifth selection probability array of each plot, and the sixth selection probability array of each plot, as follows: i[DT-VCA] p i[VGA] = (P i[Dt] + P i[RF-VCA] ) / 2, p i[VGA] = (P i[RF] + P i[SVM-VCA] ) / 2, p i[VGA] = (P i[KNN-VCA] + P(Y = 1 | Z)) / 2, p i[VGA] = (P i[KNN] + P i[LSTM-VCA] (Y = i | Z)) / 2, p i[VGA] = (P LSTM + P i[GRU-VCA] (Y = i | Z)) / 2, wherein p i[VGA] is the first selection probability array, wherein p GRU is the second selection probability array, wherein p i[DT-VCA] is the third selection probability array, wherein p i[RF-VCA]is a first selection probability array, where p i[SVM-VCA] is a third selection probability array, where p i[KNN-VCA] is a fourth selection probability array, where p i[LSTM-VCA] is a fifth selection probability array, where p i[GRU-VCA] is a sixth selection probability array, where p i[VGA] is the overall development probability of the i-th plot converting to a k-class plot, p i[DT] is a plot conversion probability array under a pre-trained DT model, P i[RF] is a plot conversion probability array under a pre-trained RF model, P i[KNN] is a plot conversion probability array under a pre-trained SVM model, P LSTM is a plot conversion probability array under a pre-trained KNN model, P GRU is a plot conversion probability array under a pre-trained LSTM model, P e roulette algorithm is performed on the first, second, third, fourth, fifth and sixth selection probability arrays of the plots respectively to obtain the converted land use categories of the plots under the pre-trained DT model, the pre-trained RF model, the pre-trained SVM model, the pre-trained KNN model, the pre-trained LSTM model and the pre-trained GRU model.

[0010] Further, the step S307 includes the following steps: according to the converted land use categories of the plots under the pre-trained DT model, the pre-trained RF model, the pre-trained SVM model, the pre-trained KNN model, the pre-trained LSTM model and the pre-trained GRU model, the quality factor, the Kappa coefficient and the overall classification accuracy of the plot simulation results under the pre-trained DT model, the pre-trained RF model, the pre-trained SVM model, the pre-trained KNN model, the pre-trained LSTM model and the pre-trained GRU model are calculated respectively, as shown in the following formula: wherein, FoM is the quality factor, Kappa is the Kappa coefficient, OA is the overall classification accuracy, A is the area error of the simulation under the condition that the actual land use category is changed while the simulated land use category remains unchanged; B is the area error of the simulation under the condition that the actual land use category is changed while the simulated land use category is correctly converted; C is the area error of the simulation under the condition that the actual land use category is changed while the simulated land use category is incorrectly converted; D is the area error of the simulation under the condition that the actual land use category is changed while the simulated land use category is converted; p0 is the consistency quantity between the observed classifier prediction value and the true value, i.e., the overall classification accuracy, p eis the consistency expectation of random classification results; N is the total number of cell units in the study area; c is the total number of categories, ConfusionMatrix is a two-dimensional confusion matrix table (cxc size) composed of true labels and predicted labels, and ConfusionMatrix(i,i) is the diagonal element of the confusion matrix, that is, the individual whose predicted classification is correct; the quality factor, the Kappa coefficient and the overall classification accuracy of the land parcel simulation results under the pre-trained DT model, the pre-trained RF model, the pre-trained SVM model, the pre-trained KNN model, the pre-trained LSTM model and the pre-trained GRU model are investigated, and the pre-trained model with a quality factor greater than 0.2 and a Kappa coefficient greater than 0.6 is selected, and the pre-trained model with the highest overall classification accuracy is the optimal time series model.

[0011] Further, the above step S3 comprises the following steps: inputting the spatial driving factor dataset into the optimal time series model to obtain the conversion probability array of each land parcel, and if the optimal time series model is a pre-trained DT model, the conversion probability array of each land parcel is obtained according to the following formula: wherein p i[DT] is the conversion probability array of each land parcel under the pre-trained DT model, L is a leaf node satisfying the propagation path of the individual decision tree of the probability sample to be calculated, N L is the number of samples contained in the leaf node L, n i is the number of i-class samples contained in the leaf node L; if the optimal time series model is a pre-trained RF model, the conversion probability array of each land parcel is obtained according to the following formula: wherein P [RF] is the conversion probability array of each land parcel under the pre-trained RF model, N is the total number of DT models trained by the pre-trained RF model, and p n is the probability vector of the nth tree; if the optimal time series model is a pre-trained SVM model, the conversion probability array of each land parcel is obtained according to the following formula: wherein P(Y=1|Z) is the conversion probability array of each land parcel under the pre-trained SVM model, Z represents the original output of the pre-trained SVM model, and A and B are calibration parameters estimated by the training data set; if the optimal time series model is a pre-trained KNN model, the conversion probability array of each land parcel under the pre-trained KNN model is obtained by predicting the category of the sample to be measured through the category of the K nearest individuals to the land parcel, and the conversion probability array of each land parcel under the pre-trained KNN model is obtained according to the following formula: wherein P i[KNN] is the conversion probability array of each land parcel under the pre-trained KNN model, w m is the weight coefficient, m is the K nearest individuals to the sample, d m is the Euclidean distance between the individual m and the predicted individual, and m Yfor the individual with the classification label of Y, y is all values of the classification label; the optimal time series model if it is a pre-trained LSTM model, as the formula: wherein P LSTM (Y = i | Z) is the conversion probability array of each plot under the pre-trained LSTM model, that is, the probability that the sample belongs to category i, z i is the ith value in the score vector, is the sum of the score vectors after exponential operation on all categories, and C is the number of plot categories; the optimal time series model if it is a pre-trained GRU model, as the formula: wherein P GRU (Y = i | Z) is the conversion probability array of each plot under the pre-trained GRU model, that is, the probability that the sample belongs to category i, z i is the ith value in the score vector, is the sum of the score vectors after exponential operation on all categories, and C is the number of plot categories; step S4 comprises the following steps: calculating the distance between the centroids of each plot before and after conversion, determining the two plots with the smallest distance as the same plot before and after conversion, and used for plot matching operation; grouping the plots before land use change according to land use categories, constructing a machine learning classification model for each group, taking the spatial auxiliary variable data corresponding to each group of plots as attribute values, and taking the land use category data after land use change as label values, and calculating the development probability Pg i of each plot through the machine learning classification model; taking each plot as a cell unit, calculating the neighborhood effect of each plot, as the formula: wherein is the neighborhood effect of the jth plot on the ith plot at the tth moment, d ij is the distance between the ith plot and the jth plot, S i represents the area of the ith plot, S j represents the area of the jth plot, d represents the field search radius, Smax represents the maximum plot area, and Smin represents the minimum plot area; constructing a vector cellular automaton model conversion rule, as the formula: wherein is the total development probability of the ith plot converted into a k-class plot, and the total development probability of the ith plot converted into a k-class plot is integrated into a vector, that is, the total development probability array P i[VGA] of each plot, is the development probability of the ith plot converted into a k-class plot, is the neighborhood effect of the k-class plot on the ith plot, The development restriction factor for converting the i-th plot to a k-type plot is 0 when the plot is in a restricted development area and 1 when the plot is in an unrestricted development area. RA is a random value. Step S5 includes the following steps: summing the overall development probability array of each plot and the conversion probability array of each plot and taking the average to obtain the final conversion probability array of each plot; performing a roulette wheel algorithm operation on the final conversion probability array to obtain the final land use category of each plot.

[0012] The present invention also provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of the above-described urban land use change simulation method.

[0013] The present invention also provides a computer device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the steps of the above-described urban land use change simulation method.

[0014] The urban land use change simulation method, medium, and equipment provided by this invention have the following beneficial effects:

[0015] The urban land use change simulation method, medium, and equipment provided by this invention, targeting the strong spatiotemporal variation characteristics of land use change, fully utilize the advantages of machine learning and deep learning models in extracting high-level features and abstract information from data through multi-level nonlinear transformations and effectively mining time series information. Based on DT, RF, SVM, KNN, LSTM, and GRU models, an urban land use change simulation method is constructed. Historical state information of different land parcels is recorded, and the land use conversion probability under each model is calculated. Then, a vector CA model is constructed using cadastral parcels as basic unit cells to calculate the land use conversion probability, which is then merged with the conversion probabilities from the aforementioned machine learning and deep learning models. Finally, a roulette wheel algorithm is used to optimize the final land use conversion probability. This method can take into account the time dependence of urban land use change, more accurately simulate the dynamic development process of cities, urban expansion, and land use change, improve the accuracy of urban land use change simulation, effectively reveal its driving factors and the mutual conversion process between land parcels, and contribute to a deeper understanding of the time dependence mechanism of urban land use, providing a reference for the formulation of regional development policies. Attached Figure Description

[0016] The present invention will be further described below with reference to the accompanying drawings and embodiments. In the accompanying drawings:

[0017] Figure 1 This is a flowchart of the urban land use change simulation method provided by the present invention;

[0018] Figure 2is a block diagram of a computer device provided by the present application;

[0019] Figure 3 is a schematic diagram of a data preprocessing process provided by the present application;

[0020] Figure 4 is a schematic diagram of urban land use change simulation, model training and evaluation provided by the present application. DETAILED DESCRIPTION

[0021] In order to have a clearer understanding of the technical features, objectives and effects of the present application, the specific embodiments of the present application will be described in detail with reference to the accompanying drawings.

[0022] Figure 1 A schematic diagram of one embodiment of the urban land use change simulation method of the present embodiment is shown. In the following, the urban land use change simulation method will be briefly described by the present application embodiment:

[0023] The urban land use change simulation method provided by the present application comprises the following steps: S1: acquiring a spatial driving factor data set, which is used for model training and conversion probability calculation; S2: constructing an optimal time series model, which is used for calculating a conversion probability array of each land block; S3: inputting the spatial driving factor data set into the optimal time series model to obtain the conversion probability array of each land block; S4: constructing a vector cellular automaton according to the spatial driving factor data set to obtain a total development probability array of each land block; S5: merging the total development probability array of each land block with the conversion probability array of each land block to obtain a final conversion probability array of each land block, and performing a roulette algorithm operation on the final conversion probability array of each land block respectively to obtain a converted land use class of each land block.

[0024] Specifically, the above-mentioned acquiring of the spatial driving factor data set comprises the following steps: acquiring land block scale data, the land block scale data comprising land block quantity, land use class of each land block, interest point data, open road network data and night light data, wherein the interest point data comprises catering, bus station and life service site information; calculating the Euclidean distance between each land block, such as the formula: d=sqrt[(x2-x1) 2 +(y2-y1) 2 ], wherein d is the Euclidean distance between two land blocks, sqrt[] is a square root function, (x1, y1) is the centroid of the first land block, and (x2, y2) is the centroid of the second land block; estimating the interest point data, open road network data and night light data of unsampled land blocks, the estimation method being spatial interpolation; calculating the point density of each land block, such as the formula D=n / A, wherein D is the point density, n is the point number of the calculation area, and A is the area of the calculation area; normalizing the interest point data, open road network data and night light data of each land block, such as the formula: X norm =(X-Xmin ) / (X max -X min ) wherein X norm is the normalized value, X is the original value, X min and X max are the minimum and maximum values of the calculated data, respectively, to obtain a spatial driving factor dataset; wherein the spatial driving factor dataset includes land use data and spatial auxiliary variables, the land use data including the number of plots, land use categories, leaf nodes of DT model sample individual decision tree propagation paths, the total number of DTs trained by the RF model, the probability vector of each tree, the original output of the SVM model, calibration parameters, the number of individuals closest to the distance sample of the KNN model, the Euclidean distance between individuals, classification labels and all their values, and the directional vector of the LSTM and GRU models; and the spatial auxiliary variables including digital elevation model, slope, distance to road, distance to railway, distance to highway, density of catering facilities, density of medical facilities, density of entertainment facilities, density of OSM road network, density of shopping facilities, density of factory facilities, density of public transport facilities, density of park facilities, density of commercial facilities, and density of night light.

[0025] Specifically, the above constructing an optimal timing model includes the following steps: S301: obtaining a spatial driving factor dataset, dividing the spatial driving factor dataset into a training set and a test set according to time sequence, and setting the training set:test set data ratio to 7:3; S302: obtaining a DT model, an RF model, an SVM model, a KNN model, an LSTM, and a GRU model, and adjusting the training weight parameters of the DT model, the RF model, the SVM model, the KNN model, the LSTM, and the GRU model, such as the formula: wherein α i is the training weight of the i-th land use category, is the sample number of all land use categories, and K is the number of land use categories; n jis the number of samples of the jth land use category; the land use category is: road or water body, unused land, industrial land, public service management land, commercial land, residential land; S303: input the training set and the test set into the DT model, the RF model, the SVM model, the KNN model, the LSTM and the GRU model respectively to obtain a pre-trained DT model, a pre-trained RF model, a pre-trained SVM model, a pre-trained KNN model, a pre-trained LSTM and a pre-trained GRU model; S304: input the training set into the pre-trained DT model, the pre-trained RF model, the pre-trained SVM model, the pre-trained KNN model, the pre-trained LSTM and the pre-trained GRU model respectively to obtain a conversion probability array of each land parcel under the DT model, the RF model, the SVM model, the KNN model, the LSTM model and the GRU model respectively; S305: construct a vector cellular automaton according to the training set to obtain a total development probability array of each land parcel; S306: merge the total development probability array of each land parcel with the conversion probability array of each land parcel under the DT model, the RF model, the SVM model, the KNN model, the LSTM model and the GRU model respectively and perform a roulette algorithm operation to obtain a conversion land use category of each land parcel under the DT model, the RF model, the SVM model, the KNN model, the LSTM model and the GRU model; S307: obtain an optimal time sequence model according to the conversion land use category of each land parcel under the DT model, the RF model, the SVM model, the KNN model, the LSTM model and the GRU model and in combination with the test set.

[0026] Specifically, the above step S304 includes the following steps: inputting the training set into the pre-trained DT model to obtain a conversion probability array of each land parcel under the pre-trained DT model, as shown in the following formula: wherein p i[DT] is the conversion probability array of each land parcel under the pre-trained DT model, L is a leaf node satisfying a propagation path of an individual decision tree of a probability sample to be calculated, N L is the number of samples contained in the leaf node L, n i is the number of i-category samples contained in the leaf node L; inputting the training set into the pre-trained RF model, wherein the RF model is an ensemble learning method composed of multiple DTs, each DT model can output the probability of each category of sample data, and a conversion probability array of each land parcel under the pre-trained RF model is obtained, as shown in the following formula: wherein P [RF] is the conversion probability array of each land parcel under the pre-trained RF model, N is the total number of DT models of the pre-trained RF model, p n is the probability vector of the nth tree; inputting the test set into the pre-trained SVM model to obtain a conversion probability array of each land parcel under the pre-trained SVM model, as shown in the following formula: Where P(Y=1|Z) is the transition probability array of each land parcel under the pre-trained SVM model, Z represents the original output of the pre-trained SVM model, and A and B are calibration parameters estimated using the training dataset. The test set is input into the pre-trained KNN model, which predicts the category of the tested sample based on the categories of the K individuals closest to the land parcel, thus obtaining the transition probability array of each land parcel under the pre-trained KNN model, as shown in the formula: Among them, P i[KNN] w represents the transition probability array for each land parcel under the pre-trained KNN model. m These are the weighting coefficients. m represents the K closest individuals to the sample, and d m Let m be the Euclidean distance from individual m to the individual being predicted. Y For an individual with classification label Y, where y represents all possible values ​​for the classification label; the test set is input into a pre-trained LSTM model to obtain the transition probability array for each land parcel under the pre-trained LSTM model, as shown in the formula: Among them, P LSTM (Y=i|Z) is the transformation probability array for each plot under the pre-trained LSTM model, which is the probability that a sample belongs to class i, z i It is the i-th value in the score vector. This is the sum of the score vectors of all categories after performing an exponential operation, where C is the number of land parcel categories; the test set is input into the pre-trained GRU model to obtain the transition probability array of each land parcel under the pre-trained GRU model, as shown in the formula: Among them, P GRU (Y=i|Z) is the transition probability array for each plot under the pre-trained GRU model, which is the probability that a sample belongs to class i, z i It is the i-th value in the score vector. It is the sum of the score vectors of all categories after performing an exponential operation, where C is the number of land parcel categories.

[0027] Specifically, step S305 above includes the following steps: calculating the distance between the centroids of each plot before and after the conversion, determining the two plots with the smallest distance as the same plot before and after the change, for plot matching operations; grouping the plots before the land use change according to land use categories, constructing a machine learning classification model for each group, using the spatial auxiliary variable data corresponding to each group of plots as attribute values, and the land category data after the land use change as label values, and calculating the development probability Pg of each plot through the machine learning classification model. i Using each plot as a unit, calculate the neighborhood effect of each plot, as shown in the formula: in, is the neighborhood effect of the jth plot on the ith plot at the tth time, d ij is the distance between the ith plot and the jth plot, S i represents the area of the ith plot, S j represents the area of the jth plot, d represents the field search radius, S max represents the maximum plot area, S min represents the minimum plot area; the conversion rule of the vector cellular automaton model is constructed as formula: wherein, is the overall development probability of the ith plot converted into the kth plot, the overall development probability of the ith plot converted into the kth plot is integrated into a vector, that is, the overall development probability array P of each plot is obtained i[VGA] , is the development probability of the ith plot converted into the kth plot, is the neighborhood effect of the kth plot on the ith plot, is the development restriction factor of the ith plot converted into the kth plot, when the plot is a restricted development area, the development restriction factor is 0, and when the plot is a non-restricted development area, the development restriction factor is 1, and RA is a random value.

[0028] Specifically, the above step S306 includes the following steps: merging the overall development probability array of each plot with the conversion probability array of each plot under the pre-trained DT model, the pre-trained RF model, the pre-trained SVM model, the pre-trained KNN model, the pre-trained LSTM model and the pre-trained GRU model, respectively, to obtain the first selection probability array of each plot, the second selection probability array of each plot, the third selection probability array of each plot, the fourth selection probability array of each plot, the fifth selection probability array of each plot, and the sixth selection probability array of each plot, respectively, as formula: i[DT-VCA] p i[VGA] = (P i[DT] + p i[RF-VCA] ) / 2, p i[VGA] = (P i[RF] + P i[SVM-VCA] ) / 2, p i[VGA] = (P i[KNN-VCA] + P(Y = 1|Z)) / 2, p i[VGA] = (P i[KNN] + P i[LSTM-VCA] ) / 2, p i[VGA] = (P LSTM + P i[GRU-VCA] (Y = i|Z)) / 2, p i[VGA] = (P GRU + P i[DT-VCA] (Y = i|Z)) / 2, wherein p i[RF-VCA]is a first selection probability array, where p i[SVM-VCA] is a third selection probability array, where p i[KNN-VCA] is a fourth selection probability array, where p i[LSTM-VCA] is a fifth selection probability array, where p i[GRU-VCA] is a sixth selection probability array, where p i[VGA] is the overall development probability of the ith plot converting to a kth land use type, p i[DT] is a first plot conversion probability array under a pre-trained DT model, P i[RF] is a second plot conversion probability array under a pre-trained RF model, P i[KNN] is a third plot conversion probability array under a pre-trained SVM model, P LSTM is a fourth plot conversion probability array under a pre-trained KNN model, P GRU is a fifth plot conversion probability array under a pre-trained LSTM model, P e is a sixth plot conversion probability array under a pre-trained GRU model; the first, second, third, fourth, fifth and sixth plot selection probability arrays are respectively subjected to a roulette algorithm operation to obtain a conversion land use type of each plot under the pre-trained DT model, the pre-trained RF model, the pre-trained SVM model, the pre-trained KNN model, the pre-trained LSTM model and the pre-trained GRU model.

[0029] Specifically, the step S307 includes the following steps: according to the conversion land use type of each plot under the pre-trained DT model, the pre-trained RF model, the pre-trained SVM model, the pre-trained KNN model, the pre-trained LSTM model and the pre-trained GRU model, respectively calculating a figure of merit, a Kappa coefficient and a classification overall accuracy of the plot simulation result under the pre-trained DT model, the pre-trained RF model, the pre-trained SVM model, the pre-trained KNN model, the pre-trained LSTM model and the pre-trained GRU model, as shown in the following formula: wherein, FoM is the figure of merit, Kappa is the Kappa coefficient, OA is the classification overall accuracy, A is an area error of simulation under the condition that the actual land use type is changed and the simulated land use type remains unchanged; B is an area error of simulation under the condition that the actual land use type is changed and the simulated land use type is correctly converted; C is an area error of simulation under the condition that the actual land use type is changed and the simulated land use type is incorrectly converted; D is an area error of simulation under the condition that the actual land use type is changed and the simulated land use type is converted; p0 is the consistency quantity between the observed classifier prediction value and the true value, i.e., the classification overall accuracy, p eis the consistency expectation of random classification results; N is the total number of cell units in the study area; c is the total number of categories, ConfusionMatrix is a two-dimensional confusion matrix table (cxc size) composed of true labels and predicted labels, and ConfusionMatrix(i,i) is the diagonal element of the confusion matrix, that is, the individual whose predicted classification is correct; the quality factor, the Kappa coefficient and the overall classification accuracy of the land parcel simulation results under the pre-trained DT model, the pre-trained RF model, the pre-trained SVM model, the pre-trained KNN model, the pre-trained LSTM model and the pre-trained GRU model are investigated, and the pre-trained model with a quality factor greater than 0.2 and a Kappa coefficient greater than 0.6 is selected, and the pre-trained model with the highest overall classification accuracy is the optimal time series model.

[0030] Specifically, the above step S3 comprises the following steps: inputting the spatial driving factor dataset into the optimal time series model to obtain the conversion probability array of each land parcel, wherein if the optimal time series model is a pre-trained DT model, the conversion probability array of each land parcel is obtained according to the following formula: wherein p i[DT] is the conversion probability array of each land parcel under the pre-trained DT model, L is a leaf node satisfying the propagation path of the individual decision tree of the probability sample to be calculated, N L is the number of samples contained in the leaf node L, n i is the number of i-class samples contained in the leaf node L; if the optimal time series model is a pre-trained RF model, the conversion probability array of each land parcel is obtained according to the following formula: wherein P [RF] is the conversion probability array of each land parcel under the pre-trained RF model, N is the total number of DT models trained by the pre-trained RF model, and p n is the probability vector of the nth tree; if the optimal time series model is a pre-trained SVM model, the conversion probability array of each land parcel is obtained according to the following formula: wherein P(Y=1|Z) is the conversion probability array of each land parcel under the pre-trained SVM model, Z represents the original output of the pre-trained SVM model, and A and B are calibration parameters estimated by the training data set; if the optimal time series model is a pre-trained KNN model, the conversion probability array of each land parcel under the pre-trained KNN model is obtained by predicting the category of the sample to be measured through the category of the K nearest individuals to the land parcel, and the conversion probability array of each land parcel under the pre-trained KNN model is obtained according to the following formula: wherein P i[KNN] is the conversion probability array of each land parcel under the pre-trained KNN model, w m is the weight coefficient, m is the K nearest individuals to the sample, d m is the Euclidean distance between the individual m and the predicted individual, and m Yfor the individual with the classification label of Y, y is all values of the classification label; the optimal time series model if it is a pre-trained LSTM model, as the formula: wherein P LSTM (Y = i | Z) is the conversion probability array of each plot under the pre-trained LSTM model, that is, the probability that the sample belongs to category i, z i is the ith value in the score vector, is the sum of the score vectors after exponential operation on all categories, and C is the number of plot categories; the optimal time series model if it is a pre-trained GRU model, as the formula: wherein P GRV (Y = i | Z) is the conversion probability array of each plot under the pre-trained GRU model, that is, the probability that the sample belongs to category i, z i is the ith value in the score vector, is the sum of the score vectors after exponential operation on all categories, and C is the number of plot categories; step S4 includes the following steps: calculating the distance between the centroids of each plot before and after conversion, determining the two plots with the smallest distance as the same plot before and after conversion, and used for plot matching operation; grouping the plots before land use change according to land use categories, constructing a machine learning classification model for each group, taking the spatial auxiliary variable data corresponding to each group of plots as attribute values, and taking the land use category data after land use change as label values, and calculating the development probability Pg i of each plot through the machine learning classification model; taking each plot as a cell unit, calculating the neighborhood effect of each plot, as the formula: wherein is the neighborhood effect of the jth plot on the ith plot at the tth moment, d ij is the distance between the ith plot and the jth plot, S i represents the area of the ith plot, S j represents the area of the jth plot, d represents the field search radius, Smax represents the maximum plot area, and Smin represents the minimum plot area; constructing a vector cellular automaton model conversion rule, as the formula: wherein is the total development probability of the ith plot converted into a k-class plot, and the total development probability of the ith plot converted into a k-class plot is integrated into a vector, that is, the total development probability array P i[VGA] of each plot, is the development probability of the ith plot converted into a k-class plot, is the neighborhood effect of the k-class plot on the ith plot, The development restriction factor of the ith plot converted into the kth plot is 0 when the plot is a restricted development area, and is 1 when the plot is a non-restricted development area, and RA is a random value; step S5 comprises the following steps: summing and averaging the total development probability array of the plots and the conversion probability array of the plots to obtain a final conversion probability array of the plots; and performing a roulette algorithm operation on the final conversion probability array to obtain a final conversion land use category of the plots.

[0031] The application further provides a computer readable storage medium, which stores a computer program, and the computer program is executed by a processor to realize the steps of the urban land use change simulation method.

[0032] The application further provides a computer device, which comprises a memory, a processor and a computer program stored in the memory and capable of running on the processor, and the processor realizes the steps of the urban land use change simulation method when executing the program.

[0033] The application further provides a computer readable storage medium, which stores a computer program, and the computer program is executed by a processor to realize the steps of the urban land use change simulation method. The storage medium can be a disc, an optical disc, a read-only memory (ROM), a random access memory (RAM), a flash memory, a hard disk drive (HDD) or a solid-state drive (SSD), etc.; and the storage medium can also comprise a combination of the above-mentioned memories.

[0034] The application further provides a computer device, which comprises a memory, a processor and a computer program stored in the memory and capable of running on the processor, and the processor realizes the steps of the urban land use change simulation method when executing the program. For example, Figure 2As shown, the computer device can include at least one processor 151, such as a CPU (Central Processing Unit), at least one communication interface 153, a memory 154, at least one communication bus 152. Among them, the communication bus 152 is used to realize the connection communication between these components. Among them, the communication interface 153 can include a display screen (Display), a keyboard (Keyboard), and the optional communication interface 153 can also include a standard wired interface, a wireless interface. The memory 154 can be a high-speed RAM memory (Random Access Memory), or a non-volatile memory, such as at least one disk memory. The memory 154 can also be at least one storage device located away from the aforementioned processor 151. Among them, the memory 154 stores an application program, and the processor 151 calls the program code stored in the memory 154 for executing any of the above method steps. Among them, the communication bus 152 can be a peripheral component interconnect (peripheral component interconnect, abbreviated as PCI) bus or an extended industry standard architecture (extended industry standard architecture, abbreviated as EISA) bus, etc. The communication bus 152 can be divided into an address bus, a data bus, a control bus, etc. For convenience of representation, Figure 2The bus 153 illustrated in FIG. 1 is only one line, but does not mean that there is only one bus or only one type of bus. Among them, the memory 154 can include volatile memory such as random-access memory (RAM); the memory can also include non-volatile memory such as flash memory, a hard disk drive (HDD) or a solid-state drive (SSD); the memory 154 can also include a combination of the above-mentioned types of memory. Among them, the processor 151 can be a central processing unit (CPU), a network processor (NP), or a combination of CPU and NP. Among them, the processor 151 can further include a hardware chip. The above-mentioned hardware chip can be an application-specific integrated circuit (ASIC), a programmable logic device (PLD), or a combination thereof. The above-mentioned PLD can be a complex programmable logic device (CPLD), a field-programmable gate array (FPGA), a generic array logic (GAL), or any combination thereof. Alternatively, the memory 154 is also used to store program instructions. The processor 151 can invoke the program instructions to implement the urban land use change simulation method as in the embodiment.

[0035] In some embodiments, the above-mentioned urban land use change simulation method can be implemented in the following manner.

[0036] The method for simulating urban land use change comprises the following steps: (1) time series model construction: this step aims to use various machine learning models to mine time series factors and their probabilities in the process of land use change based on preprocessed land use change time series data; (2) land use development probability merging and land use simulation: this step merges the prediction probabilities of the time series model with the prediction probabilities calculated by the VCA model for different types of land, and selects the final development category through the roulette algorithm; (3) urban land use prediction under different scenarios: this step designs two future scenarios of disordered development and ecological protection to simulate long-time series land use development and analyze the time series characteristics of land use; (4) model accuracy verification: this step evaluates the accuracy based on the simulated urban land use change results and the actual land use.

[0037] In some embodiments, in step (1) of the above-mentioned method for simulating urban land use change, the time series prediction model based on various machine learning models automatically extracts features from time series data and performs end-to-end training and prediction using traditional machine learning models and deep learning models. In this embodiment, the land use time series data of Shenzhen City from 2009 to 2013 is used as training data, the land use data of 2014 is used as accuracy verification data, and the training set: test set data ratio is set to 7:3. However, the data has the problem of unbalanced multi-classification samples, so this embodiment adopts the method of adjusting the training weight parameter to ensure that the model is relatively balanced in training on different classifications. The training weight parameter is calculated as follows:

[0038]

[0039] In the formula: α i represents the training weight of the i-th land use category; represents the sample size of all land use categories in the study area; n i represents the sample size of the i-th land use category. The calculated weight coefficients of each category are shown in Table 1:

[0040] Table 1 Weight coefficients of each land use category

[0041]

[0042] To fully compare the prediction performance of each time series prediction model in this multi-classification task, the embodiment selects four machine learning models, including decision tree (DT), support vector machine (SVM), random forest (RF), and K-nearest neighbor (KNN), and two deep learning models, including long short-term memory (LSTM) and gated recurrent units (GRU), aiming to mine the time series factors and their probabilities in the process of land use conversion. The embodiment calculates the land use conversion probabilities under different machine learning and deep learning models. Since the basic mathematical principles of different machine learning methods are different, the above six classifiers calculate the prediction probabilities of each category in different ways. The calculation formula of the prediction probability p i of category i in the DT model is as follows:

[0043]

[0044] where L is the leaf node that meets the propagation path of the individual decision tree of the sample to be calculated, N L is the number of samples contained in the leaf node L, and n i is the number of category i samples contained in the leaf node L. The RF model is an ensemble learning method composed of multiple DTs, each of which can output the probability of each category of sample data. The calculation formula of the prediction probability P of the RF model is as follows:

[0045]

[0046] where N is the total number of DTs trained by RF, and n is the probability vector of the nth tree. SVM is originally a binary classification model, which can be used for multi-classification problems through the one-vs-the-rest (OvR) method and uses Platt scaling to calculate the probability of each category. This can map the original output (usually a real number) of the classifier to the probability estimate of the category, so that these probabilities are more consistent with the true category probability distribution. The probability calculation formula after Platt scaling is as follows:

[0047]

[0048] where P(Y=1|Z) represents the probability that the sample belongs to the positive category, Z represents the original output of SVM, and A and B are calibration parameters estimated by the training data set. KNN predicts the category of the measured sample by the categories of the K nearest individuals, and the probability calculation formula is as follows:

[0049]

[0050] where w is the weight coefficient m The calculation formula is:

[0051]

[0052] In the formula, m is the nearest K individuals to the sample, d m is the Euclidean distance from individual m to the predicted individual, m Y is the individual with classification label Y, and y is all values of the classification label. The LSTM and GRU use the Softmax function to calculate the probability distribution of the category, and the mathematical calculation formula is as follows:

[0053]

[0054] In the formula, P(Y=i|Z) represents the probability that the sample belongs to category i, z i is the i-th value in the score vector logits, is the sum of the logits after exponential operation on all categories.

[0055] In some embodiments, in step (2) of the above urban land use change simulation method, based on the land use conversion probability calculated by the above model, the final development probability of each plot is obtained by combining the conversion probability calculated by the time series prediction model and the traditional VCA model. The calculation formula of the VCA model development probability is as follows:

[0056]

[0057] In the formula: is the overall development probability of the i-th plot converted to the k-th plot, is the development probability of each plot, is the neighborhood effect of the k-th plot on the i-th plot, is the development restriction factor, and RA is a random value. Finally, the roulette algorithm is used to determine the final category of land use development change, and the calculation method is:

[0058]

[0059] In the formula: q i is the sum of the land use category conversion probabilities after combining the time series model and the VCA model; p iSelecting probabilities for roulette. In some embodiments, in step (3) of the above-mentioned urban land use change simulation method, in the context of global environmental change and the promotion of sustainable development, urban land use change prediction is of great significance for guiding land resource management, urban planning policy making and sustainable development. To effectively predict land use change in Shenzhen, two development scenarios were designed: (1) an unordered development scenario (O1); (2) an ecological protection scenario (O2). In scenario O1, following the historical development trend, mutual conversion between various types of land is allowed; in scenario O2, all water bodies and roads within the ecological red line area are prohibited from being developed into four types of urban land: industrial land, public service and management land, commercial land, and residential land. Optimal utilization and protection of land resources can effectively balance economic development and environmental protection, and more precise land use and ecological protection policies can be developed.

[0060] In some embodiments, in step (4) of the above-mentioned urban land use change simulation method, to comprehensively evaluate the simulation effect of the constructed model, this embodiment selects evaluation indexes such as Figure-of-Merit (FoM), Kappa coefficient, and Overall Accuracy (OA) to perform multi-angle model evaluation on the model, wherein the calculation formula of FoM is as follows:

[0061]

[0062] In the formula, A represents the area error of actual land use type conversion that remains unchanged in the simulation result; B represents the area error of actual conversion that occurs correct type conversion in the simulation result; C represents the area error of actual conversion that occurs incorrect type conversion in the simulation; D represents the area error of actual non-conversion that occurs conversion in the simulation; and N represents the total number of cell units in the study area.

[0063] In some embodiments, the above-mentioned urban land use change simulation method comprises the following steps: the land plots are divided into 5 groups, the mean value of the spatial auxiliary variable corresponding to each group of land plots is defined as X, and the land use type after the change of land use is defined as Y, and a Y=f(X) model is constructed. Then, the probability of land conversion to each land use type Y in the model is mined by a machine learning algorithm as the overall development probability Pgi of the land plot. Finally, the overall development probability of the land use type without conversion is set to 0, so as to not only consider the true situation that a large number of land use types do not change, but also mine the potential of land use change of each land plot. Three machine learning algorithms (logistic regression model, neural network model and random forest model) are used to calculate the overall development probability of each land plot. Among them, the logistic regression (Logistic regression, LR) model is a multivariate analysis model, which allows to build a multiple regression relationship between independent and non-independent variables and make predictions based on variable values. The LR model can not only reveal the main factors affecting urban expansion, but also determine to what extent they promote or inhibit urban expansion, and has been widely used in urban land use simulation. The neural network (Neural network, NN) model is a widely used modeling technique with adaptive, self-organizing and self-learning capabilities, which can model land use change trends. The NN model is composed of layers and neurons that simulate the structure of the human brain, and has good fitting effect in nonlinear mapping. The NN model usually contains three types of layers: input layer, hidden layer and output layer, wherein the hidden layer is a set of neurons that apply activation functions, providing an intermediate layer between the input layer and the output layer. At the same time, the NN model can be used for model performance evaluation through k-fold cross-validation (k-fold cross-validation). For each round of cross-validation, a specified number of random restarts can be performed, so that the best network is selected after training. The random forest (Random forest, RF) model is an ensemble model that can be used to mine complex conversion rules of the CA model, which uses bagging to build many individual decision trees to achieve prediction or classification. These decision trees are built by randomly selecting a subset of the training set. When training the RF algorithm, the Bootstrap sampling method is used to obtain a sample dataset X containing ntree items (i=1, 2..., ntree) instead of the original dataset X. Other data that do not appear in the X; dataset are called "out-of-bag" (Out-of-bag, OOB) data. The RF model can use the OOB dataset for OOB prediction, estimate the OOB error and evaluate the importance of spatial variables, and determine the classification result according to the voting result of these decision trees.

[0064] The basic principle of roulette: for example, 【0.45 (first type of land), 0.25 (second type of land), 0.3 (third type of land)】 or 【0.45 (residential land), 0.25 (factory), 0.3 (park)】, the first 0.45 is not necessarily the final conversion of land type, that is, the roulette is operated on this 【0.45, 0.25, 0.3】, and the final selection of which type is the final conversion type is determined, the probability of the largest probability is generally selected, but it is also possible to select those with small probability, which needs to be considered in combination with random numbers. First, give a random number between 0 and 1, such as 0.5, 0-0.45 is the first type, 0.45-0.7 is the second type, and 0.7-1 is the third type, so the result is the second type. The specific operation is to start from 0 and add each probability in turn, 0+0.45, 0.45<0.5, so it is not yet, 0.45+0.25=0.7>0.5, which is located in the interval of 0.45-0.7. The second type can be selected. Give a label array, 【0-0+0.45】, 【0.45-0.45+0.25】, 【0.7, 0.7+0.3】. This is a probability array. It is similar to the force of your roulette, which is randomly provided by the system, and where it turns to is which one. If the random number is 0.2, it is the first type, if the random number is 0.5, it is the second; if 0.8, it is the third.

[0065] In some embodiments, the above urban land use change simulation method comprises the following steps:

[0066] (1) Data preprocessing: as shown in Figure 3 , the data sources are Baidu interest point data (Point Of Interest, POI), open road network data (Open Street Map, OSM), and night light data; the number of interest points, road network density, and night light data at the plot scale are extracted, wherein the interest point data mainly includes restaurant, bus station, and life service information, the data is processed by calculating the Euclidean distance, spatial interpolation, and point density analysis, and all variables are normalized to 0-1, and the spatial resolution is set to 30 meters. Calculate the Euclidean distance: the Euclidean distance between two points in two-dimensional or three-dimensional space is the straight line distance between the two points. In two-dimensional space, the Euclidean distance between two points (x l ,y1) and (x2,y2) is calculated by the following formula: d=sqrt[(x2-x1) 2 +(y2-y1) 2]; spatial interpolation: spatial interpolation is a technique used to predict the value of unsampled locations in geographic space; point density analysis: point density analysis is a method used to describe the distribution density of points in geographic space. A "moving window" is usually used to move on the map, calculate the number of points in the window, and then divide the number of points by the area of the window to get the point density of the location. The formula for calculating point density is: D = n / A; normalization: normalization is a method of converting data to a uniform scale, usually converting data to between 0 and 1. The formula is: X norm = (X-X min) / (X max-X min) where X norm is the normalized value, X is the original value, X min and X max are the minimum and maximum values of the data respectively; finally, the spatial driving factor data set is obtained. norm min max min min max

[0067] (2) Model training process: based on various machine learning and deep learning models, land use transition probabilities are calculated with the aid of land use time series data and spatial auxiliary variable data. The land use transition probabilities calculated by the traditional vector cellular automaton model are coupled with the probabilities calculated by the machine learning and deep learning models to obtain the overall land use transition probability. The roulette algorithm is used to determine the land use conversion category to obtain the simulation result of urban land use change. Based on different land use policy settings, future development scenarios are simulated to further simulate and predict land use, as shown in the process of Figure 4

[0068] The embodiments of the application are described above in combination with the drawings, but the application is not limited to the specific embodiments described above, and the specific embodiments described above are only illustrative and not limiting. Those skilled in the art can make many forms under the inspiration of the application without departing from the purpose of the application and the scope protected by the claims, and these are all within the protection of the application.​​​​​​​

Claims

1. A method of urban land use change simulation, characterized by, The method comprises the following steps: S1: obtaining a spatial driving factor dataset for model training and conversion probability calculation; S2: constructing an optimal time series model for calculating a conversion probability array of each plot; specifically comprising the following steps: S301: obtaining a spatial driving factor dataset, dividing the spatial driving factor dataset into a training set and a test set according to time sequence, and setting the training set:test set data ratio to 7:3; S302: Obtain the DT model, RF model, SVM model, KNN model, LSTM model, and GRU model, and adjust the training weight parameters of the DT model, RF model, SVM model, KNN model, LSTM model, and GRU model as shown in the formula: ,in, For the first Training weights for land use categories, K represents the number of samples across all land use categories, where K is the number of land use categories. For the first The number of samples for each land use category; the land use categories are: roads or water bodies, unused land, industrial land, public service management land, commercial land, and residential land; S303: inputting the training set and the test set into the DT model, the RF model, the SVM model, the KNN model, the LSTM and the GRU model respectively to obtain a pre-trained DT model, a pre-trained RF model, a pre-trained SVM model, a pre-trained KNN model, a pre-trained LSTM and a pre-trained GRU model; S304: inputting the training set into the pre-trained DT model, the pre-trained RF model, the pre-trained SVM model, the pre-trained KNN model, the pre-trained LSTM and the pre-trained GRU model respectively to obtain a conversion probability array of each plot under the DT model, the RF model, the SVM model, the KNN model, the LSTM model and the GRU model; S305: constructing a vector cellular automaton according to the training set to obtain a total development probability array of each plot; S306: merging the total development probability array of each plot with the conversion probability array of each plot under the DT model, the RF model, the SVM model, the KNN model, the LSTM model and the GRU model respectively and performing a roulette algorithm operation to obtain a converted land use category of each plot under the DT model, the RF model, the SVM model, the KNN model, the LSTM model and the GRU model; S307: obtaining an optimal time series model according to the converted land use category of each plot under the DT model, the RF model, the SVM model, the KNN model, the LSTM model and the GRU model in combination with the test set; S3: inputting the spatial driving factor dataset into the optimal time series model to obtain a conversion probability array of each plot; S4: constructing a vector cellular automaton according to the spatial driving factor dataset to obtain a total development probability array of each plot; S5: merging the total development probability array of each plot with the conversion probability array of each plot to obtain a final conversion probability array of each plot, and performing a roulette algorithm operation on the final conversion probability array of each plot to obtain a converted land use category of each plot.

2. The urban land use change simulation method according to claim 1, characterized by, The acquiring the spatial driving factor dataset comprises the following steps: acquiring patch scale data, the patch scale data comprising patch quantity, land use category of each patch, point of interest data, open road network data, and night light data, the point of interest data comprising catering, bus station, and life service site information; calculating Euclidean distance between patches, such as formula: d = sqrt[(x2-x1) 2 +(y2-y1) 2 ], wherein d is the Euclidean distance between two patches, sqrt[] is a square root function, (x1, y1) is the centroid of the first patch, and (x2, y2) is the centroid of the second patch; estimating point of interest data, open road network data, and night light data of unsampled patches, the estimation method being spatial interpolation; calculating point density of each patch, such as formula D = n / A, wherein D is point density, n is the number of points in the calculation area, and A is the area of the calculation area; normalizing point of interest data, open road network data, and night light data of each patch, such as formula: X norm =(X-X min ) / (X max -X min ) wherein X norm is the normalized value, X is the original value, X min and X max are the minimum value and the maximum value of the calculated data, respectively, to obtain the spatial driving factor dataset; the spatial driving factor dataset comprising land use data and spatial auxiliary variables, the land use data comprising patch quantity, land use category, leaf node of DT model sample individual decision tree propagation path, total number of DT trained by RF model, probability vector of each tree, original output of SVM model, calibration parameter, individual number of KNN model distance sample nearest individual, Euclidean distance between individuals, classification label and all values thereof, and split vector of LSTM and GRU model; the spatial auxiliary variables comprising digital elevation model, slope, distance to road, distance to railway, distance to highway, catering facility density, medical facility density, entertainment facility density, OSM road network density, shopping facility density, factory facility density, bus facility density, park facility density, commercial facility density, and night light density.

3. The urban land change simulation method according to claim 1, characterized by, Step S304 includes the following steps: inputting the training set into the pre-trained DT model to obtain a conversion probability array of each plot under the pre-trained DT model, as shown in the following formula: wherein is the conversion probability array of each plot under the pre-trained DT model, L is a leaf node satisfying the propagation path of the individual decision tree of the probability sample to be calculated, is the number of samples contained in the leaf node L, is the number of i-class samples contained in the leaf node L; inputting the training set into the pre-trained RF model, the RF model being an ensemble learning method composed of multiple DTs, each DT model being capable of outputting the probability of each class of sample data, to obtain a conversion probability array of each plot under the pre-trained RF model, as shown in the following formula: wherein is the conversion probability array of each plot under the pre-trained RF model, N is the total number of DT models trained by the pre-trained RF model, is the probability vector of the nth tree; inputting the test set into the pre-trained SVM model to obtain a conversion probability array of each plot under the pre-trained SVM model, as shown in the following formula: wherein, is the conversion probability array of each plot under the pre-trained SVM model, represents the original output of the pre-trained SVM model, and is a calibration parameter, which is estimated by the training data set; inputting the test set into the pre-trained KNN model, the pre-trained KNN model predicting the class of the sample to be measured by the class of the K nearest individuals of the plot to obtain a conversion probability array of each plot under the pre-trained KNN model, as shown in the following formula: wherein, is the conversion probability array of each plot under the pre-trained KNN model, is a weight coefficient, , is the K nearest individuals of the sample, is the Euclidean distance from the individual to the individual to be predicted, is the individual with the classification label Y, y is all values of the classification label; inputting the test set into the pre-trained LSTM model to obtain a conversion probability array of each plot under the pre-trained LSTM model, as shown in the following formula: wherein, is the conversion probability array of each plot under the pre-trained LSTM model, that is, the probability of the sample belonging to the class , is the ith value in the score vector, ​This is the sum of the score vectors of all categories after performing an exponential operation, where C is the number of land parcel categories; the test set is input into the pre-trained GRU model to obtain the transition probability array of each land parcel under the pre-trained GRU model, as shown in the formula: ,in, This is the transition probability array for each plot under the pre-trained GRU model, which represents the category to which the sample belongs. The probability, It is the first in the score vector One value, It is the sum of the score vectors of all categories after performing an exponential operation, where C is the number of land parcel categories.

4. The urban land use change simulation method according to claim 3, characterized by, Step S305 includes the following steps: calculating the distance between the centroids of each plot before and after the conversion, identifying the two plots with the smallest distance as the same plot before and after the change, and using them for plot matching operations; grouping the plots before the land use change according to land use categories, constructing a machine learning classification model for each group, using the spatial auxiliary variable data corresponding to each group of plots as attribute values, and the land category data after the land use change as label values, and calculating the development probability of each plot through the machine learning classification model. Using each plot as a unit, calculate the neighborhood effect of each plot, as shown in the formula: ,in, For then in the first Time of the first The plot of land is for the first The neighborhood effect of individual plots Let be the distance between the i-th plot and the j-th plot. Indicates the first The area of ​​each plot of land, Indicates the first The area of ​​each plot, d represents the search radius, S max S represents the area of ​​the largest plot of land. min Represents the minimum plot area; construct the vector cellular automata model transformation rule, as shown in the formula: ,in, For the first One plot of land was converted The overall development probability of land parcels will be the first One plot of land was converted The overall development probabilities of land parcels are combined into a vector, thus obtaining the overall development probability array for each land parcel. , For the first One plot of land was converted The development probability of land parcels. No. Class 1 land parcel for the first The neighborhood effect of individual plots For the first One plot of land was converted The development restriction factor for a land parcel is 0 when the parcel is in a restricted development area and 1 when the parcel is in an unrestricted development area. It is a random value.

5. The urban land change simulation method according to claim 4, characterized by, Step S306 includes the following steps: merging the overall development probability array of each land parcel with the transformation probability arrays of each land parcel under the pre-trained DT model, pre-trained RF model, pre-trained SVM model, pre-trained KNN model, pre-trained LSTM model, and pre-trained GRU model, respectively, to obtain the first selection probability array, the second selection probability array, the third selection probability array, the fourth selection probability array, the fifth selection probability array, and the sixth selection probability array for each land parcel, as shown in the formula: , , , , , ,in For the first selection probability array, where For the second selection probability array, where For the third selection probability array, where For the 4th selection probability array, where For the 5th selection probability array, where Choose the probability array for the 6th option. For the first One plot of land was converted The overall development probability of land parcels, This is the transformation probability array for each land parcel under the pre-trained DT model. This is the transition probability array for each land parcel under the pre-trained RF model. This represents the transition probability array for each land parcel under the pre-trained SVM model. This is the transition probability array for each land parcel under the pre-trained KNN model. This is the transition probability array for each land parcel under the pre-trained LSTM model. The conversion probability arrays of each plot under the pre-trained GRU model are used; the roulette wheel algorithm is applied to the 1st, 2nd, 3rd, 4th, 5th and 6th selection probability arrays of each plot to obtain the converted land use categories of each plot under the pre-trained DT model, pre-trained RF model, pre-trained SVM model, pre-trained KNN model, pre-trained LSTM model and pre-trained GRU model.

6. The urban land use change simulation method according to claim 5, characterized by, Step S307 comprises the following steps: according to the converted land use class of each plot under the pre-trained DT model, the pre-trained RF model, the pre-trained SVM model, the pre-trained KNN model, the pre-trained LSTM model and the pre-trained GRU model, the quality factor, the Kappa coefficient and the overall classification accuracy of the plot simulation result under the pre-trained DT model, the pre-trained RF model, the pre-trained SVM model, the pre-trained KNN model, the pre-trained LSTM model and the pre-trained GRU model are calculated respectively, as shown in the following formula: , , wherein, FoM is the quality factor, Kappa is the Kappa coefficient, OA is the overall classification accuracy, A is the area error of simulation under the condition that the actual land use class is changed and the simulated land use class remains unchanged; B is the area error of simulation under the condition that the actual land use class is changed and the simulated land use class is correctly changed; C is the area error of simulation under the condition that the actual land use class is changed and the simulated land use class is incorrectly changed; D is the area error of simulation under the condition that the actual land use class is changed and the simulated land use class is changed; is the consistency between the observed classifier prediction value and the true value, that is, the overall classification accuracy, is the consistency expectation of random classification result; N is the total number of cell units in the study area; is the total number of classes, is a two-dimensional confusion matrix table composed of true labels and predicted labels, is the diagonal element of the confusion matrix, that is, the individual with correct prediction classification; the quality factor, the Kappa coefficient and the overall classification accuracy of the plot simulation result under the pre-trained DT model, the pre-trained RF model, the pre-trained SVM model, the pre-trained KNN model, the pre-trained LSTM model and the pre-trained GRU model are investigated, and the pre-trained model with the quality factor greater than 0.2 and the Kappa coefficient greater than 0.6 is selected, and the pre-trained model with the highest overall classification accuracy is the optimal time series model.

7. The urban land change simulation method according to claim 6, characterized by, Step S3 includes the following steps: inputting the spatial driving factor dataset into the optimal time series model to obtain the transformation probability array for each land parcel. If the optimal time series model is a pre-trained DT model, then it is as shown in the formula: ,in Let L be the transformation probability array for each plot under the pre-trained DT model, and L be the leaf nodes that satisfy the decision tree propagation path of the individual sample whose probability is to be calculated. Let L be the number of samples contained in the leaf node. Let L be the number of samples of class i contained in the leaf node L; if the optimal time series model is a pre-trained RF model, then it is as follows: ,in Let N be the transition probability array for each land parcel under the pre-trained RF model, and let N be the total number of DT models trained by the pre-trained RF model. Let be the probability vector of the nth tree; if the optimal time series model is a pre-trained SVM model, then it is as follows: ,in, This represents the transition probability array for each land parcel under the pre-trained SVM model. This represents the raw output of the pre-trained SVM model. and These are calibration parameters, estimated using the training dataset; if the optimal time-series model is a pre-trained KNN model, the pre-trained KNN model predicts the category of the tested sample by the categories of the K individuals closest to the plot, obtaining the transition probability array of each plot under the pre-trained KNN model, as shown in the formula: ,in, This is the transformation probability array for each land parcel under the pre-trained KNN model. These are the weighting coefficients. , The K individuals closest to the sample. For individuals Euclidean distance to the individual being predicted For an individual with classification label Y, where y represents all possible values ​​for the classification label; if the optimal time series model is a pre-trained LSTM model, then it is as follows: ,in, This is the transformation probability array for each plot of land under the pre-trained LSTM model, which represents the category to which the sample belongs. The probability, It is the first in the score vector One value, It is the sum of the score vectors of all categories after performing an exponential operation, where C is the number of land parcel categories; if the optimal time series model is a pre-trained GRU model, then it is as follows: ,in, is the probability array of each patch under the pre-trained GRU model, that is, the probability of the sample belonging to the category , is the th value in the score vector, is the sum of the exponential operation on the score vector of all categories, and C is the number of patch categories; step S4 includes the following steps: calculating the distance between the centroids of each patch before and after conversion, determining the two patches with the smallest distance as the same patch before and after conversion, and using them for patch matching operation; grouping the patches before land use conversion according to land use categories, constructing a machine learning classification model for each group, taking the spatial auxiliary variable data corresponding to each group of patches as attribute values, and taking the land category data after land use conversion as label values, and calculating the development probability of each patch through the machine learning classification model ; taking each patch as a cell unit, calculating the neighborhood effect of each patch, such as formula: , is the neighborhood effect of the th patch on the th patch at the th moment, is the distance between the i th patch and the j th patch, represents the area of the th patch, represents the area of the th patch, d represents the domain search radius, Smax represents the maximum patch area, and Smin represents the minimum patch area; construct the conversion rule of vector cellular automaton model, such as formula: , is the overall development probability of the th patch converting to the th category patch, and the overall development probability of the th patch converting to the th category patch is integrated into a vector, that is, the overall development probability array of each patch is obtained , is the development probability of the th patch converting to the th category patch, the neighborhood effect of the th category patch on the th patch, is the development restriction factor of the th patch converting to the th category patch, which is 0 when the patch is a restricted development area and 1 when the patch is a non-restricted development area, The random value is obtained; step S5 comprises the following steps: summing and averaging the total development probability array of each land block and the conversion probability array of each land block to obtain a final conversion probability array of each land block; and performing a roulette algorithm operation on the final conversion probability array of each land block respectively to obtain a land use class conversion of each land block.

8. A computer-readable storage medium having stored thereon a computer program, characterized in that, The computer program is executed by a processor to implement the steps of the urban land use change simulation method according to any one of claims 1-7.

9. A computer device comprising a memory, a processor, and a computer program stored on the memory and executable on the processor, characterized in that, The processor executes the program to implement the steps of the urban land use change simulation method according to any one of claims 1-7.

Images