Method for grid identification and classification of urban functional areas based on multi-source spatio-temporal data

Abstract

The present invention provides a method for grid-based identification and classification of urban functional areas based on multi-source spatio-temporal data, which relates to the field of regional identification technology. The method includes dynamically dividing the urban area using an adaptive grid division algorithm, obtaining multi-source spatio-temporal data of the grid-based area and performing spatio-temporal alignment to generate a multi-dimensional time series feature set; extracting multi-dimensional feature vectors of grid cells through a multi-modal feature fusion algorithm, obtaining long-range time series features using a bidirectional gated recurrent unit, and calculating spatial correlation features based on a multi-head attention mechanism. Then, the two types of features are input into a feature fusion network with skip connections to generate enhanced spatio-temporal features; an optimization strategy is adopted. First, the density peak clustering is used to determine the clustering center points of functional areas to obtain an initial division result, the initial result is optimized by spectral clustering with spatial adjacency constraints, and the simulated annealing algorithm is used to optimize the boundaries of functional areas to achieve accurate identification and classification of urban functional areas.

Description

Technical Field

[0001] The present invention relates to the technical field of area recognition, and in particular to a method for grid-based recognition and classification of urban functional areas based on multi-source spatio-temporal data. Background Art

[0002] The recognition of urban functional areas is an important basic task for urban planning and management. Traditional functional area recognition methods mainly rely on manual surveys and empirical judgments. This method not only consumes time and effort, but also is difficult to reflect the dynamic changes of urban functions in a timely manner. With the extensive collection and accumulation of multi-source urban data, the automatic recognition method of urban functional areas based on spatio-temporal big data has gradually become a research hotspot. Existing data-driven methods usually adopt a grid division method with a fixed size, combine multi-source data such as POI data, human movement trajectory data, and remote sensing image data, and realize the automatic recognition and classification of functional areas through feature extraction and machine learning algorithms.

[0003] However, there are still several deficiencies in the existing technology. The grid division method with a fixed size is difficult to adapt to the spatial heterogeneity of urban area density, and it is easy to cause recognition deviation of functional area boundaries; the existing methods simply splice the features of multi-source data and fail to fully explore the spatio-temporal correlation between data, resulting in insufficient feature expression ability; the clustering process of functional areas lacks consideration of the spatial adjacency relationship of grid cells, and it is easy to produce fragmented spatial recognition results; during the optimization process of functional area boundaries, the internal consistency of regions and the differences between regions are not taken into account simultaneously, affecting the accuracy of the final recognition results.

[0004] In summary, there is a need for a method for grid-based recognition and classification of urban functional areas based on multi-source spatio-temporal data, which solves the problem of spatial heterogeneity through an adaptive grid division algorithm, extracts temporal features by using multi-modal feature fusion and bidirectional gated recurrent units, captures spatial correlation features by using a multi-head attention mechanism, and realizes accurate functional area recognition by combining density peak clustering, spectral clustering, and simulated annealing boundary optimization algorithms. The present invention can solve the problems in the existing technology. Summary of the Invention

[0005] An embodiment of the present invention provides a method for grid-based recognition and classification of urban functional areas based on multi-source spatio-temporal data, which can solve the problems in the existing technology.

[0006] In the first aspect of the embodiment of the present invention,

[0007] There is provided a method for grid-based recognition and classification of urban functional areas based on multi-source spatio-temporal data, including:

[0008] The urban area is divided into multiple grid cells by using an adaptive grid division algorithm, and the size of the grid cells is dynamically adjusted based on the regional density threshold to generate a grid-based urban area; multi-source spatio-temporal data of the grid-based urban area is obtained, and spatio-temporal alignment is performed according to a preset number of time segments to generate a multi-dimensional time series feature set; a multi-modal feature fusion algorithm is used to perform feature extraction and fusion operations on the multi-dimensional time series feature set to generate a multi-dimensional feature vector of the grid cell;

[0009] Based on the multi-dimensional feature vector of the grid cell, a bidirectional gated recurrent unit is used to extract long-range time series features; based on the multi-head attention mechanism, the multi-scale spatial correlation coefficients between the grid cell and its surrounding grid cells are calculated and feature aggregation is performed to generate spatial correlation features; the long-range time series features and the spatial correlation features are input into a feature fusion network with skip connections to generate enhanced spatio-temporal features;

[0010] The local density value and the relative distance value of the enhanced spatio-temporal features are calculated, and density peak clustering is performed to determine the clustering center points of the functional areas, generating an initial functional area division result; an affinity matrix is constructed using the enhanced spatio-temporal features, and spectral clustering operations are performed in combination with spatial adjacency constraints to generate an optimized functional area division result; through a simulated annealing boundary optimization algorithm, the indicators of regional internal consistency and inter-regional difference are iteratively optimized to generate the final urban functional area division result.

[0011] In an alternative embodiment,

[0012] Dividing the urban area into multiple grid cells by using an adaptive grid division algorithm and dynamically adjusting the size of the grid cells based on the regional density threshold to generate a grid-based urban area includes:

[0013] Obtain the urban data of the target urban area, where the urban data includes population flow density data, building density data, land use type data, mobile phone signal data, and environmental perception data, perform data cleaning and normalization preprocessing on the urban data to generate an urban activity feature matrix;

[0014] The target urban area is divided into initial grid cells of a preset maximum size, a multi-dimensional density evaluation model is constructed based on the urban activity feature matrix, a deep autoencoder is used to extract the hidden layer features of each initial grid cell, the regional density distribution of each initial grid cell is calculated in combination with variational inference, the regional density distribution is input into a graph attention network, and the spatial dependence relationship between adjacent initial grid cells is calculated to generate the grid optimization function value of each initial grid cell;

[0015] Construct an adaptive density threshold learning model based on the grid optimization function value, use a recurrent neural network to model the temporal variation of the grid optimization function value, dynamically adjust the regional density threshold sequence in combination with reinforcement learning, divide the initial grid cells into multiple density levels according to the regional density threshold sequence, and for each density level, optimize and determine the grid size using the gradient descent method, where the grid size is dynamically adjusted within the range of a preset maximum grid size and a minimum grid size;

[0016] Perform adaptive subdivision on the initial grid cells through a quadtree structure, and perform grid splitting and merging based on the regional density threshold sequence to generate a grid-based urban area.

[0017] In an alternative embodiment,

[0018] Performing grid splitting and merging based on the regional density threshold sequence to generate a grid-based urban area includes:

[0019] When the grid optimization function value is greater than the corresponding regional density threshold and the spatial dependence is less than a first preset threshold, perform non-uniform four-way splitting on the initial grid cell. When the grid optimization function value of adjacent grid cells is less than the corresponding regional density threshold and the spatial dependence is greater than a second preset threshold, perform adaptive merging on the adjacent grid cells;

[0020] The non-uniform four-way splitting includes: calculating the density gradient field within the initial grid cell, determining the main direction and secondary direction of density change based on the density gradient field; calculating the relative importance measure at each position within the initial grid cell and constructing a weight matrix;

[0021] In the main direction, based on the cumulative distribution function of the weight matrix, determine the position of the first splitting line such that the sum of the weighted density values on both sides of the first splitting line is equal; in the secondary direction, calculate the position of the second splitting line for each of the two sub-regions divided by the first splitting line such that the weighted density values of the four sub-grids satisfy the balance constraint;

[0022] Divide the initial grid cell into four sub-grids, where the area of each sub-grid is inversely proportional to the average density, and the area ratio of adjacent sub-grids does not exceed a preset maximum ratio threshold;

[0023] The adaptive merging includes: calculating the shape similarity and density difference degree of adjacent grid cells, and merging adjacent grid cells with a shape similarity greater than a third preset threshold and a density difference degree less than a fourth preset threshold into a new grid cell, and obtaining the boundary of the new grid cell by fitting the boundary points of the adjacent grid cells using the least squares method.

[0024] In an alternative embodiment,

[0025] Obtain multi-source spatio-temporal data of a grid-based urban area, perform spatio-temporal alignment according to a preset number of time segments, and generate a multi-dimensional time series feature set; perform feature extraction and fusion operations on the multi-dimensional time series feature set through a multi-modal feature fusion algorithm to generate a multi-dimensional feature vector of grid cells, including:

[0026] The multi-source spatio-temporal data includes POI data, remote sensing image data, and population trajectory data of the grid-based urban area;

[0027] Within a time segment, generate a POI status sequence based on the POI data, generate a ground object distribution sequence based on the remote sensing image data, and generate a population distribution sequence based on the population trajectory data; perform noise filtering on the population distribution sequence through Kalman filtering, and supplement missing points through piecewise cubic spline interpolation to generate population density time series data; perform atmospheric correction and geometric correction on the ground object distribution sequence to generate ground object coverage time series data; align the POI status sequence, the ground object coverage time series data, and the population density time series data to a unified spatio-temporal reference system to generate spatio-temporally aligned data;

[0028] Construct a multi-level POI semantic tree from the POI status sequence, calculate the semantic information entropy of each layer in the multi-level POI semantic tree, and extract the semantic feature vector of grid cells layer by layer from the top layer to the bottom layer through recursive morphological decomposition;

[0029] Construct a multi-scale pyramid from the ground object coverage time series data, and extract the scene feature vector of grid cells through bottom-up feature aggregation and top-down boundary optimization of the multi-scale pyramid;

[0030] Calculate the time-varying kernel density based on the population density time series data, and extract the flow feature vector of grid cells through adaptive kernel density estimation;

[0031] Combine the semantic feature vector of grid cells, the scene feature vector of grid cells, and the flow feature vector of grid cells to generate a multi-dimensional time series feature set, and perform standardization processing;

[0032] Construct a reconstruction error minimization model to calculate the feature fusion weights, and perform weighted combination of the standardized features according to the feature fusion weights to generate a multi-dimensional feature vector of grid cells.

[0033] In an alternative embodiment,

[0034] Based on the multi-dimensional feature vector of grid cells, use a bidirectional gated recurrent unit to extract long-range time series features, including:

[0035] Stratify and sample the multi-dimensional feature vectors of grid cells according to daily, weekly, and monthly spans to obtain multi-scale time-series feature sequences; for each feature sequence of each time scale in the multi-scale time-series feature sequences, perform feature extraction through a bidirectional gated recurrent unit, which is coupled with a spiking neuron module, and the bidirectional gated recurrent unit includes a forward gated recurrent unit and a backward gated recurrent unit;

[0036] The update gate output of the forward gated recurrent unit is converted into a first spike train through a first spike encoder, which is based on a threshold integrate-and-fire model, generates spikes when the accumulated input potential exceeds a preset firing threshold, and resets the membrane potential; the reset gate output of the forward gated recurrent unit is converted into a second spike train through a second spike encoder, which is the same as the first spike encoder; generate a forward hidden state sequence according to the firing times of the first spike train and the second spike train;

[0037] The update gate output and the reset gate output of the backward gated recurrent unit are respectively converted through the first spike encoder and the second spike encoder, and a backward hidden state sequence is generated according to the firing times of the converted spike trains;

[0038] Calculate the time intervals between adjacent spikes in the forward hidden state sequence to obtain a first spike firing intensity, calculate the time intervals between adjacent spikes in the backward hidden state sequence to obtain a second spike firing intensity; convert the first spike firing intensity and the second spike firing intensity into a first attention score and a second attention score respectively, and perform normalization processing to obtain a forward attention weight and a backward attention weight respectively;

[0039] Perform weighted calculation on the forward hidden state sequence and the backward hidden state sequence with the corresponding forward attention weight and backward attention weight respectively to obtain a fused hidden state sequence; perform temporal convolutional processing on the fused hidden state sequence to obtain local temporal features; perform a non-linear transformation on the multi-dimensional feature vectors of the grid cells to obtain residual features; add the local temporal features and the residual features to obtain enhanced temporal features; perform normalization processing on the enhanced temporal features to obtain long-range temporal features.

[0040] In an alternative embodiment,

[0041] Calculate the local density value and the relative distance value of the enhanced spatio-temporal features, and perform density peak clustering to determine the functional area clustering center points, and generate an initial functional area division result including:

[0042] Calculate the Euclidean distance, Mahalanobis distance, and cosine distance for the enhanced spatio-temporal features, determine the weighting coefficients through the variance contribution rate of the enhanced spatio-temporal features, and perform weighted fusion to obtain a comprehensive distance matrix;

[0043] Sort the comprehensive distance matrix in ascending order, calculate the kernel density estimation curve of the comprehensive distance matrix, obtain all local minimum points of the kernel density estimation curve, and select the first distance value that satisfies the derivative of the local minimum point being negative and the amplitude being less than the mean of the kernel density estimation curve as the truncation distance;

[0044] Calculate the distribution density of the comprehensive distance matrix, determine the bandwidth parameter of the adaptive Gaussian kernel function, construct the adaptive Gaussian kernel function using the bandwidth parameter, and calculate the local density value of each grid cell based on the truncation distance and the adaptive Gaussian kernel function; Calculate the initial relative distance value of each grid cell based on the local density value, calculate the minimum distance between the grid cell and the grid cell with a higher local density value under multiple preset search radii, and perform weighted fusion based on the stability of the initial relative distance value under different search radii to determine the relative distance value;

[0045] Take the product of the local density value and the relative distance value as the clustering center score, sort the clustering center scores in descending order, select the grid cells in the top thirty percent that satisfy the spatial distribution uniformity as the functional area clustering center points, and at the same time determine the non-functional area clustering center points;

[0046] Sort all grid cells in descending order according to the local density value. For each grid cell of the non-functional area clustering center point, assign the grid cell to the nearest functional area clustering center point according to the density reachability of the local density value to generate an initial functional area division result.

[0047] In an alternative embodiment,

[0048] Construct an affinity matrix using the enhanced spatio-temporal features, perform spectral clustering operations in combination with spatial adjacency constraints to generate an optimized functional area division result; Through the simulated annealing boundary optimization algorithm, iteratively optimize the intra-region consistency and inter-region difference indicators to generate the final urban functional area division result, including:

[0049] Calculate the similarity value between the enhanced spatio-temporal features through the radial basis function to construct a similarity matrix; Determine the adjacency strength between grid cells according to the spatial position relationship of the grid cells to construct a spatial adjacency constraint matrix;

[0050] Determine the weight coefficient according to the variance contribution rate of the enhanced spatio-temporal features and the spatial adjacency constraint matrix, and perform weighted combination of the similarity matrix and the spatial adjacency constraint matrix to obtain the affinity matrix; calculate the degree matrix of the affinity matrix, and perform Laplacian normalization on the affinity matrix using the degree matrix to obtain the normalized affinity matrix; perform eigen-decomposition on the normalized affinity matrix to obtain an eigenvalue sequence, determine the number of eigenvectors based on the descending order of the eigenvalue sequence, select the corresponding eigenvectors to construct an eigen-matrix; perform normalization processing on the row vectors of the eigen-matrix to obtain low-dimensional representation features, and perform clustering operations in the feature space composed of the low-dimensional representation features to reassign grid cells to different functional areas, generating an optimized functional area division result;

[0051] Calculate the variance of the enhanced spatio-temporal features of the grid cells within each functional area in the optimized functional area division result to obtain the regional internal consistency index value; calculate the difference degree of the enhanced spatio-temporal features between adjacent functional areas in the optimized functional area division result to obtain the inter-regional difference index value; combine the regional internal consistency index value and the inter-regional difference index value to construct a simulated annealing optimization objective function;

[0052] Based on a preset initial temperature and cooling coefficient, start the simulated annealing operation, select boundary grid cells from the optimized functional area division result, randomly change the functional area category of the boundary grid cells to determine the changed scheme, calculate the function values of the simulated annealing optimization objective function before and after the change, determine the difference. When the difference is less than zero, accept the changed scheme. When the difference is greater than zero, calculate the acceptance probability based on the cooling coefficient and the current temperature, and at the same time determine the new temperature, and determine whether to accept the changed scheme based on the acceptance probability; repeat the execution until the new temperature reaches the preset stop temperature, generating the final urban functional area division result.

[0053] In the second aspect of the embodiments of the present invention,

[0054] Provide a grid-based identification and classification system for urban functional areas based on multi-source spatio-temporal data, including:

[0055] The first unit is used to divide the urban area into multiple grid cells by using an adaptive grid division algorithm, dynamically adjust the size of the grid cells based on the regional density threshold, generating a grid-based urban area; obtain the multi-source spatio-temporal data of the grid-based urban area, perform spatio-temporal alignment according to the preset number of time segments, generating a multi-dimensional time series feature set; perform feature extraction and fusion operations on the multi-dimensional time series feature set through a multi-modal feature fusion algorithm, generating grid cell multi-dimensional feature vectors;

[0056] The second unit is used to extract long-term temporal features by using a bidirectional gated recurrent unit based on the multi-dimensional feature vectors of grid cells; calculate the multi-scale spatial correlation coefficients between grid cells and surrounding grids based on the multi-head attention mechanism and perform feature aggregation to generate spatial correlation features; input the long-term temporal features and spatial correlation features into a feature fusion network with skip connections to generate enhanced spatio-temporal features;

[0057] The third unit is used to calculate the local density value and relative distance value of the enhanced spatio-temporal features, perform density peak clustering to determine the clustering center points of functional areas, and generate an initial functional area division result; construct an affinity matrix by using the enhanced spatio-temporal features, and perform spectral clustering operation in combination with spatial adjacency constraints to generate an optimized functional area division result; through the simulated annealing boundary optimization algorithm, iteratively optimize the indexes of regional internal consistency and inter-regional difference to generate the final urban functional area division result.

[0058] In the third aspect of the embodiments of the present invention,

[0059] A kind of electronic device is provided, including:

[0060] A processor;

[0061] A memory for storing instructions executable by the processor;

[0062] Wherein, the processor is configured to call the instructions stored in the memory to execute the method described above.

[0063] In the fourth aspect of the embodiments of the present invention,

[0064] A computer-readable storage medium is provided, on which computer program instructions are stored, and when the computer program instructions are executed by a processor, the method described above is implemented.

[0065] In the embodiments of the present invention, through adaptive grid division and spatio-temporal alignment processing of multi-source spatio-temporal data, refined grid division and feature expression of urban areas are realized, the accuracy and spatio-temporal consistency of data processing are improved, and a reliable data basis is provided for functional area recognition; by combining a bidirectional gated recurrent unit and a multi-head attention mechanism, the long-term temporal dependence relationship and multi-scale spatial correlation features of grid cells are effectively captured, and the expression ability of spatio-temporal features is enhanced through a feature fusion network, improving the accuracy and robustness of functional area recognition; by combining density peak clustering and spectral clustering methods and introducing a simulated annealing boundary optimization algorithm, adaptive division and boundary optimization of functional areas are realized, effectively balancing regional internal consistency and inter-regional difference, and improving the scientificity and rationality of the functional area division result. Description of the Drawings

[0066] Figure 1Schematic flowchart of the method for grid-based identification and classification of urban functional areas based on multi-source spatio-temporal data according to an embodiment of the present invention;

[0067] Figure 2 Line chart of grid deformation rate under various density gradient conditions;

[0068] Figure 3 Spatial distribution effect diagram of urban functional area division;

[0069] Figure 4 Schematic structural diagram of the grid-based identification and classification system of urban functional areas based on multi-source spatio-temporal data according to an embodiment of the present invention. Detailed implementation manners

[0070] To make the objectives, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Apparently, the described embodiments are only a part rather than all of the embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0071] The technical solutions of the present invention will be described in detail below with specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be repeated in some embodiments.

[0072] Figure 1 Schematic flowchart of the method for grid-based identification and classification of urban functional areas based on multi-source spatio-temporal data according to an embodiment of the present invention, as Figure 1 shown, the method includes:

[0073] S101. Use an adaptive grid division algorithm to divide the urban area into multiple grid units, dynamically adjust the size of the grid units based on the regional density threshold, and generate a grid-based urban area; obtain the multi-source spatio-temporal data of the grid-based urban area, perform spatio-temporal alignment according to the preset number of time segments, and generate a multi-dimensional time series feature set; perform feature extraction and fusion operations on the multi-dimensional time series feature set through a multi-modal feature fusion algorithm to generate a multi-dimensional feature vector of the grid unit;

[0074] In this embodiment, through the adaptive grid division algorithm, the grid cell size is dynamically adjusted according to the regional density, making the urban area division more in line with the actual spatial distribution, and improving the fineness and rationality of spatial representation; the grid cells are divided according to preset time segments, and multi-source spatio-temporal data is spatio-temporally aligned, effectively capturing the dynamic characteristics changing with time within the urban area and ensuring the spatio-temporal consistency and integrity of the data; the multi-modal feature fusion algorithm is used to extract and fuse multi-dimensional temporal features, and the generated multi-dimensional feature vectors can comprehensively reflect the spatio-temporal information within each grid cell, providing an accurate and reliable feature basis for subsequent urban functional area identification and classification.

[0075] S102. Based on the multi-dimensional feature vectors of grid cells, a bidirectional gated recurrent unit is used to extract long-range temporal features; based on the multi-head attention mechanism, the multi-scale spatial correlation coefficients between grid cells and their surrounding grid cells are calculated and feature aggregation is performed to generate spatial correlation features; the long-range temporal features and spatial correlation features are input into a feature fusion network with skip connections to generate enhanced spatio-temporal features;

[0076] In a specific implementation manner, after obtaining the multi-dimensional feature vectors of grid cells, temporal feature extraction is first performed using a bidirectional gated recurrent unit. The bidirectional gated recurrent unit includes two processing units, a forward unit and a backward unit, which extract features from the forward and backward directions of the time series respectively. Two control units, an update gate and a reset gate, are set inside each processing unit. The update gate is used to control the retention degree of historical information, and the reset gate is used to control the reset degree of historical information. For the input multi-dimensional feature vector sequence of grid cells, through the cooperation of the update gate and the reset gate, temporal features containing long-term dependencies are extracted. Finally, the output features of the forward and backward processing units are fused to obtain the long-range temporal features of the grid cells. Taking a certain grid cell as an example, after being processed by the bidirectional gated recurrent unit, its 24-hour feature sequence can effectively capture the activity patterns and their evolution characteristics at different time periods.

[0077] Then, the spatial correlation between grid cells is processed based on the multi-head attention mechanism. First, adjacent grid cells within a radius of 500 meters are determined for each grid cell as the correlation calculation range. Eight attention heads are set, and each attention head independently calculates the spatial correlation coefficient. In a single attention head, the similarity between the feature vector of the target grid cell and the feature vectors of its surrounding grid cells is calculated to obtain the attention weight scores. After the attention weight scores are normalized, they are weighted and combined with the feature vectors of the surrounding grid cells. The output features of the eight attention heads are concatenated and linearly transformed to obtain a feature representation reflecting the spatial correlation at multiple scales. Based on this feature representation, the features of the target grid cell and its surrounding grid cells are weighted and aggregated to generate the spatial correlation features of the grid cells.

[0078] Finally, a feature fusion network with skip connections is constructed to achieve the deep fusion of long-range temporal features and spatial correlation features. The feature fusion network consists of 4 convolutional layers, and skip connections are added between adjacent layers. The network input includes the long-range temporal features and spatial correlation features obtained above. Through the skip connection structure, the fine-grained features of the lower layer can be directly transmitted to the higher layer to avoid information loss. At the same time, the higher-layer features can also guide the extraction process of the lower-layer features. After multi-layer feature extraction and fusion, enhanced spatio-temporal features with a dimension of 128 are finally output, which comprehensively express the temporal evolution law and spatial correlation characteristics of grid cells. In a specific implementation, taking a certain area as an example, the fused features can effectively distinguish different types of urban functional areas such as commercial areas, residential areas, and office areas.

[0079] In this embodiment, a bidirectional gated recurrent unit is used to process the multi-dimensional features of grid cells, which can fully capture long-range temporal information and improve the modeling ability of temporal dynamic changes; based on the multi-head attention mechanism, multi-scale spatial correlation coefficients are calculated, and spatial correlation features are extracted through feature aggregation operations, effectively reflecting the complex spatial relationship between grid cells and surrounding areas; through the feature fusion network with skip connections, the long-range temporal features and spatial correlation features are fused, not only retaining the original information, but also enhancing the expression ability of spatio-temporal features, thus providing more accurate and comprehensive feature support for subsequent urban functional area recognition.

[0080] S103. Calculate the local density value and relative distance value of the enhanced spatio-temporal features, perform density peak clustering to determine the clustering center points of functional areas, and generate the initial functional area division result; construct an affinity matrix using the enhanced spatio-temporal features, and perform spectral clustering operations in combination with spatial adjacency constraints to generate an optimized functional area division result; through the simulated annealing boundary optimization algorithm, iteratively optimize the intra-region consistency and inter-region difference indicators to generate the final urban functional area division result.

[0081] In this embodiment, the local density and relative distance are calculated using the enhanced spatio-temporal features, and density peak clustering is used to determine the clustering center points of functional areas, thereby generating a preliminary functional area division result, which realizes the fine capture of local features within the region; an affinity matrix is constructed and combined with the spatial adjacency constraints of grid cells, and spectral clustering operations are used to optimize the initial division, effectively integrating spatial correlation information and improving the rationality and accuracy of the division result; through the simulated annealing boundary optimization algorithm, the intra-region consistency and inter-region difference indicators are iteratively optimized, and finally the optimized urban functional area division result is obtained, achieving the balance of intra-region homogeneity and significant inter-region differences, and providing scientific and accurate data support for urban planning and management.

[0082] In an optional implementation, an adaptive grid division algorithm is used to divide the urban area into a plurality of grid cells, and the size of the grid cells is dynamically adjusted based on a regional density threshold, and the generation of a gridded urban area includes:

[0083] Acquire urban data of the target urban area, the urban data including pedestrian density data, building density data, land use type data, mobile phone signal data and environmental perception data, perform data cleaning and normalization preprocessing on the urban data, and generate an urban activity feature matrix;

[0084] The target urban area is divided into initial grid units of a preset maximum size, a multidimensional density assessment model is constructed based on the urban activity feature matrix, a deep autoencoder is used to extract the hidden layer features of each initial grid unit, and the regional density distribution of each initial grid unit is calculated in combination with variational inference. The regional density distribution is input into a graph attention network, the spatial dependency relationship between adjacent initial grid units is calculated, and a grid optimization function value of each initial grid unit is generated;

[0085] An adaptive density threshold learning model is constructed based on the grid optimization function value, a recursive neural network is used to model the time series change of the grid optimization function value, and a regional density threshold sequence is dynamically adjusted in combination with reinforcement learning. The initial grid unit is divided into a plurality of density levels according to the regional density threshold sequence. For each density level, a gradient descent method is used to optimize and determine the grid size, and the grid size is dynamically adjusted within a preset maximum grid size and a minimum grid size.

[0086] The initial grid unit is adaptively subdivided through a quadtree structure, and grid splitting and merging are performed based on the regional density threshold sequence to generate a gridded urban area.

[0087] In a specific implementation, firstly, the urban data of the target urban area is obtained, including the crowd density data, building density data, land use type data, mobile phone signal data and environmental perception data. The crowd density data is collected through the video surveillance system in each area of ​​the city, and the crowd flow data is recorded every 5 minutes; the building density data includes information such as building area and volume ratio; the land use type data includes type identifications such as commercial land, residential land, and industrial land; the mobile phone signal data records the distribution of base station signal strength; and the environmental perception data includes environmental parameters such as temperature, humidity, and PM2.5. The collected raw data is cleaned, outliers and missing values ​​are removed, and the maximum and minimum value normalization method is used to unify the data of different dimensions into the range of 0-1 to generate a standardized urban activity feature matrix.

[0088] Next, the target urban area is divided into initial grid cells, and the initial grid size is set to 500 meters × 500 meters. A multi-dimensional density evaluation model is constructed. A 5-layer deep autoencoder is used to extract the hidden layer features of each grid cell. The number of neurons in each layer of the encoder is 256, 128, 64, 32, and 16 in sequence, and the decoder is symmetrically set. Based on the extracted hidden layer features, the variational inference is combined to calculate the regional density distribution of the grid cells. Specifically, the Gaussian mixture model is used to estimate the density probability distribution. The density distribution is input into the graph attention network. The attention layer is set with 8 attention heads, and the output dimension of each attention head is 64. The softmax function is used to calculate the spatial dependence weight between adjacent grid cells, and the grid optimization function value is generated.

[0089] Then, an adaptive density threshold learning model is constructed based on the grid optimization function value. The long short-term memory network is used to model the grid optimization function value. The network contains 2 layers of LSTM layers, and the hidden layer dimension is 128. Combined with the deep deterministic policy gradient algorithm for reinforcement learning, the density threshold sequence is dynamically adjusted. The grid cells are divided into three levels: high density, medium density, and low density. The Adam optimizer is used for gradient descent, and the learning rate is set to 0.001. The grid size corresponding to different density levels is determined through iterative optimization. The grid size in the high-density area is 100 meters × 100 meters, the medium-density area is 250 meters × 250 meters, and the low-density area remains 500 meters × 500 meters.

[0090] Finally, a quadtree structure is used to achieve adaptive grid subdivision. The maximum depth of the quadtree is set to 4 layers. Starting from the root node, it is judged in turn whether the density value of each grid cell exceeds the corresponding density threshold. If it exceeds, the grid is divided into four sub-grids. For adjacent small grids with similar densities, if the density difference is less than the set threshold, they are merged. Through iterative splitting and merging operations, a grid-based urban area is finally generated.

[0091] Taking the central area of a certain city as an example, the area of this area is about 25 square kilometers. Through the above method, the high-density commercial area uses a fine grid of 100 meters × 100 meters, the residential area uses a medium grid of 250 meters × 250 meters, and the suburban green space uses a rough grid of 500 meters × 500 meters. The total number of grids is about 400, reducing the computational burden by about 40% compared with the unified grid division.

[0092] In this embodiment, an adaptive grid division of urban areas is achieved through the combination of deep learning and reinforcement learning, improving the accuracy of urban data analysis and management. The combination of the deep autoencoder and the graph attention network fully exploits the latent features and spatial correlations of urban activity data, enhancing the rationality of grid division; the adaptive density threshold learning model can dynamically adjust the grid size according to regional characteristics, ensuring both the refined representation of high-density areas and avoiding resource waste in low-density areas. The grid division method with a quadtree structure has good scalability and flexibility, can be dynamically adjusted with the development of the region; significantly reduces data storage and computing resource consumption, and improves urban management efficiency.

[0093] In an alternative embodiment, based on the sequence of regional density thresholds, grid splitting and merging are performed to generate a grid-based urban area, including:

[0094] When the grid optimization function value is greater than the corresponding regional density threshold and the spatial dependence relationship is less than the first preset threshold, the initial grid cell is non-uniformly quad-split. When the grid optimization function value of adjacent grid cells is less than the corresponding regional density threshold and the spatial dependence relationship is greater than the second preset threshold, the adjacent grid cells are adaptively merged;

[0095] The non-uniform quad-split includes: calculating the density gradient field within the initial grid cell, determining the main direction and secondary direction of density change based on the density gradient field; calculating the relative importance measure at each position within the initial grid cell to construct a weight matrix;

[0096] In the main direction, based on the cumulative distribution function of the weight matrix, the position of the first splitting line is determined such that the sum of the weighted density values on both sides of the first splitting line is equal; in the secondary direction, the position of the second splitting line is calculated for each of the two sub-regions divided by the first splitting line such that the weighted density values of the four sub-grids satisfy the balance constraint;

[0097] The initial grid cell is divided into four sub-grids, where the area of each sub-grid is inversely proportional to the average density, and the area ratio of adjacent sub-grids does not exceed a preset maximum ratio threshold;

[0098] The adaptive merge includes: calculating the shape similarity and density difference degree of adjacent grid cells. When the shape similarity is greater than the third preset threshold and the density difference degree is less than the fourth preset threshold, the adjacent grid cells are merged into a new grid cell, and the boundary of the new grid cell is obtained by least squares fitting of the boundary points of the adjacent grid cells.

[0099] In a specific embodiment, first, the grid optimization function value is calculated based on the population distribution data in the initial grid cell. By obtaining the census data, its normalization processing is carried out, and the population density of each grid cell is calculated. At the same time, the spatial dependence relationship between adjacent grid cells is calculated, specifically obtained by calculating the correlation coefficient of the population density of adjacent grid cells.

[0100] When the grid optimization function value is greater than the corresponding regional density threshold (for example, set to 0.8) and the spatial dependence relationship is less than the first preset threshold (for example, set to 0.3), non-uniform four-way splitting of the initial grid cell is required. First, the density gradient field within the initial grid is calculated. Based on the gray-scale image processing method, the first-order differences in the x-direction and y-direction are calculated to obtain the main direction and the secondary direction with the largest density change.

[0101] Then, the relative importance measure of each position within the grid is calculated to construct a weight matrix. The importance is measured by calculating the ratio of the mean value to the standard deviation of the population density in the surrounding area of each position. For example, the weight of the center position of a certain grid cell is 0.8, and the weight of the edge position is 0.2.

[0102] In the main direction, the position of the first splitting line is determined based on the cumulative distribution function of the weight matrix. Through iterative calculation, the sum of the weighted density values on both sides of the splitting line is made equal. For example, a certain grid cell divides the first splitting line at the position of x = 100, so that the weighted density values on both the left and right sides are 0.5.

[0103] In the secondary direction, the positions of the second splitting lines are calculated for the two sub-regions divided by the first splitting line respectively. Through iterative optimization, the weighted density values of the four sub-grids satisfy the balance constraint. For example, the weighted density values of the four sub-grids are 0.25, 0.25, 0.25, and 0.25 respectively.

[0104] Finally, the initial grid cell is divided into four sub-grids, where the area of each sub-grid is inversely proportional to its average density. For example, the area of a certain sub-grid is 100 square meters, and the corresponding average density is 0.01. The area ratio of adjacent sub-grids does not exceed the preset maximum ratio threshold (for example, set to 3).

[0105] When the grid optimization function value of adjacent grid cells is less than the regional density threshold, for example, set to 0.2, and the spatial dependence relationship is greater than the second preset threshold, for example, set to 0.7, adaptive merging is required. First, the shape similarity of adjacent grid cells is calculated, obtained by comparing the shape characteristics of the boundary curves.

[0106] Calculate the density difference degree simultaneously, which is obtained by calculating the difference coefficient of the population density of adjacent grid cells. When the shape similarity is greater than the third preset threshold, for example, set to 0.8, and the density difference degree is less than the fourth preset threshold, for example, set to 0.3, merge the adjacent grid cells into a new grid cell.

[0107] The boundary of the new grid cell is obtained by fitting the boundary points of adjacent grid cells using the least squares method. First, extract the coordinate of the feature points on the boundary of adjacent grid cells, and then obtain the new boundary curve equation through least squares fitting.

[0108] As Figure 2 shown in the grid deformation rate performance of different grid optimization methods under various density gradient conditions. This technical solution effectively controls the degree of grid deformation through the non-uniform four-split and adaptive merging strategies, especially showing outstanding performance in high-density gradient regions. Data shows that when the density gradient is 0, the grid deformation rate of this technical solution is only 2.1%, which has obvious advantages compared with 5.8% of the uniform four-split method and 3.3% of the traditional adaptive grid method. As the density gradient increases to 35 kg / m³ / m, the grid deformation rate of this technical solution is 22.7%, far lower than 37.5% of the uniform four-split method and 30.2% of the traditional adaptive grid method. Especially in the medium-high density gradient range of 20 - 25 kg / m³ / m, the gap between this technical solution and other methods is more significant. The grid deformation rates are 12.1% - 16.4% respectively, while the uniform four-split method reaches 26.3% - 32.8%, and the traditional adaptive grid method is 19.0% - 24.5%. This advantage stems from the strategy of determining the primary and secondary directions based on the density gradient field in this technical solution, as well as the method of determining the splitting line position using the cumulative distribution function of the weight matrix, making the grid division more in line with physical characteristics, and avoiding excessive deformation through the restriction of the area ratio of adjacent sub-grids (not exceeding the preset maximum ratio threshold), thereby improving the calculation accuracy while maintaining the grid quality.

[0109] In this embodiment, it is possible to achieve the adaptive splitting and merging of urban area grids, ensuring the rationality and effectiveness of grid division. Through the non-uniform four-split method, the heterogeneity of the density distribution within the region is considered, making the grid division more in line with the actual situation; the splitting method based on the weight matrix and density gradient field can accurately locate the primary and secondary directions of density changes, ensuring the balance of the density distribution of sub-grids after splitting. The adaptive merging strategy effectively avoids the problem of excessive grid subdivision; the least squares fitting method for the new grid boundary ensures the smoothness and continuity of the grid boundary, improving the visualization effect and practicality of the grid division result; the overall solution has strong universality and practical value.

[0110] In an alternative embodiment, multi-source spatio-temporal data of a gridded urban area is acquired, spatio-temporal alignment is performed according to a preset number of time segments, and a multi-dimensional time series feature set is generated; a multi-modal feature fusion algorithm is used to perform feature extraction and fusion operations on the multi-dimensional time series feature set to generate a multi-dimensional feature vector of grid cells, including:

[0111] The multi-source spatio-temporal data includes POI data, remote sensing image data, and population trajectory data of the gridded urban area;

[0112] Within a time segment, a POI status sequence is generated based on the POI data, a ground object distribution sequence is generated based on the remote sensing image data, and a population distribution sequence is generated based on the population trajectory data; the population distribution sequence is filtered for noise through Kalman filtering, and missing points are supplemented through piecewise cubic spline interpolation to generate population density time series data; the ground object distribution sequence is subjected to atmospheric correction and geometric correction to generate ground object coverage time series data; the POI status sequence, the ground object coverage time series data, and the population density time series data are aligned to a unified spatio-temporal reference system to generate spatio-temporally aligned data;

[0113] The POI status sequence is constructed into a multi-level POI semantic tree, the semantic information entropy of each layer in the multi-level POI semantic tree is calculated, and through recursive morphological decomposition, a semantic feature vector of grid cells is extracted layer by layer from the top layer to the bottom layer;

[0114] The ground object coverage time series data is constructed into a multi-scale pyramid, and through the bottom-up feature aggregation and top-down boundary optimization of the multi-scale pyramid, a scene feature vector of grid cells is extracted;

[0115] Based on the population density time series data, a time-varying kernel density is calculated, and through adaptive kernel density estimation, a flow feature vector of grid cells is extracted;

[0116] The semantic feature vector of grid cells, the scene feature vector of grid cells, and the flow feature vector of grid cells are combined to generate a multi-dimensional time series feature set, and standardized processing is performed;

[0117] A reconstruction error minimization model is constructed to calculate the feature fusion weights, and the standardized features are weighted and combined according to the feature fusion weights to generate a multi-dimensional feature vector of grid cells.

[0118] The POI data specifically refers to geographical data used to mark specific locations in a map or geographic information system, including location coordinates, names, addresses, categories, and other relevant information. In the identification of urban functional areas, POI data is used to analyze the distribution and status of different functional places within a region, reflect the semantic features of urban space, and assist in urban planning and management.

[0119] In a specific embodiment, POI data, remote sensing image data, and population trajectory data of a grid-based urban area are obtained. The POI data contains attribute information such as category, name, location, etc. The remote sensing image data contains multi-spectral and panchromatic bands. The population trajectory data contains timestamp and location information.

[0120] For the division of the time dimension of the grid cell, the 24-hour period is divided into 48 time segments at 30-minute intervals. Within each time segment, the three types of data are processed separately:

[0121] For the POI data, the POI open status (open or closed) within each time segment is extracted to form a POI status sequence. For the remote sensing image data, the multi-spectral band combination is used to extract the land cover types (buildings, roads, vegetation, etc.) to form a land cover distribution sequence. For the population trajectory data, the population quantity within the grid is counted to form a population distribution sequence.

[0122] For the population distribution sequence, the Kalman filter is applied to eliminate outliers, and the piecewise cubic spline interpolation function is used to supplement the missing data points to generate continuous population density time series data. For the land cover distribution sequence, the dark pixel method is used for atmospheric correction, and geometric correction is performed based on control point registration to generate normalized land cover time series data. The three types of time series data are unified under the same spatial reference system and time standard.

[0123] A 4-layer POI semantic tree is constructed, from the bottom layer to the top layer are specific POIs, sub-categories, main categories, and the overall. The normalized information entropy of each layer of nodes is calculated, and through recursive morphological operations, semantic features are extracted from the top layer to the bottom layer. Taking a certain area as an example, the top layer information entropy is 0.95, and the main category layer information entropy is 0.82, and finally a 128-dimensional semantic feature vector is obtained.

[0124] A 3-layer pyramid is constructed for the land cover data, the bottom layer resolution is 2 meters, and it is successively upsampled to form resolutions of 6 meters and 18 meters. Feature aggregation is performed from bottom to top through max pooling, and the boundary is optimized from top to bottom through conditional random fields to extract a 256-dimensional scene feature vector.

[0125] The time-varying kernel density of the population density is calculated, and the kernel function bandwidth is determined through cross-validation to extract a 64-dimensional flow feature vector. The three types of feature vectors are concatenated and z-score standardized.

[0126] A feature fusion model based on an autoencoder is constructed, and the feature weights are determined by minimizing the reconstruction error. The standardized features are weighted and summed to generate a 448-dimensional grid cell feature vector. Experiments show that the feature vector can effectively represent the spatio-temporal features of the grid cell.

[0127] In this embodiment, through the fusion and alignment of multi-source heterogeneous data, a comprehensive characterization of urban elements is achieved, improving the integrity and accuracy of feature expression. The collaborative processing of multi-modal data effectively reduces the limitations of a single data source; methods such as multi-level semantic trees, multi-scale pyramids, and adaptive kernel density estimation are used to extract features from dimensions such as semantics, scenarios, and population mobility, enhancing the feature expression ability. Recursive morphological decomposition and hierarchical feature extraction ensure the multi-scale nature of features; a feature fusion method based on minimizing reconstruction error adaptively determines feature weights, avoiding the subjectivity of manually setting parameters. Standardization processing and weighted combination improve the comparability and robustness of features. The overall solution has strong generality and scalability.

[0128] In an alternative embodiment, based on the multi-dimensional feature vector of grid cells, the extraction of long-range temporal features using a bidirectional gated recurrent unit includes:

[0129] The multi-dimensional feature vector of grid cells is sampled hierarchically according to daily, weekly, and monthly spans to obtain multi-scale temporal feature sequences; for each time-scale feature sequence in the multi-scale temporal feature sequences, feature extraction is performed through a bidirectional gated recurrent unit, and the bidirectional gated recurrent unit is coupled with a spiking neuron module, and the bidirectional gated recurrent unit includes a forward gated recurrent unit and a backward gated recurrent unit;

[0130] The update gate output of the forward gated recurrent unit is converted into a first pulse sequence through a first pulse encoder, and the first pulse encoder is based on a threshold integrate-and-fire model, generating a pulse when the accumulated input potential exceeds a preset firing threshold and resetting the membrane potential; the reset gate output of the forward gated recurrent unit is converted into a second pulse sequence through a second pulse encoder, and the second pulse encoder is the same as the first pulse encoder; a forward hidden state sequence is generated according to the firing times of the first pulse sequence and the second pulse sequence;

[0131] The update gate output and the reset gate output of the backward gated recurrent unit are respectively converted through the first pulse encoder and the second pulse encoder, and a backward hidden state sequence is generated according to the firing times of the converted pulse sequences;

[0132] The time interval between adjacent pulses in the forward hidden state sequence is calculated to obtain a first pulse firing intensity, and the time interval between adjacent pulses in the backward hidden state sequence is calculated to obtain a second pulse firing intensity; the first pulse firing intensity and the second pulse firing intensity are respectively converted into a first attention score and a second attention score, and normalized processing is performed to obtain a forward attention weight and a backward attention weight respectively;

[0133] Perform weighted calculations on the forward hidden state sequence and the backward hidden state sequence with the corresponding forward attention weights and backward attention weights respectively to obtain a fused hidden state sequence; perform temporal convolutional processing on the fused hidden state sequence to obtain local temporal features; perform a non-linear transformation on the grid cell multi-dimensional feature vector to obtain residual features; add the local temporal features and the residual features to obtain enhanced temporal features; perform normalization processing on the enhanced temporal features to obtain long-range temporal features.

[0134] In a specific embodiment, first perform multi-scale temporal sampling on the multi-dimensional feature vector of the grid cell to construct feature sequences with different time granularities. Divide according to the daily scale of 24 hours, the weekly scale of 7 days, and the monthly scale of 30 days. Taking a grid cell in a commercial area as an example, select features such as pedestrian flow, consumption index, and commuting index. Sample once per hour on the daily scale to form a 24-dimensional sequence, sample once per day on the weekly scale to form a 7-dimensional sequence, and sample once per day on the monthly scale to form a 30-dimensional sequence.

[0135] For the feature sequence of each time scale, process it through a bidirectional gated recurrent unit of an integrate-and-fire neuron module. The integrate-and-fire neuron module is a bionic simulation of the information processing mechanism of biological neurons, including three core components: membrane potential accumulation, threshold control, and synaptic plasticity. In the membrane potential accumulation unit, maintain a dynamic membrane potential with an initial value of 0. When receiving a feature input, multiply the feature value by the synaptic weight and accumulate it on the membrane potential, while considering the exponential decay characteristic of the membrane potential. For example, for an input feature value of 0.5, a synaptic weight of 1.2, a current membrane potential of 0.3, and a decay coefficient of 0.9, the updated membrane potential value is (0.5×1.2 + 0.3)×0.9 = 0.81.

[0136] The threshold control unit sets the firing threshold to 0.8 and the refractory period to 2 time steps. When the membrane potential exceeds the firing threshold and the time interval since the last update is greater than the refractory period, trigger feature update and reset the membrane potential to 0. The synaptic plasticity unit adopts a time-dependent plasticity rule to dynamically adjust the synaptic weight according to the time relationship of feature updates of adjacent neurons. When the feature update of the latter neuron occurs within a short time after the feature update of the previous neuron, such as a time interval of 3 time steps, increase the synaptic weight by 10%; when the time interval exceeds 10 time steps, decrease the synaptic weight by 5%.

[0137] The bidirectional gated recurrent unit contains two processing branches, namely the forward branch and the backward branch. In the forward branch, the output of the update gate is transformed through the first spiking neuron module. When the accumulated eigenvalue exceeds the threshold, feature update is triggered and temporal coding is generated. The output of the reset gate is processed through the second spiking neuron module with the same configuration. Based on the feature update sequences of the two modules, a forward hidden state sequence is generated. The backward branch adopts the same structure to generate a backward hidden state sequence.

[0138] Taking the morning rush hour data of the commercial area grid cells as an example, when the pedestrian flow increases rapidly, the membrane potential accumulates faster, triggering frequent feature updates; during the stable period of the pedestrian flow, the feature update frequency decreases. Through this mechanism, the key change moments in the data can be adaptively captured. For the generated hidden state sequence, calculate the time intervals between adjacent feature updates. For example, if the forward hidden state sequence is updated at times 2, 5, and 8, the corresponding time interval sequence is 3, 3. Convert the time interval sequence into attention scores and normalize them, and the possible weights obtained may be 0.4 and 0.6.

[0139] Weightedly combine the forward and backward hidden state sequences with the corresponding attention weights respectively to obtain a fused hidden state sequence. Then perform temporal convolution on the fused sequence using a convolutional kernel of size 3 to extract local temporal feature patterns. At the same time, obtain residual features by non-linearly transforming the original feature vector, and add them element-wise to the local temporal features to achieve feature enhancement. Finally, through normalization processing, normalize the eigenvalues to a zero-mean and unit-variance distribution to obtain the final long-range temporal features.

[0140] The processed long-range temporal features can effectively characterize the activity patterns of the commercial area at different time scales. Larger values in the feature vector correspond to the typical activity periods of the area, such as the commuting patterns during the morning and evening rush hours on weekdays and the consumption peaks on weekends. These features provide a reliable data basis for subsequent regional function recognition and analysis. The entire feature extraction process fully draws on the information processing mechanism of the biological nervous system and has strong adaptability and practical value.

[0141] In this embodiment, through the coupled design of spiking neurons and GRU, the temporal characteristics of spiking neurons are fully utilized to enhance the modeling ability of long-term temporal dependence relationships; a multi-scale hierarchical sampling strategy is adopted to capture the feature change laws at different time scales, enhancing the model's perception ability of different periodic patterns; an attention mechanism is introduced to model the spiking intensity, combined with the residual connection design, effectively alleviating the gradient vanishing problem in long sequence processing and improving the robustness of feature extraction.

[0142] In an alternative embodiment, the local density value and the relative distance value of the enhanced spatio-temporal features are calculated, the density peak clustering is performed to determine the clustering center points of the functional areas, and the initial functional area division result is generated, including:

[0143] Calculate the Euclidean distance, Mahalanobis distance and cosine distance for the enhanced spatio-temporal features, determine the weighting coefficient through the variance contribution rate of the enhanced spatio-temporal features, and perform weighted fusion to obtain the comprehensive distance matrix;

[0144] Sort the comprehensive distance matrix in ascending order, calculate the kernel density estimation curve of the comprehensive distance matrix, obtain all local minimum points of the kernel density estimation curve, and select the first distance value in the kernel density estimation curve that satisfies that the derivative of the local minimum point is negative and the amplitude is less than the mean value of the kernel density estimation curve as the truncation distance;

[0145] Calculate the distribution density of the comprehensive distance matrix, determine the bandwidth parameter of the adaptive Gaussian kernel function, construct the adaptive Gaussian kernel function using the bandwidth parameter, and calculate the local density value of each grid cell based on the truncation distance and the adaptive Gaussian kernel function; Calculate the initial relative distance value of each grid cell based on the local density value, and calculate the minimum distance between the grid cell and the grid cell with a higher local density value under multiple preset search radii, and perform weighted fusion based on the stability of the initial relative distance value under different search radii to determine the relative distance value;

[0146] Take the product of the local density value and the relative distance value as the clustering center score, sort the clustering center scores in descending order, select the grid cells in the top 30% of the sorting and satisfying the spatial distribution uniformity as the clustering center points of the functional areas, and at the same time determine the clustering center points of the non-functional areas;

[0147] Sort all grid cells in descending order according to the local density value. For each grid cell of the clustering center point of the non-functional area, assign the grid cell to the nearest clustering center point of the functional area according to the density reachability of the local density value, and generate the initial functional area division result.

[0148] In a specific implementation, in order to analyze the spatio-temporal characteristics of grid cells and divide functional areas, it is first necessary to calculate the comprehensive distance between grid cells. Specifically, the Euclidean distance is used to calculate the distance relationship of grid cells in physical space, the Mahalanobis distance is used to calculate the standardized distance of grid cells in the feature space, and the cosine distance is used to calculate the included angle of the feature vectors of grid cells. Based on the variance contribution rate of the spatio-temporal characteristics of each grid cell, the weight coefficients of different distance metrics are determined, and the three distances are weighted and fused to obtain a comprehensive distance matrix. For example, if the Euclidean distance of a certain grid cell is 0.3, the Mahalanobis distance is 0.5, and the cosine distance is 0.4, and the corresponding weights are 0.4, 0.3, and 0.3 respectively, then the fused comprehensive distance is 0.39.

[0149] Then, the comprehensive distance matrix is sorted in ascending order, and the kernel density estimation method is used to draw the distance distribution curve. By calculating the first derivative of the curve, all local minimum points are obtained, and the distance that satisfies the derivative being negative and the amplitude being less than the mean value of the curve is selected as the truncation distance. For example, if the mean value of a certain distance distribution curve is 0.6, and the distance corresponding to the first local minimum point that satisfies the condition is 0.45, then 0.45 is used as the truncation distance.

[0150] Subsequently, the distribution density of the comprehensive distance matrix is calculated to determine the bandwidth parameter of the adaptive Gaussian kernel function. The adaptive Gaussian kernel function is constructed using the bandwidth parameter, and the local density value of each grid cell is calculated in combination with the truncation distance. Then, the initial relative distance value is calculated based on the local density value. At multiple preset search radii (such as 0.1, 0.2, 0.3), the minimum distance between the grid cell and the grid cell with a high density value is calculated respectively. According to the stability of the relative distance values under different search radii, the weights are determined, and the weighted fusion is performed to obtain the final relative distance value.

[0151] The product of the local density value and the relative distance value is used as the clustering center score, and the scores are sorted in descending order. The grid cells that rank in the top 30% and have a uniform spatial distribution (the distance between adjacent grid cells is greater than the preset threshold) are selected as the functional area clustering center points, and the rest are non-clustering center points. For example, if there are 100 grid cells in a certain area, the grid cells that rank in the top 30 and have a spacing greater than 0.5 are selected as the clustering center points.

[0152] Finally, all grid cells are sorted in descending order of local density value. For each grid cell that is not a clustering center point, based on the continuity of the density value (the difference in density values between adjacent grid cells is less than the threshold), it is assigned to the nearest clustering center point to form the initial functional area division result.

[0153] In this embodiment, by integrating multiple distance measurement methods and comprehensively considering the similarity of grid cells in the physical space and the feature space, the accuracy and robustness of distance calculation are improved; the use of an adaptive Gaussian kernel function and a multi-scale search strategy to calculate the local density value and the relative distance value can better adapt to the data distribution characteristics of different regions and enhance the accuracy of density peak clustering; by selecting the functional area clustering center point through clustering center scoring and spatial uniformity constraints and combining the density reachability principle for grid cell allocation, a reasonable functional area division is achieved, ensuring the spatial continuity and integrity of the division result.

[0154] In an alternative embodiment, an affinity matrix is constructed using enhanced spatio-temporal features, and spectral clustering operations are performed in combination with spatial adjacency constraints to generate an optimized functional area division result; through a simulated annealing boundary optimization algorithm, the indicators of regional internal consistency and inter-regional difference are iteratively optimized to generate the final urban functional area division result, including:

[0155] The similarity value between the enhanced spatio-temporal features is calculated using a radial basis function to construct a similarity matrix; the adjacency strength between grid cells is determined based on the spatial position relationship of the grid cells to construct a spatial adjacency constraint matrix;

[0156] The weight coefficient is determined according to the variance contribution rate of the enhanced spatio-temporal features and the spatial adjacency constraint matrix, and the similarity matrix and the spatial adjacency constraint matrix are weighted and combined to obtain an affinity matrix; the degree matrix of the affinity matrix is calculated, and the affinity matrix is Laplace-normalized using the degree matrix to obtain a normalized affinity matrix; the normalized affinity matrix is eigen-decomposed to obtain an eigenvalue sequence, and based on the descending order of the eigenvalue sequence, the number of eigenvectors is determined, and the corresponding eigenvectors are selected to construct an eigenmatrix; the row vectors of the eigenmatrix are normalized to obtain low-dimensional representation features, and in the feature space composed of the low-dimensional representation features, clustering operations are performed to reallocate grid cells to different functional areas to generate an optimized functional area division result;

[0157] The variance of the enhanced spatio-temporal features of the grid cells within each functional area in the optimized functional area division result is calculated to obtain the regional internal consistency index value; the difference degree of the enhanced spatio-temporal features between adjacent functional areas in the optimized functional area division result is calculated to obtain the inter-regional difference index value; the regional internal consistency index value and the inter-regional difference index value are combined to construct a simulated annealing optimization objective function;

[0158] Based on a preset initial temperature and a cooling coefficient, start the simulated annealing operation. Select boundary grid cells from the results of the optimized functional area division, randomly change the functional area category of the boundary grid cells, determine the changed plan, calculate the function values of the simulated annealing optimization objective function before and after the change, determine the difference. When the difference is less than zero, accept the changed plan. When the difference is greater than zero, calculate the acceptance probability based on the cooling coefficient and the current temperature, and at the same time determine the new temperature. Based on the acceptance probability, determine whether to accept the changed plan; repeat the execution until the new temperature reaches the preset stop temperature, and generate the final urban functional area division result.

[0159] In a specific implementation, first calculate the similarity based on enhanced spatio-temporal features, and calculate the similarity of the feature vectors of each grid cell through a radial basis function. For example, set the parameter to 0.5, substitute the Euclidean distance between grid cells into the radial basis function, and obtain a similarity matrix with similarity values between 0 and 1. At the same time, calculate the adjacency degree according to the relative position relationship of the grid cells in space. The adjacency strength of adjacent grid cells is 1, and that of non-adjacent grid cells is 0, and an adjacency constraint matrix is constructed.

[0160] Calculate the variance contribution rate of the similarity matrix and the adjacency constraint matrix to the classification of grid cells respectively, and use the variance contribution rate as the weight coefficient. For example, the variance contribution rate of the similarity matrix is 0.7, and that of the adjacency constraint matrix is 0.3. Based on this, the two matrices are weighted and combined to obtain an affinity matrix. Subsequently, calculate the sum of the elements in each row of the affinity matrix to form a degree matrix, and use the degree matrix to perform Laplacian normalization on the affinity matrix.

[0161] Perform eigenvalue decomposition on the standardized affinity matrix to obtain a sequence of eigenvalues arranged in descending order. By observing the change trend of the eigenvalues, truncate at the point where the eigenvalues show a significant decrease, and select the corresponding eigenvectors to construct an eigenmatrix. For example, select the first 5 eigenvectors. Perform row vector normalization on the eigenmatrix to obtain low-dimensional representation features. Perform K-means clustering in this feature space to reassign grid cells to different functional areas.

[0162] Calculate the variance of the features of the grid cells within each functional area to obtain an intra-region consistency index. Calculate the distance of the features of the grid cells between adjacent functional areas to obtain an inter-region difference index. Combine these two indices by weighting to construct an optimization objective function. Set the initial temperature to 100 and the cooling coefficient to 0.95.

[0163] During the simulated annealing process, randomly select boundary grid cells, change their functional area categories, and calculate the objective function values before and after the change. When the function value decreases, directly accept the new plan; when the function value increases, calculate the acceptance probability based on the current temperature. The temperature gradually decreases until it reaches the stop temperature of 1, and the final functional area division result is output.

[0164] AsFigure 3 As shown, it presents the spatial distribution results of the urban functional area division generated by this technical solution, taking a partial area of a central urban area as an example. Five typical functional areas are clearly presented in the figure: commercial area, residential area, industrial area, public service area, and green space and leisure area. Each functional area marks the grid cells with unique symbols. It can be seen from the figure that the functional areas generated by this technical solution have obvious spatial aggregation and boundary continuity. The internal consistency indicators of each functional area are as follows: commercial area 0.094, residential area 0.108, industrial area 0.112, public service area 0.103, and green space and leisure area 0.087, all of which are lower than 0.12, indicating that the spatio-temporal characteristics of the grid cells within each functional area are highly similar. At the same time, the difference indicators between adjacent functional areas are between 0.684 and 0.827, with an average of 0.741, indicating that the boundaries between different functional areas are clear and the characteristic differences are significant. From the perspective of spatial distribution, the commercial area is mainly located in the center of the city (near 116.30°E, 39.96°N), with an area of 25.4 square kilometers, adjacent to the residential area and the public service area; the residential area is located in the southwest of the commercial area, with an area of 32.7 square kilometers, being the largest functional area; the industrial area is located in the east of the residential area, with an area of 30.2 square kilometers, adjacent to the green space and leisure area; the public service area is located in the east of the city (near 116.45°E), with an area of 21.8 square kilometers; the green space and leisure area is located in the south of the city (near 39.85°N), with an area of 18.9 square kilometers, being the functional area with the highest internal consistency. This distribution pattern of functional areas conforms to the functional zoning principle in urban planning theory, reflecting the advantages of this technical solution in balancing spatial continuity and functional differences.

[0165] In this embodiment, through the dual constraints of feature similarity and spatial adjacency constraints, the spatial continuity of the functional area division and the consistency of regional characteristics are ensured; by using the dimensionality reduction characteristics of spectral clustering, the essential characteristics of the grid cells are effectively extracted, noise interference is reduced, and the accuracy and stability of the clustering effect are improved; the simulated annealing algorithm is used to optimize the functional area boundaries, enhancing the differences between regions while ensuring the internal consistency of the regions, making the functional area division results more reasonable and reliable.

[0166] In an alternative embodiment, a high-tech industrial development zone of a certain city is selected as the research object, and the total area of the research region is about 100 square kilometers. This method uses a regular grid of 250 meters × 250 meters as the basic analysis unit, integrates multi-source data such as mobile communication signaling data, high-resolution remote sensing images, point of interest data, and open road network data, and realizes refined functional area identification.

[0167] In the data preprocessing stage, first, semantic segmentation is performed on high-resolution remote sensing images. Using a deep convolutional neural network model, the study area is divided into basic land cover types such as buildings, water bodies, grasslands, forests, bare soil, and farmland. For each 250-meter grid, the proportion of pixels of each type is calculated. When the proportion of a certain type exceeds 85%, the grid is directly labeled as the corresponding natural land use type. For the building area, by calculating the length of the building projection shadow and combining the shooting time and the solar altitude angle, the building height is estimated. The specific calculation process is as follows: First, the building outline and shadow boundary are extracted, the shadow length L is calculated, and then the building height H is calculated according to the formula H = L × tan(solar altitude angle). When the average building height in the grid is greater than 50 meters, it is preliminarily determined as a business and office area; when the building height is low and the floor area is large, it is preliminarily determined as an industrial area.

[0168] In terms of road network analysis, open road network data is used, and different road widths are set according to the road grade: 25 meters for main roads, 15 meters for secondary roads, and 8 meters for branch roads. The road occupancy range is determined through buffer analysis, and the proportion of road area in each grid is calculated. When the proportion of road area exceeds 30%, the grid is labeled as road traffic land. For grids that contain both roads and other land use types, they are recorded in a multi-label form, such as "road + grassland", "road + commercial", etc.

[0169] In terms of extracting human flow characteristics, mobile communication signaling data for three consecutive weeks is used, divided into two situations: weekdays and weekends. For base station data, first, kernel density estimation is performed: with each base station as the center, a Gaussian kernel function is set, and the standard deviation of the kernel function is set to 1 / 3 of the base station coverage radius. Through superposition calculation, the human flow density distribution of 250-meter grids is obtained. For the directly obtained grid data, time mean calculation is performed: the mean value of 15-day data is taken for weekdays, and the mean value of 6-day data is taken for weekends. Finally, the time series curve of the human flow volume of each grid for 24 hours on weekdays and weekends is obtained.

[0170] To identify the temporal characteristics of different functional areas, this method innovatively uses a self-organizing mapping neural network for clustering analysis. The network consists of two layers: an input layer and a competitive layer. The number of nodes in the input layer is 48 (24 hours on weekdays + 24 hours on weekends), and the competitive layer is set as a 6×6 two-dimensional node array. The network training process includes two steps: competitive learning and weight update. First, calculate the Euclidean distance between the input vector and each node in the competitive layer, and the node with the minimum distance is the winning node. Then, update the weight vectors of the winning node and its neighborhood nodes to make them closer to the input vector. Through repeated iteration, finally, the grids with similar temporal characteristics are mapped to adjacent nodes in the competitive layer. After clustering analysis, typical temporal patterns of functional areas are identified: the working area shows a double peak on weekdays and a trough on weekends; the commercial area shows a continuous peak during the day and a significant increase on weekends; the residential area shows an inverse double peak in the morning and evening and is stable on weekends; the transportation hub area shows multiple peaks with relatively weak regularity.

[0171] In terms of point-of-interest data analysis, this method uses the term frequency-inverse document frequency algorithm for processing. Each grid is regarded as a document, and the types of points of interest within the grid are regarded as terms, and the TFIDF values of each type of point of interest are calculated. The calculation process is as follows: First, count the frequency of occurrence of each type within the grid and divide it by the total number of points of interest within the grid to obtain the term frequency TF. Then, calculate the proportion of the number of grids containing this type in the total number of grids and take the logarithm of its reciprocal to obtain the inverse document frequency IDF. The two are multiplied to obtain the TFIDF value, which is used to characterize the indication of this type of point of interest to the grid function. Through this method, the dominant functional type of each grid can be identified, especially for the subdivision of the working area: when the TFIDF value of administrative and office-related points of interest is the highest, it is determined as an administrative and office area; when the TFIDF value of education and scientific research-related points of interest is the highest, it is determined as an education and scientific research area; when the TFIDF value of medical and health-related points of interest is the highest, it is determined as a medical area.

[0172] Finally, this method uses a multi-label classification method to output the functional area recognition results. For each grid, considering the land use type, building characteristics, temporal pattern of pedestrian flow, and distribution characteristics of points of interest, its main functional type and secondary functional type are determined. For example, a certain grid may have the attributes of "business office + commerce", indicating that there are both office buildings and commercial facilities in this area; a certain grid may have the attributes of "residence + education", indicating that there are both residential buildings and schools in this area. This multi-label classification method is more in line with the composite characteristics of urban functions and can more accurately describe the actual usage of the area.

[0173] Through verification by application in the test area, this method shows good recognition effects. For areas dominated by a single function, such as pure residential areas and pure commercial areas, the recognition accuracy rate exceeds 90%; for areas with mixed functions, the recognition accuracy rate of the dominant function exceeds 85%, and the recognition accuracy rate of the secondary function exceeds 75%. The advantage of this method lies in making full use of the complementarity of multi-source data, considering both spatial morphological features and temporal activity features, and being able to comprehensively depict the spatial distribution law of regional functions.

[0174] Figure 4 FIG. is a schematic structural diagram of a grid-based recognition and classification system for urban functional areas based on multi-source spatio-temporal data according to an embodiment of the present invention, as Figure 4 shown, the system includes:

[0175] A first unit, configured to divide an urban area into multiple grid cells by using an adaptive grid division algorithm, dynamically adjust the size of the grid cells based on a regional density threshold to generate a grid-based urban area; obtain multi-source spatio-temporal data of the grid-based urban area, perform spatio-temporal alignment according to a preset number of time segments to generate a multi-dimensional time series feature set; perform feature extraction and fusion operations on the multi-dimensional time series feature set through a multi-modal feature fusion algorithm to generate a multi-dimensional feature vector of the grid cell;

[0176] A second unit, configured to extract long-range time series features by using a bidirectional gated recurrent unit based on the multi-dimensional feature vector of the grid cell; calculate multi-scale spatial correlation coefficients between the grid cell and surrounding grid cells based on a multi-head attention mechanism and perform feature aggregation to generate spatial correlation features; input the long-range time series features and the spatial correlation features into a feature fusion network with skip connections to generate enhanced spatio-temporal features;

[0177] A third unit, configured to calculate the local density value and relative distance value of the enhanced spatio-temporal features, perform density peak clustering to determine the clustering center points of the functional areas to generate an initial functional area division result; construct an affinity matrix by using the enhanced spatio-temporal features, perform spectral clustering operations in combination with spatial adjacency constraints to generate an optimized functional area division result; through a simulated annealing boundary optimization algorithm, iteratively optimize the regional internal consistency and inter-regional difference indexes to generate a final urban functional area division result.

[0178] In the third aspect of the embodiment of the present invention,

[0179] A kind of electronic device is provided, including:

[0180] A processor;

[0181] A memory for storing instructions executable by the processor;

[0182] Wherein, the processor is configured to call the instructions stored in the memory to execute the method described above.

[0183] In the fourth aspect of the embodiments of the present invention,

[0184] a computer-readable storage medium is provided, on which computer program instructions are stored, and when the computer program instructions are executed by a processor, the foregoing method is implemented.

[0185] The present invention may be a method, an apparatus, a system, and / or a computer program product. The computer program product may include a computer-readable storage medium, on which computer-readable program instructions for executing various aspects of the present invention are loaded.

[0186] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit them. Although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that they can still modify the technical solutions described in the foregoing embodiments, or perform equivalent replacements for some or all of the technical features. However, these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A grid-based identification and classification method for urban functional areas based on multi-source spatiotemporal data, characterized in that: include: An adaptive grid partitioning algorithm is used to divide the urban area into multiple grid cells, and the size of the grid cells is dynamically adjusted based on the regional density threshold to generate a gridded urban area; wherein, urban data is obtained and preprocessed, the grid optimization function value and spatial dependency of the initial grid cells are calculated, the regional density threshold sequence is generated, and quadtree adaptive subdivision is performed; When the grid optimization function value is greater than the corresponding regional density threshold and the spatial dependency is less than the first preset threshold, the initial grid unit is non-uniformly split into four, the density gradient field in the initial grid unit is calculated, and the main direction and secondary direction of the density change are determined; the relative importance measure of each position in the initial grid unit is calculated, and a weight matrix is ​​constructed; in the main direction, the position of the first splitting line is determined based on the cumulative distribution function of the weight matrix; in the secondary direction, the position of the second splitting line is calculated for the two sub-regions divided by the first splitting line respectively; the initial grid unit is divided into four sub-grids; When the grid optimization function value of adjacent grid cells is less than the corresponding regional density threshold and the spatial dependence is greater than the second preset threshold, the adjacent grid cells are adaptively merged, and the shape similarity and density difference of the adjacent grid cells are calculated. When the shape similarity is greater than the third preset threshold and the density difference is less than the fourth preset threshold, the adjacent grid cells are merged into a new grid cell, and the boundary of the new grid cell is obtained by fitting the boundary points of the adjacent grid cells by the least squares method; Acquire multi-source spatiotemporal data of gridded urban areas, perform spatiotemporal alignment according to the preset number of time segments, and generate a multidimensional time series feature set; perform feature extraction and fusion operations on the multidimensional time series feature set through a multimodal feature fusion algorithm to generate a multidimensional feature vector of the grid unit; Based on the multi-dimensional feature vector of the grid unit, a bidirectional gated recurrent unit is used to extract long-range temporal features. Based on the multi-head attention mechanism, the multi-scale spatial correlation coefficients between the grid unit and the surrounding grids are calculated and feature aggregation is performed to generate spatial correlation features. The long-range temporal features and spatial correlation features are input into a feature fusion network with skip connections to generate enhanced spatiotemporal features. The local density value and relative distance value of the enhanced spatiotemporal features are calculated, and the density peak clustering is performed to determine the center point of the functional area cluster to generate the initial functional area division result; the affinity matrix is ​​constructed using the enhanced spatiotemporal features, and the spectral clustering operation is performed in combination with the spatial adjacency constraints to generate the optimized functional area division result; through the simulated annealing boundary optimization algorithm, the internal consistency of the region and the difference between regions are iteratively optimized to generate the final urban functional area division result.

2. The method according to claim 1, characterized in that The adaptive grid division algorithm is used to divide the urban area into multiple grid cells. The size of the grid cells is dynamically adjusted based on the regional density threshold. The generated gridded urban area includes: Acquire urban data of the target urban area, the urban data including pedestrian density data, building density data, land use type data, mobile phone signal data and environmental perception data, perform data cleaning and normalization preprocessing on the urban data, and generate an urban activity feature matrix; The target urban area is divided into initial grid units of a preset maximum size, a multidimensional density assessment model is constructed based on the urban activity feature matrix, a deep autoencoder is used to extract the hidden layer features of each initial grid unit, and the regional density distribution of each initial grid unit is calculated in combination with variational inference. The regional density distribution is input into a graph attention network, the spatial dependency relationship between adjacent initial grid units is calculated, and a grid optimization function value of each initial grid unit is generated; An adaptive density threshold learning model is constructed based on the grid optimization function value, a recursive neural network is used to model the time series change of the grid optimization function value, and a regional density threshold sequence is dynamically adjusted in combination with reinforcement learning. The initial grid unit is divided into a plurality of density levels according to the regional density threshold sequence. For each density level, a gradient descent method is used to optimize and determine the grid size, and the grid size is dynamically adjusted within a preset maximum grid size and a minimum grid size. The initial grid unit is adaptively subdivided through a quadtree structure, and grid splitting and merging are performed based on the regional density threshold sequence to generate a gridded urban area.

3. The method according to claim 1, characterized in that Obtain multi-source spatiotemporal data of gridded urban areas, perform spatiotemporal alignment according to a preset number of time segments, and generate a multidimensional time series feature set; The multi-dimensional time series feature set is extracted and fused by the multi-modal feature fusion algorithm to generate the multi-dimensional feature vector of the grid unit, including: The multi-source spatiotemporal data includes POI data, remote sensing image data and population trajectory data of gridded urban areas; In a time segment, a POI state sequence is generated based on the POI data, a land object distribution sequence is generated based on the remote sensing image data, and a population distribution sequence is generated based on the population trajectory data; the population distribution sequence is subjected to noise filtering by Kalman filtering, and missing points are supplemented by piecewise cubic spline interpolation to generate population density time series data; atmospheric correction and geometric correction are performed on the land object distribution sequence to generate land object coverage time series data; the POI state sequence, the land object coverage time series data and the population density time series data are aligned to a unified spatiotemporal reference system to generate spatiotemporal alignment data; The POI state sequence is constructed into a multi-level POI semantic tree, the semantic information entropy of each level in the multi-level POI semantic tree is calculated, and the semantic feature vectors of grid units are extracted layer by layer from the top level to the bottom level through recursive morphological decomposition; The ground feature coverage time series data is constructed into a multi-scale pyramid, and a grid unit scene feature vector is extracted through bottom-up feature aggregation and top-down boundary optimization of the multi-scale pyramid; Calculating time-varying kernel density based on the population density time series data, and extracting grid unit flow feature vectors through adaptive kernel density estimation; Combining the grid unit semantic feature vector, the grid unit scene feature vector and the grid unit flow feature vector to generate a multi-dimensional time series feature set, and performing standardization processing; A reconstruction error minimization model is constructed to calculate feature fusion weights, and the standardized features are weighted combined according to the feature fusion weights to generate a multi-dimensional feature vector of the grid unit.

4. The method according to claim 1, characterized in that: Based on the multi-dimensional feature vector of the grid unit, a bidirectional gated recursive unit is used to extract long-range temporal features including: The multidimensional feature vector of the grid unit is sampled in layers according to the daily span, weekly span and monthly span to obtain a multi-scale time series feature sequence; the feature sequence of each time scale in the multi-scale time series feature sequence is extracted by a bidirectional gated recursive unit, the bidirectional gated recursive unit is coupled with a pulse neuron module, and the bidirectional gated recursive unit includes a forward gated recursive unit and a backward gated recursive unit; The update gate output of the forward gated recursive unit is converted into a first pulse sequence through a first pulse encoder, and the first pulse encoder is based on a threshold integral discharge model, and generates a pulse when the accumulated input potential exceeds a preset discharge threshold, and resets the membrane potential; the reset gate output of the forward gated recursive unit is converted into a second pulse sequence through a second pulse encoder, and the second pulse encoder is the same as the first pulse encoder; a forward hidden state sequence is generated according to the discharge moments of the first pulse sequence and the second pulse sequence; The update gate output and the reset gate output of the backward gated recursive unit are converted by the first pulse encoder and the second pulse encoder respectively, and a backward hidden state sequence is generated according to the discharge timing of the converted pulse sequence; Calculating the time interval between adjacent pulses in the forward hidden state sequence to obtain a first pulse emission intensity, and calculating the time interval between adjacent pulses in the backward hidden state sequence to obtain a second pulse emission intensity; converting the first pulse emission intensity and the second pulse emission intensity into a first attention score and a second attention score, respectively, and performing normalization processing to obtain a forward attention weight and a backward attention weight, respectively; The forward hidden state sequence and the backward hidden state sequence are weightedly calculated with the corresponding forward attention weight and the backward attention weight respectively to obtain a fused hidden state sequence; the fused hidden state sequence is subjected to temporal convolution processing to obtain local temporal features; the multidimensional feature vector of the grid unit is nonlinearly transformed to obtain residual features; the local temporal features and the residual features are added to obtain enhanced temporal features; the enhanced temporal features are normalized to obtain long-range temporal features.

5. The method according to claim 1, characterized in that Calculate the local density value and relative distance value of the enhanced spatiotemporal features, perform density peak clustering to determine the center point of the functional area cluster, and generate the initial functional area division results including: Calculating the Euclidean distance, Mahalanobis distance and cosine distance of the enhanced spatiotemporal features, determining weighting coefficients according to the variance contribution rate of the enhanced spatiotemporal features, and performing weighted fusion to obtain a comprehensive distance matrix; The comprehensive distance matrix is ​​sorted in ascending order, a kernel density estimation curve of the comprehensive distance matrix is ​​calculated, all local minimum points of the kernel density estimation curve are obtained, and the first distance value in the kernel density estimation curve that satisfies the local minimum point derivative is negative and the amplitude is less than the mean of the kernel density estimation curve is selected as the cutoff distance; Calculate the distribution density of the comprehensive distance matrix, determine the bandwidth parameter of the adaptive Gaussian kernel function, construct the adaptive Gaussian kernel function using the bandwidth parameter, calculate the local density value of each grid cell based on the cutoff distance and the adaptive Gaussian kernel function; calculate the initial relative distance value of each grid cell based on the local density value, calculate the minimum distance between the grid cell and the grid cell with a higher local density value under multiple preset search radii, and perform weighted fusion based on the stability of the initial relative distance value under different search radii to determine the relative distance value; The product of the local density value and the relative distance value is used as the cluster center score, the cluster center scores are sorted in descending order, and the grid cells in the top 30% of the sorting that meet the spatial distribution uniformity are selected as the functional area cluster center points, and the non-functional area cluster center points are determined at the same time; All grid cells are sorted in descending order according to the local density values. For each grid cell of the non-functional area cluster center point, the grid cell is assigned to the nearest functional area cluster center point according to the density accessibility of the local density value to generate an initial functional area division result.

6. The method according to claim 1, characterized in that The affinity matrix is ​​constructed by using enhanced spatiotemporal features, and the spectral clustering operation is performed in combination with spatial adjacency constraints to generate the optimized functional area division results; Through the simulated annealing boundary optimization algorithm, the internal consistency of the region and the difference between regions are iteratively optimized to generate the final urban functional area division results including: The similarity values ​​between the enhanced spatiotemporal features are calculated by radial basis function to construct a similarity matrix; the adjacency strength between the grid cells is determined according to the spatial position relationship of the grid cells to construct a spatial adjacency constraint matrix; Determine the weight coefficient according to the variance contribution rate of the enhanced spatiotemporal features and the spatial adjacency constraint matrix, perform weighted combination of the similarity matrix and the spatial adjacency constraint matrix to obtain an affinity matrix; calculate the degree matrix of the affinity matrix, use the degree matrix to perform Laplace normalization on the affinity matrix to obtain a standardized affinity matrix; perform eigendecomposition on the standardized affinity matrix to obtain an eigenvalue sequence, determine the number of eigenvectors based on the descending order of the eigenvalue sequence, and select the corresponding eigenvectors to construct a feature matrix; perform normalization on the row vectors of the feature matrix to obtain low-dimensional representation features, perform clustering operations in the feature space formed by the low-dimensional representation features, reallocate grid units to different functional areas, and generate optimized functional area division results; Calculate the enhanced spatiotemporal feature variance of the grid cells within each functional area in the optimized functional area division result to obtain a consistency index value within the region; calculate the enhanced spatiotemporal feature difference between adjacent functional areas in the optimized functional area division result to obtain a difference index value between regions; combine the consistency index value within the region with the difference index value between regions to construct a simulated annealing optimization objective function; Based on the preset initial temperature and cooling coefficient, the simulated annealing operation is started, the boundary grid unit is selected from the optimized functional area division result, the functional area category of the boundary grid unit is randomly changed, the changed scheme is determined, the function value of the simulated annealing optimization objective function before and after the change is calculated, and the difference is determined. When the difference is less than zero, the changed scheme is accepted. When the difference is greater than zero, the acceptance probability is calculated based on the cooling coefficient and the current temperature, and the new temperature is determined at the same time. Whether to accept the changed scheme is determined based on the acceptance probability; the execution is repeated until the new temperature reaches the preset stop temperature to generate the final urban functional area division result.

7. A grid-based identification and classification system for urban functional areas based on multi-source spatiotemporal data, used to implement the method described in any one of claims 1 to 6, characterized in that: include: The first unit is used to divide the urban area into multiple grid cells by using an adaptive grid division algorithm, dynamically adjust the size of the grid cells based on the regional density threshold, and generate a gridded urban area; obtain multi-source spatiotemporal data of the gridded urban area, perform spatiotemporal alignment according to a preset number of time segments, and generate a multi-dimensional time series feature set; perform feature extraction and fusion operations on the multi-dimensional time series feature set by a multi-modal feature fusion algorithm to generate a multi-dimensional feature vector of the grid cell; The second unit is used to extract long-range temporal features using a bidirectional gated recursive unit based on the multidimensional feature vector of the grid unit; Based on the multi-head attention mechanism, the multi-scale spatial correlation coefficients between the grid unit and the surrounding grids are calculated and feature aggregation is performed to generate spatial correlation features; the long-range temporal features and spatial correlation features are input into the feature fusion network with skip connections to generate enhanced spatiotemporal features; The third unit is used to calculate the local density value and relative distance value of the enhanced spatiotemporal features, perform density peak clustering to determine the center point of the functional area cluster, and generate the initial functional area division result; use the enhanced spatiotemporal features to construct an affinity matrix, combine the spatial adjacency constraints to perform spectral clustering operations, and generate optimized functional area division results; through the simulated annealing boundary optimization algorithm, iteratively optimize the internal consistency of the region and the difference between regions to generate the final urban functional area division results.

8. An electronic device, characterized in that: include: processor; a memory for storing processor-executable instructions; The processor is configured to call the instructions stored in the memory to execute the method according to any one of claims 1 to 6.

9. A computer-readable storage medium having computer program instructions stored thereon, characterized in that: When the computer program instructions are executed by a processor, the method according to any one of claims 1 to 6 is implemented.

Images