A method for identifying potential boundaries of urban development for national land planning

Abstract

The present invention provides a method for identifying potential boundaries of urban development for national land planning, comprising the following steps: step 1, determining the geographic coordinates of the core points of the city based on multi-year nighttime light remote sensing data; step 2, digitizing the historical map based on remote sensing images and geographic information system technology, and identifying the historical city boundaries; step 3, calculating the distance between the core points of the city and the historical city boundaries; step 4, extracting the extended distance between the radiation lines and the historical city boundaries over many years, and generating a data set; step 5, constructing a radial basis function network; step 6, training the data set based on the radial basis function network, and obtaining the predicted boundary expansion distance by inputting the relevant attributes of the predicted expansion. The present invention can simplify the computational complexity of identifying boundaries, improve the efficiency of identifying urban boundaries, and more accurately detect and simulate urban development boundaries with less data.

Description

Technical Field

[0001] The invention relates to the field of geography and remote sensing science and technology, and in particular to a method for identifying potential boundaries of urban development for land planning. Background Art

[0002] The exploration of urban development boundaries is an important part of national land space planning, and plays an important role in curbing disorderly urban expansion and ensuring the sustainability of economic development. However, there are many factors that affect urban development boundaries, involving different angles and dimensions. Consideration and processing of these factors are the difficulties in identifying urban development boundaries.

[0003] The current methods for determining urban development boundaries can be summarized into two categories: one is qualitative methods, the most common of which is the method based on expert experience judgment, that is, on the basis of comprehensive consideration of the knowledge of experts (or decision makers), the urban development boundaries are determined by artificial vision and manual digital marking; the other is quantitative methods, the most common of which is the combination of multi-agent systems and statistical regression to determine the potential development boundaries of the city.

[0004] However, common qualitative methods may be susceptible to subjective bias, and most quantitative methods have certain requirements for data scale. In summary, the present invention proposes a new method for identifying urban development boundaries for national land planning. Compared with traditional methods, the new method of this patent invention can improve data utilization efficiency, and can more flexibly reflect the impact of various attribute factors on urban boundaries with a smaller sample size, thereby making more effective integrated learning and prediction. Secondly, the operation is relatively simple. By using machine learning algorithms, the need for complex parameter adjustments is reduced, making the model easier to understand and use. The method has certain flexibility and adaptability, and can be flexibly adjusted according to different urban characteristics and development models to adapt to different urban environments and planning needs. Summary of the invention

[0005] The purpose of the present invention is to provide a method for identifying urban development boundaries for land planning. The method utilizes a new method for exploring urban development boundaries, which can simplify the computational complexity of identifying boundaries, improve the efficiency of identifying urban boundaries, and more accurately detect and simulate urban development boundaries with less data.

[0006] This patent proposes a method for identifying potential boundaries of urban development for national land planning, which includes the following steps:

[0007] Step 1: Determine the geographic coordinates of the core points of the city based on years of nighttime light remote sensing data;

[0008] Step 2: digitize historical maps based on remote sensing images and geographic information system technology and identify historical city boundaries;

[0009] Step 3, calculate the distance between the city core point and the historical city boundary;

[0010] Step 4, extract the extended distance between the radiation lines and the historical city boundaries over many years to generate a data set;

[0011] Step 5, constructing a radial basis function network;

[0012] Step 6: Train the data set based on the radial basis function network, and obtain the predicted boundary expansion distance by inputting the relevant attributes of the expansion to be predicted.

[0013] Furthermore, in step 1, the night light data of previous years are superimposed, and the location where the maximum value of the superposition of the night light data of previous years is located is taken as the core point of the city;

[0014] The superposition process is as follows:

[0015]

[0016] In the formula, a represents the starting year, n represents the time span, Count x,y It represents the accumulated value of nighttime light data for n consecutive years starting from the starting year a for the element at the xth row and the yth column. Represents the night light data value of the pixel in the i-th year at the element located in the x-th row and the y-th column, where i = a, a+1, a+2, a+3,…, a+n.

[0017] Furthermore, in step 2, the following steps are also included:

[0018] Step 21, binarization processing is performed on the historical map;

[0019] Step 22, using the regional grouping method, spatially cluster these binary urban cells to form independent urban clusters;

[0020] Step 23, by calculating the coverage area of ​​each urban cluster, select the largest urban cluster as the main urban area, and perform edge extraction.

[0021] Furthermore, in step 21, the binarization process is:

[0022]

[0023] U represents the threshold of urban construction land in land use data, Binary mod(x,y) represents the binarized map, and Pixel(x,y) represents the value of the urban construction land corresponding to the coordinates on the map.

[0024] Furthermore, in step 22, the spatial clustering is:

[0025]

[0026] Cluster(Q k ) represents the city cell Q k The clustering result of its adjacent cells is the urban patch; Q L Indicates that Q k Adjacent city cells, L is Adj(Q k ) is an index for traversal, Adj(Q k ) represents the city cell Q k The set of all adjacent city cells.

[0027] Furthermore, in step 23, the edge extraction is performed as follows:

[0028]

[0029] M represents the selected city cluster, Edge mod (M) and ΔM are the gradient of the cluster, that is, the measure of the degree of change of the function in the urban cluster, which is used to identify the boundary pixels. and Represent the partial derivatives of the urban cluster along the x and y directions respectively.

[0030] Furthermore, in step 3, the point closest to the current origin is the only intersection between the radial line and the boundary, and the distance between the historical city boundary and the city core point is:

[0031]

[0032] In the formula, x c and c is the coordinate of the city core point, x i and i represents the coordinates of the intersection of each azimuth radial line and the city boundary in the i-th year, Distance is the distance between the historical city boundary and the city core point; the average annual extension distance is:

[0033]

[0034] Where, Distance m is the average annual extension distance from the qth year to the pth year, that is, the urban expansion speed, Distance pIndicates the distance between the city boundary and the city core in year p, Distance q It represents the distance between the city boundary and the city core in the qth year, and pq is the time interval.

[0035] Furthermore, in step 4, all data are made to have a uniform spatial resolution and spatial geographic coordinate system using resampling and projection conversion methods, and then normalized:

[0036] s=(s x,y -s min *90%) / (s max *110%-s min *90%

[0037]

[0038] In the formula, s is the attribute grid data, s x,y is the pixel value of each point in the attribute raster data, s min and max are the minimum and maximum pixel values ​​in the raster data, respectively. bq is the average value of a single attribute at the intersection of the historical city boundaries in the previous q years in the direction b, s be is the attribute value in the e-th year in the b direction.

[0039] Furthermore, in step 5, the radial basis function network includes an input layer, a hidden layer and an output layer; each neuron in the hidden layer maps a local area in the input space, which is realized by a Gaussian radial basis function, and the formula is:

[0040]

[0041] h j (x) is the output of the hidden layer neuron, x is the input vector, c j and σ j are the center and width of the jth neuron, respectively, ||xc j || represents the Euclidean distance between the input vector and the neuron center;

[0042] The mathematical representation of the output layer is:

[0043]

[0044] Among them, y′ is the output value of the output layer, w j is the weight and d is the bias value.

[0045] Furthermore, in step 6, for the training set weights, the mean square error is used as the loss function to measure the gap between the predicted value and the true value during the training process:

[0046]

[0047] MSE represents the mean square error as the loss function, y k is the true value of the kth sample, is the model's predicted value for the kth sample, where k represents the index of the sample, ranging from 1 to n;

[0048] In order to minimize the loss, the AdamW optimizer is used, and the formula is as follows:

[0049]

[0050] θ represents the model parameters, and are bias corrections for the first-order and second-order moment estimates, ∈ is a small constant to prevent division by zero, and θ nextstep Indicates the model parameter value after the next update, learning rate indicates the learning rate, and weightdecay indicates the weight decay. is the second momentum decay rate parameter of the AdamW optimizer.

[0051] The beneficial effects achieved by the present invention are:

[0052] The method provided by the present invention uses a new method for exploring urban development boundaries, which can simplify the computational complexity of identifying boundaries, improve the efficiency of identifying urban boundaries, and more accurately detect and simulate urban development boundaries with less data. The algorithm model derived from this can provide data support and decision-making thinking for the decision-making of relevant industry departments. BRIEF DESCRIPTION OF THE DRAWINGS

[0053] Figure 1 .It is a technical roadmap for identifying potential boundaries of urban development for national land planning;

[0054] Figure 2 . Schematic diagram of urban development boundary coordinate information;

[0055] Figure 3 . Comparison of the distance and expansion rate at each azimuth from the city core to the city boundary in 2030;

[0056] Figure 4 .Schematic diagram of city boundaries in previous years and predicted years. DETAILED DESCRIPTION

[0057] The present invention will be further described below in conjunction with specific embodiments, and the advantages and features of the present invention will become clearer as the description proceeds. However, these embodiments are exemplary only and do not constitute any limitation to the scope of the present invention. It should be understood by those skilled in the art that the details and forms of the technical solution of the present invention may be modified or replaced without departing from the spirit and scope of the present invention, but these modifications and replacements all fall within the scope of protection of the present invention.

[0058] As attached Figure 1 As shown, the present invention provides a method for identifying potential boundaries of urban development for land planning, comprising the following steps:

[0059] Step 1: Determine the urban core based on multi-year nighttime light remote sensing data.

[0060] First, the resampling and projection conversion methods are used to make all night light data have a unified spatial resolution and spatial geographic coordinate system; secondly, the night light data of previous years are superimposed as follows:

[0061]

[0062] Where Count x,y Represents the accumulated value of nighttime light data for n consecutive years starting from the starting year a for the element at row x and column y. i Represents the night light data value of the pixel in the element located in the xth row and the yth column in the ith year, where i = a, a+1, a+2, a+3,…, a+n. Take Count x,y The location of the maximum value is taken as the city core.

[0063] The brightness of nighttime lights is often closely related to the density and intensity of human activities. Lights in urban core areas are often more continuous and concentrated in spatial distribution, which easily forms a clear contrast with surrounding areas. By accumulating light data from previous years, these long-standing high brightness points can be captured more accurately, thereby scientifically locating the urban core. The accumulation process enhances the continuity and consistency of the signal, reduces the impact of accidental or temporary changes, and makes the identified urban core points more representative and accurate.

[0064] Step 2: Determination of historical city boundaries.

[0065] Step 21, the historical map is digitized based on remote sensing images and geographic information system (GIS) technology, using improved image binarization technology. In this process, urban areas are marked as 1 (i.e., urban units), while other areas are marked as 0 (non-urban units). This conversion process can be expressed by the following formula, where U represents the threshold of urban construction land in land use data:

[0066]

[0067] Step 22: Use regional grouping to perform spatial clustering on these binary urban cells to form independent urban clusters. This process uses the adjacency principle to merge adjacent urban cells into independent urban patches. This process can be expressed by the following formula, where Q k Represents a city cell, Adj(Q k ) represents and Q k Adjacent cell sets:

[0068]

[0069] Cluster(Q k ) represents the city cell Q k The clustering result of its adjacent cells is the urban patch. L Indicates that Q k Adjacent city cells, L is Adj(Q k ) is an index for traversal. k ) represents the city cell Q k The set of all adjacent city cells.

[0070] Step 23, by calculating the coverage area of ​​each urban cluster, select the largest urban cluster as the main urban area. The city boundary is determined by extracting the edge of this largest cluster using an edge detection algorithm. The edge detection process can be summarized by the following formula:

[0071]

[0072] M represents the selected city cluster, Edge mod (M) and ΔM are the gradient of the cluster, that is, the measure of the degree of change of the function in the urban cluster, which is used to identify the boundary pixels. and Represent the partial derivatives of the city cluster along the x-direction and the y-direction respectively. Through this method, we improve the accuracy of city boundary identification.

[0073] Step 3, the distance between the urban core and the historical urban boundary.

[0074] The radial method is used to describe the urban expansion at a given orientation. That is, taking the urban core as the origin, a separate radial line is generated from each orientation (1° apart, i.e. 360 radial lines).

[0075] Figure 2The parameters involved are the latitude and longitude coordinates of the city core, which are start_x and start_y. The azimuth starts from 0° due north, and radial lines are generated clockwise with each radial line separated by 1°. The distance value can be greater than the city boundary.

[0076] Where the radial line intersects the boundary, multiple intersection points may be generated. The present invention sets a rule that the point closest to the current origin is the only intersection point between the radial line and the boundary.

[0077]

[0078] In the formula, x c and c is the coordinate of the city core point, x i and i Represents the coordinates of the intersection of each azimuth radial line and the city boundary in the i-th year. Distance is the distance between the historical city boundary and the city core. The average annual extension distance is:

[0079]

[0080] Where, Distance m is the average annual extension distance from the qth year to the pth year, that is, the urban expansion speed, Distance p Indicates the distance between the city boundary and the city core in year p, Distance q It represents the distance between the city boundary and the city core in the qth year, and pq is the time interval.

[0081] Step 4: Data extraction and preprocessing.

[0082] Extract the properties of the intersection of the radial line and the historical city boundary (such as altitude, slope, distance between the point and the water body or road, etc.). The average value of the location attribute value of the historical city boundary intersection in each direction in the previous q years is used as the input attribute of the feedforward neural network, and the average annual extension distance Distance from the qth year to the pth year in each direction is used as the input attribute of the feedforward neural network. m As the output target and label, a total of 360 samples of the data set can be obtained.

[0083] Resampling and projection conversion methods are used to make all data have a uniform spatial resolution and spatial geographic coordinate system, and then normalization is performed to make them more comparable and characteristic.

[0084] s=(s x,y -s min *90%) / (s max *110%-s min *90%

[0085]

[0086] In the formula, s is the attribute grid data, s x,y is the pixel value of each point in the attribute raster data, s min and max are the minimum and maximum pixel values ​​in the raster data, respectively. bq is the average value of a single attribute at the intersection of the historical city boundaries in the previous q years in the direction b, s be is the attribute value in the e-th year in the b direction.

[0087] Step 5: Construction of radial basis function (RBF) network.

[0088] The RBF network consists of an input layer, a radial basis hidden layer, and a linear output layer. Each neuron in the hidden layer maps a local area in the input space, which is implemented by the Gaussian radial basis function, and its formula is:

[0089]

[0090] h j (x) is the output of the hidden layer neuron, x is the input vector, c j and σ j are the center and width of the jth neuron, respectively, ||xc j || represents the Euclidean distance between the input vector and the neuron center;

[0091] The mathematical representation of the output layer is:

[0092]

[0093] Among them, y′ is the output value of the output layer, w j is the weight and d is the bias value.

[0094] Step 6: Training and Validation

[0095] The present invention divides the data set into 80% training set and 20% validation set. For the training set weights, the mean square error (MSE) is used as the loss function to measure the gap between the predicted value and the true value during the training process:

[0096]

[0097] In order to minimize the loss, we chose the AdamW optimizer, which adds weight decay to the Adam optimizer. The formula is as follows:

[0098]

[0099] θ represents the model parameters, and are bias corrections for the first-order and second-order moment estimates, respectively, and ∈ is a small constant to prevent division by zero. next step Indicates the model parameter value after the next update, that is, the updated weight parameter. Learning rate indicates the learning rate, which is used to control the step size of parameter update. Weight decay indicates weight decay, which is used to reduce the overfitting of the model by adding an additional term to the loss function to limit the larger weight value. The second momentum decay rate parameter of the AdamW optimizer is used to calculate the exponentially weighted moving average of the second-order moment estimate. This parameter is fixed and represents the update rate used to adjust the second-order moment estimate.

[0100] The above network is iteratively trained, and the network parameters are adjusted according to different results. After the training is completed, the intersection attributes such as the altitude, slope, and distance between the point and the water body or road to be predicted are input, and the predicted boundary expansion distance ( Figure 3 ).

[0101] Based on the above method steps, the urban expansion boundary of a typical city in 2030 ( Figure 4 ), using the attribute data from 2000 to 2010 as input data and the average extended distance from 2000 to 2020 (Distance m ) as the output for model training and verification, and then the attribute data from 2000 to 2020 are imported into the model to obtain the average annual extension distance from 2000 to 2030 (Distance m ), which is the urban expansion speed at each azimuth from 2000 to 2030. The urban expansion speed from 2000 to 2030 multiplied by the interval year (30) equals the urban expansion distance at each azimuth from 2000 to 2030.

[0102] The present invention is not limited to the above-mentioned specific implementation modes. A person skilled in the art may implement the present invention in various other specific implementation modes according to the embodiments and the disclosure of the drawings. Therefore, any design that adopts the design structure and concept of the present invention and makes some simple transformations or changes shall fall within the scope of protection of the present invention.

Claims

1. A method for identifying potential boundaries of urban development for land planning, characterized in that: The method for identifying potential boundaries of urban development for national land planning comprises the following steps: Step 1: Determine the geographic coordinates of the core points of the city based on years of nighttime light remote sensing data; Step 2: digitize historical maps based on remote sensing images and geographic information system technology and identify historical city boundaries; Step 3, calculate the distance between the city core point and the historical city boundary; Step 4, extract the extended distance between the radiation lines and the historical city boundaries over many years to generate a data set; Step 5, constructing a radial basis function network; Step 6: Train the data set based on the radial basis function network, and input the relevant attributes of the predicted expansion to obtain the predicted boundary expansion distance; In step 3, the point closest to the current origin is the only intersection between the radial line and the boundary. The distance between the historical city boundary and the city core point is: In the formula, x c and c is the coordinate of the city core point, x i and i represents the coordinates of the intersection of each azimuth radial line and the city boundary in the i-th year, Distance is the distance between the historical city boundary and the city core point; the average annual extension distance is: Where, Distance m is the average annual extension distance from the qth year to the pth year, that is, the urban expansion speed, Distance p Indicates the distance between the city boundary and the city core in year p, Distance q represents the distance between the city boundary and the city core in the qth year, and pq is the time interval; In step 4, the intersection attributes of the radial line and the historical city boundary are extracted; the intersection attributes include altitude, slope, and the distance between the point and the water body or road; the average value of the location attribute value of the historical city boundary intersection in each direction in the previous q years is used as the input attribute of the feedforward neural network, and the average annual extension distance Distance from the qth year to the pth year in each direction is used as the input attribute of the feedforward neural network. m As output targets and labels, build a dataset; Resampling and projection conversion methods are used to make all data have a uniform spatial resolution and spatial geographic coordinate system, and then normalization is performed: s=(s x,y -s min *90%) / (s max *110%-s min *90%) In the formula, s is the attribute grid data, s x,y is the pixel value of each point in the attribute raster data, s min and max are the minimum and maximum pixel values ​​in the raster data, respectively. bq is the average value of a single attribute at the intersection of the historical city boundaries in the previous q years in the direction b, s be is the attribute value in the e-th year in the b direction.

2. The method for identifying potential boundaries of urban development for national land planning according to claim 1 is characterized in that: In step 1, the night light data of previous years are superimposed, and the location where the maximum value of the superposition of the night light data of previous years is located is taken as the core point of the city; The superposition process is as follows: In the formula, a represents the starting year, n represents the time span, Count x,y It represents the accumulated value of nighttime light data for n consecutive years starting from the starting year a for the element at the xth row and the yth column. Represents the night light data value of the pixel in the i-th year at the element located in the x-th row and the y-th column, where i = a, a+1, a+2, a+3,…, a+n.

3. The method for identifying potential boundaries of urban development for national land planning according to claim 1 is characterized in that: In step 2, the following steps are also included: Step 21, binarization processing is performed on the historical map; Step 22, using the regional grouping method, spatially cluster these binary urban cells to form independent urban clusters; Step 23, by calculating the coverage area of ​​each urban cluster, select the largest urban cluster as the main urban area, and perform edge extraction.

4. The method for identifying potential boundaries of urban development for national land planning according to claim 3 is characterized in that: In step 21, the binarization process is: U represents the threshold of urban construction land in land use data, Binary mod (x,y) represents the binarized map, and Pixel(x,y) represents the value of the urban construction land corresponding to the coordinates on the map.

5. The method for identifying potential boundaries of urban development for national land planning according to claim 3 is characterized in that: In step 22, the spatial clustering is: Cluster(Q k ) represents the city cell Q k The clustering result of its adjacent cells is the urban patch; Q L Indicates that Q k Adjacent city cells, L is Adj(Q k ) is an index for traversal, Adj(Q k ) represents the city cell Q k The set of all adjacent city cells.

6. The method for identifying potential boundaries of urban development for national land planning according to claim 3 is characterized in that: In step 23, the edge extraction is performed as follows: M represents the selected city cluster, Edge mod (M) and ΔM are the gradient of the cluster, that is, the measure of the degree of change of the function in the urban cluster, which is used to identify the boundary pixels. and Represent the partial derivatives of the urban cluster along the x and y directions respectively.

7. The method for identifying potential boundaries of urban development for national land planning according to claim 1 is characterized in that: In step 5, the radial basis function network includes an input layer, a hidden layer and an output layer; each neuron in the hidden layer maps a local area in the input space, which is realized by a Gaussian radial basis function, and its formula is: h j (x) is the output of the hidden layer neuron, x is the input vector, c j and σ j are the center and width of the jth neuron, respectively, ||xc j || represents the Euclidean distance between the input vector and the neuron center; The mathematical representation of the output layer is: Among them, y′ is the output value of the output layer, w j is the weight and d is the bias value.

8. The method for identifying potential boundaries of urban development for national land planning according to claim 1 is characterized in that: In step 6, for the training set weights, the mean square error is used as the loss function to measure the gap between the predicted value and the true value during the training process: MSE represents the mean square error as the loss function, y k is the true value of the kth sample, is the model's predicted value for the kth sample, where k represents the index of the sample, ranging from 1 to n; In order to minimize the loss, the AdamW optimizer is used, and the formula is as follows: θ represents the model parameters, and are bias corrections for the first-order and second-order moment estimates, ∈ is a small constant to prevent division by zero, and θ nextstep Indicates the model parameter value after the next update, learning rate indicates the learning rate, and weight decay indicates the weight decay. is the second momentum decay rate parameter of the AdamW optimizer.

Images