Urban functional area gridding identification and classification method based on multi-source spatio-temporal data

Abstract

The invention provides an urban functional area gridding identification and classification method based on multi-source spatio-temporal data, and relates to the technical field of area identification, and the method comprises the steps: carrying out the dynamic grid division of an urban area through employing a self-adaptive grid division algorithm, obtaining the multi-source spatio-temporal data of a gridding area, carrying out the space-time alignment, and generating a multi-dimensional time sequence feature set; extracting multi-dimensional feature vectors of the grid units through a multi-modal feature fusion algorithm, acquiring long-range time sequence features by using a bidirectional gating recursion unit, calculating spatial correlation features based on a multi-head attention mechanism, and inputting the two types of features into a feature fusion network with jump connection to generate enhanced spatial-temporal features; and an optimization strategy is adopted, a functional area clustering center point is determined through density peak value clustering in sequence to obtain an initial division result, the initial result is optimized by combining spectral clustering of spatial adjacency constraints, and the boundary of the functional area is optimized by using a simulated annealing algorithm, so that accurate identification and classification of the urban functional area are realized.

Description

Technical Field

[0001] The present invention relates to the technical field of region identification, and in particular to a method for gridding identification and classification of urban functional areas based on multi-source spatiotemporal data. Background Art

[0002] Identification of urban functional areas is an important basic work for urban planning and management. Traditional methods for identifying functional areas mainly rely on manual investigation and empirical judgment, which is not only time-consuming and labor-intensive, but also 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 identification method of urban functional areas based on spatiotemporal big data has gradually become a research hotspot. Existing data-driven methods usually adopt a fixed-size grid division method, combined with multi-source data such as POI data, human movement trajectory data, remote sensing image data, etc., to achieve automatic identification and classification of functional areas through feature extraction and machine learning algorithms.

[0003] However, the existing technology still has several shortcomings. The fixed-size grid division method is difficult to adapt to the spatial heterogeneity of urban area density, which easily causes the recognition deviation of functional area boundaries; the existing method simply splices the features of multi-source data, fails to fully explore the spatiotemporal correlation between data, resulting in insufficient feature expression ability; the clustering process of functional areas lacks consideration of the spatial adjacency relationship of grid units, which easily produces spatially fragmented recognition results; the optimization process of functional area boundaries fails to take into account both the internal consistency of the region and the differences between regions, affecting the accuracy of the final recognition results.

[0004] In summary, a method for gridding identification and classification of urban functional areas based on multi-source spatiotemporal data is proposed, which solves the problem of spatial heterogeneity through an adaptive grid partitioning algorithm, extracts temporal features using multimodal feature fusion and bidirectional gated recursive units, captures spatial correlation features using a multi-head attention mechanism, and combines density peak clustering, spectral clustering and simulated annealing boundary optimization algorithms to achieve accurate functional area identification. The present invention can solve the problems in the prior art. Summary of the invention

[0005] The embodiment of the present invention provides a grid-based identification and classification method for urban functional areas based on multi-source spatiotemporal data, which can solve the problems in the prior art.

[0006] According to a first aspect of the embodiments of the present invention, A method for identifying and classifying urban functional areas based on multi-source spatiotemporal data is provided, including: Adopting adaptive grid division algorithm to divide urban area into multiple grid units, dynamically adjusting the size of grid units based on regional density threshold, and generating gridded urban area; obtaining multi-source spatiotemporal data of gridded urban area, performing spatiotemporal alignment according to the preset number of time segments, and generating multi-dimensional time series feature set; performing feature extraction and fusion operation on multi-dimensional time series feature set through multi-modal feature fusion algorithm, and generating multi-dimensional feature vector of 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.

[0007] In an optional embodiment, 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.

[0008] In an optional embodiment, Based on the regional density threshold sequence, grid splitting and merging are performed to generate a gridded urban area, including: When the grid optimization function value is greater than the corresponding regional density threshold and the spatial dependency is less than a first preset threshold, the initial grid unit is non-uniformly split into four, and when the grid optimization function value of an adjacent grid unit is less than the corresponding regional density threshold and the spatial dependency is greater than a second preset threshold, the adjacent grid units are adaptively merged; The non-uniform quadripartitioning includes: calculating the density gradient field in the initial grid unit, determining the main direction and the secondary direction of the density change based on the density gradient field; calculating the relative importance measure of each position in the initial grid unit, and constructing a weight matrix; In the main direction, based on the cumulative distribution function of the weight matrix, the position of the first splitting line is determined so that the sum of the weighted density values ​​on both sides of the first splitting line is equal; in the secondary direction, the positions of the second splitting line are calculated for the two sub-regions divided by the first splitting line respectively so that the weighted density values ​​of the four sub-grids satisfy the balance constraint; Dividing the initial grid unit into four sub-grids, wherein 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; The adaptive merging includes: calculating the shape similarity and density difference of adjacent grid cells, and merging adjacent grid cells whose shape similarity is greater than a third preset threshold and whose density difference is less than a fourth preset threshold into a new grid cell, wherein the boundary of the new grid cell is obtained by fitting the boundary points of the adjacent grid cells using the least squares method.

[0009] In an optional embodiment, Obtain 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, 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.

[0010] In an optional embodiment, 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.

[0011] In an optional embodiment, 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.

[0012] In an optional embodiment, 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.

[0013] According to a second aspect of the embodiments of the present invention, A grid-based identification and classification system for urban functional areas based on multi-source spatiotemporal data is provided, including: 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 based on the multi-dimensional feature vector of the grid unit using a bidirectional gated recursive unit; calculate the multi-scale spatial correlation coefficient between the grid unit and the surrounding grids based on the multi-head attention mechanism and perform feature aggregation to generate spatial correlation features; input the long-range temporal features and spatial correlation features into a 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.

[0014] According to a third aspect of the embodiments of the present invention, An electronic device is provided, comprising: processor; a memory for storing processor-executable instructions; The processor is configured to call the instructions stored in the memory to execute the aforementioned method.

[0015] A fourth aspect of the embodiments of the present invention is: A computer-readable storage medium is provided, on which computer program instructions are stored. When the computer program instructions are executed by a processor, the aforementioned method is implemented.

[0016] In an embodiment of the present invention, through adaptive grid division and spatiotemporal alignment processing of multi-source spatiotemporal data, refined grid division and feature expression of urban areas are achieved, the accuracy and spatiotemporal consistency of data processing are improved, and a reliable data basis is provided for functional area identification; a combination of bidirectional gated recursive units and multi-head attention mechanisms is adopted to effectively capture the long-range temporal dependencies and multi-scale spatial correlation features of grid cells, and the expression ability of spatiotemporal features is enhanced through a feature fusion network, thereby improving the accuracy and robustness of functional area identification; combined with density peak clustering and spectral clustering methods, and the introduction of a simulated annealing boundary optimization algorithm, adaptive division and boundary optimization of functional areas are achieved, effectively balancing the internal consistency of the region and the differences between regions, and improving the scientificity and rationality of the functional area division results. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 It is a flow chart of a method for gridding identification and classification of urban functional areas based on multi-source spatiotemporal data according to an embodiment of the present invention; Figure 2 It is a line graph of mesh deformation rate under various density gradient conditions; Figure 3 Provide spatial distribution effect diagram for urban functional areas; Figure 4 It is a structural schematic diagram of the urban functional area grid identification and classification system based on multi-source spatiotemporal data according to an embodiment of the present invention. DETAILED DESCRIPTION

[0018] In order to make the purpose, technical solution and advantages of the embodiments of the present invention clearer, the technical solution in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0019] The technical solution of the present invention is described in detail with specific embodiments below. The following specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described in detail in some embodiments.

[0020] Figure 1 FIG. 1 is a flow chart of a method for identifying and classifying urban functional areas based on multi-source spatiotemporal data according to an embodiment of the present invention. Figure 1 As shown, the method includes: S101. Use an adaptive grid division algorithm to divide the urban area into multiple grid units, dynamically adjust the size of the grid unit 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 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; In this embodiment, an adaptive grid division algorithm is used to dynamically adjust the grid unit size according to the regional density, so that the urban area division is more in line with the actual spatial distribution, and the precision and rationality of the spatial representation are improved; the grid units are divided according to preset time segments, and the multi-source spatiotemporal data are aligned in time and space, which effectively captures the dynamic characteristics of the urban area that change over time and ensures the spatiotemporal consistency and integrity of the data; a multimodal feature fusion algorithm is used to extract and fuse multidimensional time series features, and the generated multidimensional feature vector can comprehensively reflect the spatiotemporal information in each grid unit, providing an accurate and reliable feature basis for the subsequent identification and classification of urban functional areas.

[0021] S102. Based on the multi-dimensional feature vector of the grid unit, a bidirectional gated recursive unit is used to extract long-range temporal features; based on the multi-head attention mechanism, the multi-scale spatial correlation coefficients of 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; In a specific implementation, after obtaining the multi-dimensional feature vector of the grid unit, a bidirectional gated recursive unit is first used to extract time series features. The bidirectional gated recursive unit includes two processing units, forward and backward, which extract features from the forward and reverse 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 degree of retention of historical information, and the reset gate is used to control the degree of reset of historical information. For the input multi-dimensional feature vector sequence of the grid unit, the time series features containing long-term dependencies are extracted through the cooperation of the update gate and the reset gate. Finally, the output features of the forward and backward processing units are fused to obtain the long-range time series features of the grid unit. Taking a certain grid unit as an example, its 24-hour feature sequence can effectively capture the activity patterns and their evolution characteristics of different time periods after being processed by the bidirectional gated recursive unit.

[0022] Then, the spatial correlation between grid cells is processed based on the multi-head attention mechanism. First, for each grid cell, the adjacent grid cells within a radius of 500 meters are determined 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 feature vector of the target grid cell is similar to the feature vector of the surrounding grid cells to obtain the attention weight score. After normalization, the attention weight score is weighted and combined with the feature vectors of the surrounding grid cells. The output features of the 8 attention heads are spliced ​​and linearly transformed to obtain a feature representation that reflects the spatial correlation of multiple scales. Based on this feature representation, the features of the target grid cell and its surrounding grid cells are weightedly aggregated to generate the spatial correlation features of the grid cell.

[0023] Finally, a feature fusion network with skip connections is constructed to achieve 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 high-level features can also guide the extraction process of the lower-level features. After multi-layer feature extraction and fusion, the enhanced spatiotemporal features with a dimension of 128 are finally output. This feature comprehensively expresses the temporal evolution law and spatial correlation characteristics of the grid unit. In the 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.

[0024] In this embodiment, bidirectional gated recursive units are used to process the multi-dimensional features of grid units, which can fully capture long-range temporal information and improve the ability to model dynamic changes in temporal series. Multi-scale spatial correlation coefficients are calculated based on a multi-head attention mechanism, and spatial correlation features are extracted through feature aggregation operations, which effectively reflects the complex spatial relationship between grid units and surrounding areas. Through a feature fusion network with jump connections, long-range temporal features and spatial correlation features are fused, which not only retains the original information, but also enhances the expression ability of spatiotemporal features, thereby providing more accurate and comprehensive feature support for subsequent identification of urban functional areas.

[0025] S103. 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.

[0026] In this embodiment, enhanced spatiotemporal features are used to calculate local density and relative distance, and density peak clustering is used to determine the center points of functional area clusters, thereby generating preliminary functional area division results, achieving fine capture of local features within the region; an affinity matrix is ​​constructed and combined with the spatial adjacency constraints of grid units, 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 results; the simulated annealing boundary optimization algorithm is used to iteratively optimize the internal consistency of the region and the difference between regions. Finally, the optimized urban functional area division results are obtained, achieving a balance between internal homogeneity within the region and significant differences between regions, providing scientific and accurate data support for urban planning and management.

[0027] 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: 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.

[0028] 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.

[0029] Then the target urban area is divided into initial grid cells, and the initial grid size is set to 500 meters × 500 meters. A multidimensional density assessment model is constructed, and 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, respectively, and the decoder is set symmetrically. Based on the extracted hidden layer features, the regional density distribution of the grid cells is calculated in combination with variational reasoning, and the density probability distribution is estimated by the Gaussian mixture model. The density distribution is input into the graph attention network, and the attention layer is set with 8 attention heads. The output dimension of each attention head is 64. The spatial dependency weights between adjacent grid cells are calculated by the softmax function to generate the grid optimization function value.

[0030] Then, an adaptive density threshold learning model is constructed based on the grid optimization function value. The grid optimization function value is modeled using a long short-term memory network. The network contains two LSTM layers with a hidden dimension of 128. The deep deterministic policy gradient algorithm is combined for reinforcement learning to dynamically adjust the density threshold sequence. The grid cells are divided into three levels: high density, medium density, and low density. The Adam optimizer is used for gradient descent, the learning rate is set to 0.001, and the grid size corresponding to different density levels is determined by 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 is kept at 500 meters × 500 meters.

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

[0032] Take the central area of ​​a city as an example, the 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 coarse grid of 500 meters × 500 meters. The total number of grids is about 400, which reduces the computational burden by about 40% compared with the unified grid division.

[0033] In this embodiment, adaptive grid division of urban areas is achieved by combining deep learning and reinforcement learning, which improves the accuracy of urban data analysis and management. The combination of deep autoencoders and graph attention networks fully explores the potential characteristics and spatial correlation of urban activity data, and improves the rationality of grid division; the adaptive density threshold learning model can dynamically adjust the grid size according to regional characteristics, which not only ensures the refined expression of high-density areas, but also avoids the waste of resources in low-density areas. The grid division method of the quadtree structure has good scalability and flexibility, and can be dynamically adjusted with regional development; it significantly reduces data storage and computing resource consumption, and improves urban management efficiency.

[0034] In an optional implementation, performing grid splitting and merging based on the region density threshold sequence to generate a gridded urban area includes: When the grid optimization function value is greater than the corresponding regional density threshold and the spatial dependency is less than a first preset threshold, the initial grid unit is non-uniformly split into four, and when the grid optimization function value of an adjacent grid unit is less than the corresponding regional density threshold and the spatial dependency is greater than a second preset threshold, the adjacent grid units are adaptively merged; The non-uniform quadripartitioning includes: calculating the density gradient field in the initial grid unit, determining the main direction and the secondary direction of the density change based on the density gradient field; calculating the relative importance measure of each position in the initial grid unit, and constructing a weight matrix; In the main direction, based on the cumulative distribution function of the weight matrix, the position of the first splitting line is determined so that the sum of the weighted density values ​​on both sides of the first splitting line is equal; in the secondary direction, the positions of the second splitting line are calculated for the two sub-regions divided by the first splitting line respectively so that the weighted density values ​​of the four sub-grids satisfy the balance constraint; Dividing the initial grid unit into four sub-grids, wherein 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; The adaptive merging includes: calculating the shape similarity and density difference of adjacent grid cells, and merging adjacent grid cells whose shape similarity is greater than a third preset threshold and whose density difference is less than a fourth preset threshold into a new grid cell, wherein the boundary of the new grid cell is obtained by fitting the boundary points of the adjacent grid cells using the least squares method.

[0035] In a specific implementation, first, the grid optimization function value is calculated based on the population distribution data in the initial grid cells. The population density of each grid cell is calculated by obtaining the census data and normalizing it. At the same time, the spatial dependency between adjacent grid cells is calculated, which is specifically obtained by calculating the correlation coefficient of the population density of adjacent grid cells.

[0036] When the grid optimization function value is greater than the corresponding regional density threshold (for example, set to 0.8) and the spatial dependency is less than the first preset threshold (for example, set to 0.3), the initial grid unit needs to be non-uniformly split into four. First, the density gradient field in the initial grid is calculated, and based on the grayscale image processing method, the first-order difference in the x-direction and y-direction is calculated to obtain the main direction and secondary direction with the largest density change.

[0037] Next, the relative importance of each location in the grid is calculated and a weight matrix is ​​constructed. The importance is measured by calculating the ratio of the mean and standard deviation of the population density in the surrounding area of ​​each location. For example, the weight of the center of a grid cell is 0.8, and the weight of the edge location is 0.2.

[0038] In the main direction, the first splitting line position 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 equal. For example, a grid cell divides the first splitting line at x=100, so that the weighted density values ​​on both sides are 0.5.

[0039] In the secondary direction, the positions of the second splitting lines are calculated for the two sub-regions divided by the first splitting line. 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.

[0040] Finally, the initial grid unit 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 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).

[0041] 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 dependency 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, which is obtained by comparing the shape features of the boundary curve.

[0042] At the same time, the density difference is calculated by calculating the difference coefficient of population density of adjacent grid cells. When the shape similarity is greater than a third preset threshold, such as 0.8, and the density difference is less than a fourth preset threshold, such as 0.3, the adjacent grid cells are merged into a new grid cell.

[0043] The boundary of the new grid cell is obtained by fitting the boundary points of the adjacent grid cells using the least squares method. First, the coordinates of the feature points on the boundary of the adjacent grid cells are extracted, and then the new boundary curve equation is obtained by least squares fitting.

[0044] like Figure 2 The grid deformation rate performance of different grid optimization methods under various density gradient conditions is shown. This technical solution effectively controls the degree of grid deformation through non-uniform four-splitting and adaptive merging strategies, especially in high-density gradient areas. The data shows that when the density gradient is 0, the grid deformation rate of this technical solution is only 2.1%, which is significantly better than 5.8% of the uniform four-splitting 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%, which is much lower than 37.5% of the uniform four-splitting method and 30.2% of the traditional adaptive grid method. Especially in the medium and high density gradient range of 20-25 kg / m³ / m, the gap between this technical solution and other methods is more significant, with grid deformation rates of 12.1%-16.4%, while the uniform four-splitting 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 and the method of determining the splitting line position based on the cumulative distribution function of the weight matrix, which makes the mesh division more in line with physical properties and avoids excessive deformation by limiting the area ratio of adjacent sub-grids (not exceeding the preset maximum ratio threshold), thereby improving the calculation accuracy while maintaining the mesh quality.

[0045] In this embodiment, adaptive splitting and merging of urban area grids can be achieved, ensuring the rationality and effectiveness of grid division. The heterogeneity of density distribution within the region is taken into account through the non-uniform four-splitting method, 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 change, ensuring the balance of sub-grid density distribution after splitting. The adaptive merging strategy is adopted to effectively avoid the problem of excessive grid subdivision; the least squares fitting method of the new grid boundary ensures the smoothness and continuity of the grid boundary, and improves the visualization and practicality of the grid division results; the overall solution has strong universality and practical value.

[0046] In an optional implementation, multi-source spatiotemporal data of a gridded urban area is obtained, spatiotemporal alignment is performed according to a preset number of time segments, and a multidimensional time series feature set is generated; and feature extraction and fusion operations are performed on the multidimensional time series feature set through a multimodal feature fusion algorithm to generate a multidimensional feature vector of a 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.

[0047] The POI data specifically refers to geographic data used to mark specific locations in maps or geographic information systems, 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 in the area, reflect the semantic characteristics of urban space, and assist in urban planning and management.

[0048] In a specific implementation, POI data, remote sensing image data and population trajectory data of a gridded urban area are obtained. The POI data includes attribute information such as category, name, and location, the remote sensing image data includes multispectral and panchromatic bands, and the population trajectory data includes timestamp and location information.

[0049] For the division of the time dimension of the grid unit, 24 hours are divided into 48 time segments at intervals of 30 minutes. In each time segment, three types of data are processed respectively: For POI data, the POI open status (open or closed) in each time segment is extracted to form a POI status sequence. For remote sensing image data, multispectral band combination is used to extract the type of ground objects (buildings, roads, vegetation, etc.) to form a ground object distribution sequence. For population trajectory data, the population in the grid is counted to form a population distribution sequence.

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

[0051] A 4-layer POI semantic tree is constructed, from the bottom to the top, which are specific POI, subcategory, main category and overall. The normalized information entropy of each layer of nodes is calculated, and semantic features are extracted from the top to the bottom through recursive morphological operations. 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.

[0052] A three-layer pyramid is constructed for the ground coverage data, with a bottom-level resolution of 2 meters, and then upsampled to 6 and 18 meters resolutions. From the bottom up, features are aggregated through maximum pooling, and from the top down, the boundary is optimized through conditional random fields to extract a 256-dimensional scene feature vector.

[0053] The time-varying kernel density of population density was calculated, the kernel function bandwidth was determined by cross-validation, and a 64-dimensional flow feature vector was extracted. The three types of feature vectors were concatenated and z-score normalized.

[0054] A feature fusion model based on autoencoder is constructed, and the feature weight is determined by minimizing the reconstruction error. The standardized features are weighted and summed to generate a 448-dimensional grid unit feature vector. Experiments show that the feature vector can effectively represent the spatiotemporal characteristics of the grid unit.

[0055] In this embodiment, through the fusion and alignment of multi-source heterogeneous data, a comprehensive characterization of urban elements is achieved, and the integrity and accuracy of feature expression are improved. The collaborative processing of multimodal data effectively reduces the limitations of a single data source; using methods such as multi-level semantic trees, multi-scale pyramids, and adaptive kernel density estimation, features are extracted from dimensions such as semantics, scenes, and population flow, enhancing the expressive power of features. Recursive morphological decomposition and hierarchical feature extraction ensure the multi-scale nature of features; based on the feature fusion method of minimizing reconstruction errors, feature weights are adaptively determined to avoid the subjectivity of manually setting parameters. Standardization processing and weighted combination improve the comparability and robustness of features. The overall solution has strong versatility and scalability.

[0056] In an optional implementation, based on the multi-dimensional feature vector of the grid unit, using a bidirectional gated recursive unit to extract long-range temporal features includes: 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.

[0057] In a specific implementation, multi-scale time series sampling is first performed on the multi-dimensional feature vector of the grid unit to construct feature sequences of different time granularities. It is divided according to the daily scale of 24 hours, the weekly scale of 7 days, and the monthly scale of 30 days. Taking a grid unit in a commercial district as an example, features such as passenger flow, consumption index, and commuting index are selected, and samples are taken once an hour on the daily scale to form a 24-dimensional sequence, once a day on the weekly scale to form a 7-dimensional sequence, and once a day on the monthly scale to form a 30-dimensional sequence.

[0058] The feature sequence of each time scale is processed by a bidirectional gated recursive unit that integrates a spiking neuron module. The spiking neuron module is a bionic simulation of the information processing mechanism of biological neurons, and contains three core components: membrane potential accumulation, threshold control, and synaptic plasticity. In the membrane potential accumulation unit, a dynamic membrane potential with an initial value of 0 is maintained. When a feature input is received, the feature value is multiplied by the synaptic weight and added to the membrane potential, while considering the exponential decay characteristics 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.

[0059] 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 from the last update is greater than the refractory period, the feature update is triggered and the membrane potential is reset to 0. The synaptic plasticity unit uses the timing-dependent plasticity rule to dynamically adjust the synaptic weight according to the time relationship of the feature updates of adjacent neurons. When the feature update of the latter neuron occurs in a short time after the previous neuron, such as when the time interval is 3 time steps, the synaptic weight is increased by 10%; when the time interval exceeds 10 time steps, the synaptic weight is reduced by 5%.

[0060] The bidirectional gated recurrent unit contains two processing branches: forward and backward. In the forward branch, the output of the update gate is transformed by the first spiking neuron module. When the accumulated feature value exceeds the threshold, the feature update is triggered and the timing code is generated. The output of the reset gate is processed by the second spiking neuron module with the same configuration. Based on the feature update sequence of the two modules, the forward hidden state sequence is generated. The backward branch uses the same structure to generate the backward hidden state sequence.

[0061] Taking the morning rush hour data of commercial district grid cells as an example, when the flow of people increases rapidly, the accumulation speed of membrane potential accelerates, triggering frequent feature updates; during the stable period of flow of people, the frequency of feature updates decreases. Through this mechanism, key moments of change in the data can be adaptively captured. For the generated hidden state sequence, the time interval between adjacent feature updates is calculated. For example, if the forward hidden state sequence is updated at time 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 are 0.4 and 0.6.

[0062] The forward and backward hidden state sequences are weighted and combined with the corresponding attention weights to obtain the fused hidden state sequence. Then, a convolution kernel of size 3 is used to perform temporal convolution on the fused sequence to extract the local temporal feature pattern. At the same time, the original feature vector is transformed nonlinearly to obtain the residual feature, which is added element by element to the local temporal feature to achieve feature enhancement. Finally, the feature value is normalized to zero mean and unit variance distribution through standardization to obtain the final long-range temporal feature.

[0063] The processed long-range time series features can effectively characterize the activity patterns of commercial areas at different time scales. The larger values ​​in the feature vector correspond to the typical activity periods of the area, such as the commuting mode during the morning and evening rush hours on weekdays, and the consumption peak on weekends. These features provide a reliable data basis for subsequent regional function identification 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.

[0064] In this embodiment, through the coupling design of pulse neurons and GRU, the timing characteristics of pulse neurons are fully utilized to improve the modeling ability of long-term temporal dependencies; a multi-scale layered sampling strategy is adopted to capture the law of feature changes at different time scales, thereby enhancing the model's perception of different periodic patterns; an attention mechanism is introduced to model the intensity of pulse emission, and combined with the residual connection design, the gradient vanishing problem in long sequence processing is effectively alleviated, thereby improving the robustness of feature extraction.

[0065] In an optional implementation, calculating the local density value and relative distance value of the enhanced spatiotemporal feature, performing density peak clustering to determine the functional area cluster center point, and generating the initial functional area division result includes: 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.

[0066] In a specific implementation, in order to realize the spatiotemporal feature analysis and functional area division of grid units, it is first necessary to calculate the comprehensive distance between grid units. Specifically, the distance relationship of grid units in physical space is calculated by Euclidean distance, the standardized distance of grid units in feature space is calculated by Mahalanobis distance, and the angle of grid unit feature vectors is calculated by cosine distance. The weight coefficients of different distance metrics are determined based on the variance contribution rate of the spatiotemporal features of each grid unit, and the three distances are weightedly fused to obtain a comprehensive distance matrix. For example, if the Euclidean distance of a grid unit 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.

[0067] Then, the comprehensive distance matrix is ​​sorted in ascending order, and the distance distribution curve is drawn using the kernel density estimation method. All local minimum points are obtained by calculating the first-order derivative of the curve, and the first distance that satisfies the condition that the derivative is negative and the amplitude is less than the mean of the curve is selected as the cutoff distance. For example, if the mean of a distance distribution curve is 0.6, and the corresponding distance of the first local minimum point that meets the conditions is 0.45, then 0.45 is used as the cutoff distance.

[0068] Then, 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 cutoff distance. The initial relative distance value is then calculated based on the local density value, and the minimum distance between the grid cell and the high-density value grid cell is calculated under multiple preset search radii (such as 0.1, 0.2, and 0.3). The weight is determined according to the stability of the relative distance value under different search radii, and the final relative distance value is obtained by weighted fusion.

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

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

[0071] In this embodiment, by integrating multiple distance measurement methods and comprehensively considering the similarities of grid units in physical space and feature space, the accuracy and robustness of distance calculation are improved; the adaptive Gaussian kernel function and multi-scale search strategy are used to calculate local density values ​​and relative distance values, which can better adapt to the data distribution characteristics of different regions and improve the accuracy of density peak clustering; the functional area clustering center points are selected through cluster center scoring and spatial uniformity constraints, and the grid unit allocation is carried out in combination with the density accessibility principle, which realizes a reasonable functional area division and ensures the spatial continuity and integrity of the division results.

[0072] In an optional implementation, an affinity matrix is ​​constructed using enhanced spatiotemporal features, and a spectral clustering operation is performed in combination with spatial adjacency constraints to generate an optimized functional area division result; a simulated annealing boundary optimization algorithm is used to iteratively optimize the regional internal consistency and inter-regional difference indicators to generate the final urban functional area division result, 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.

[0073] In a specific implementation, similarity is first calculated based on enhanced spatiotemporal features, and similarity is calculated for the feature vectors of each grid unit using a radial basis function. For example, the parameter is set to 0.5, and the Euclidean distance between grid units is substituted into the radial basis function to obtain a similarity matrix with similarity values ​​between 0 and 1. At the same time, the adjacency is calculated based on the relative position relationship of the grid units in space, with the adjacency strength of adjacent grid units being 1 and that of non-adjacent grid units being 0, to construct an adjacency constraint matrix.

[0074] The variance contribution rate of the similarity matrix and the adjacency constraint matrix to the grid unit classification is calculated respectively, and the variance contribution rate is used 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 the affinity matrix. Then, the sum of the elements in each row of the affinity matrix is ​​calculated to form the degree matrix, and the affinity matrix is ​​Laplace normalized using the degree matrix.

[0075] Perform eigendecomposition on the standardized affinity matrix to obtain a descending sequence of eigenvalues. By observing the trend of eigenvalue changes, cut off at the point where the eigenvalue drops significantly, and select the corresponding eigenvectors to construct the feature matrix. For example, select the first 5 eigenvectors. Normalize the row vectors of the feature matrix to obtain low-dimensional representation features. Perform K-means clustering in this feature space to redistribute the grid cells to different functional areas.

[0076] Calculate the variance of the grid cell characteristics within each functional area to obtain the consistency index within the region. Calculate the distance between the grid cell characteristics of adjacent functional areas to obtain the difference index between regions. The weighted combination of these two indicators constructs the optimization objective function. Set the initial temperature to 100 and the cooling coefficient to 0.95.

[0077] 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, the new solution is directly accepted; when the function value increases, the acceptance probability is calculated based on the current temperature. The temperature is gradually reduced until it reaches the stop temperature 1, and the final functional area division result is output.

[0078] like Figure 3 As shown in the figure, the spatial distribution results of the urban functional area division generated by this technical solution are shown, taking part of a central urban area as an example. The figure clearly presents five typical functional areas: commercial area, residential area, industrial area, public service area and green space leisure area, and each functional area is marked with a unique symbol grid unit. As can be seen from the figure, the functional areas generated by this technical solution have obvious spatial aggregation and boundary continuity. The internal consistency indexes of each functional area are: 0.094 for commercial area, 0.108 for residential area, 0.112 for industrial area, 0.103 for public service area and 0.087 for green space leisure area, all of which are lower than 0.12, indicating that the spatiotemporal characteristics of the grid units within each functional area are highly similar. At the same time, the difference index between adjacent functional areas is between 0.684 and 0.827, and the average reaches 0.741, indicating that the boundaries between different functional areas are clear and the characteristics are significantly different. From the perspective of spatial distribution, the commercial area is mainly located in the city center (near 116.30°E, 39.96°N), with an area of ​​25.4 square kilometers, adjacent to the residential area and public service area; the residential area is located in the southwest of the commercial area, with an area of ​​32.7 square kilometers, which is 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 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 leisure area is located in the south of the city (near 39.85°N), with an area of ​​18.9 square kilometers, which is the functional area with the highest internal consistency. This functional area distribution diagram conforms to the functional zoning principle in urban planning theory, and reflects the advantages of this technical solution in taking into account both spatial continuity and functional differences.

[0079] In this embodiment, the spatial continuity of the functional area division and the consistency of regional characteristics are guaranteed by the dual constraints of feature similarity and spatial adjacency constraints; the dimensionality reduction characteristics of spectral clustering are utilized to effectively extract the essential characteristics of the grid unit, reduce noise interference, and improve the accuracy and stability of the clustering effect; the simulated annealing algorithm is used to optimize the functional area boundaries, which ensures the internal consistency of the region while enhancing the differences between regions, making the functional area division results more reasonable and reliable.

[0080] In an optional embodiment, a high-tech industrial development zone in a certain city is selected as the research object, and the total area of ​​the research area is about 100 square kilometers. The method uses a regular grid of 250 meters by 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.

[0081] In the data preprocessing stage, the high-resolution remote sensing images are firstly processed by semantic segmentation. The deep convolutional neural network model is used to divide the study area into basic landform types such as buildings, water bodies, grasslands, woods, 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 marked as the corresponding natural land type. For the building area, the height of the building is estimated by calculating the length of the building's shadow projection, combining the shooting time and the solar altitude angle. The specific calculation process is: first extract the building outline and shadow boundary, calculate the shadow length L, and then calculate the building height H according to the formula H=L×tan(sun altitude angle). When the average height of the buildings in the grid is greater than 50 meters, it is preliminarily determined to be a business office area; when the building height is low and the area is large, it is preliminarily determined to be an industrial area.

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

[0083] In terms of crowd flow feature extraction, three consecutive weeks of mobile communication signaling data are used, divided into two situations: weekdays and weekends. For base station data, kernel density estimation is first 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. The crowd flow density distribution of a 250-meter grid is obtained by superposition calculation. For directly obtained grid data, time mean calculation is performed: the mean of 15 days of data is taken on weekdays, and the mean of 6 days of data is taken on weekends. Finally, the time series curve of the crowd flow of each grid for 24 hours on weekdays and weekends is obtained.

[0084] In order to identify the temporal characteristics of different functional areas, this method innovatively uses a self-organizing map neural network for cluster analysis. The network consists of two layers: an input layer and a competition layer. The number of nodes in the input layer is 48 (24-hour working days + 24-hour weekends), and the competition layer is set to a 6×6 two-dimensional node array. The network training process includes two steps: competitive learning and weight update. First, the Euclidean distance between the input vector and each competition layer node is calculated, and the node with the smallest distance is the winning node; then the weight vector of the winning node and its neighboring nodes is updated to make it closer to the input vector. Through repeated iterations, grids with similar temporal characteristics are eventually mapped to adjacent competition layer nodes. After cluster analysis, typical functional area temporal patterns are identified: the work area shows double peaks on weekdays and troughs on weekends; the commercial area shows continuous peaks during the day and significant increases on weekends; the residential area shows double peaks in the morning and evening, and is stable on weekends; the transportation hub area shows multiple peaks and weak regularity.

[0085] In terms of POI data analysis, this method uses the term frequency-inverse document frequency algorithm for processing. Each grid is regarded as a document, and the POI type in the grid is regarded as a term, and the TFIDF value of each type of POI is calculated. The calculation process is as follows: first, the frequency of occurrence of each type in the grid is counted, and the term frequency TF is obtained by dividing it by the total number of POIs in the grid; then the proportion of the number of grids containing this type to the total number of grids is calculated, and the inverse document frequency IDF is obtained by taking the logarithm of the reciprocal; the two are multiplied to obtain the TFIDF value, which is used to characterize the indicative effect of this type of POI on the grid function. Through this method, the dominant functional type of each grid can be identified, especially for the subdivision of the work area: when the administrative office type POI has the highest TFIDF value, it is determined to be an administrative office area; when the educational and scientific research type has the highest TFIDF value, it is determined to be an educational and scientific research area; when the medical and health type has the highest TFIDF value, it is determined to be a medical area.

[0086] Finally, this method uses a multi-label classification method to output the functional area recognition results. For each grid, the land use type, building characteristics, flow time series pattern and interest point distribution characteristics are comprehensively considered to determine its main function type and secondary function type. For example, a grid may have the attributes of "business office + commercial" at the same time, indicating that the area has both office buildings and commercial facilities; a grid may have the attributes of "residence + education", indicating that the area has both residential buildings and schools. This multi-label classification method is more in line with the complex characteristics of urban functions and can more accurately describe the actual use of the area.

[0087] Through application verification in the test area, the method shows good recognition effect. 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 is that it fully utilizes the complementarity of multi-source data, taking into account both spatial morphological characteristics and temporal activity characteristics, and can more comprehensively characterize the spatial distribution law of regional functions.

[0088] Figure 4 FIG. 1 is a schematic diagram of the structure of a grid-based urban functional area identification and classification system based on multi-source spatiotemporal data according to an embodiment of the present invention. Figure 4 As shown, the system comprises: 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 based on the multi-dimensional feature vector of the grid unit using a bidirectional gated recursive unit; calculate the multi-scale spatial correlation coefficient between the grid unit and the surrounding grids based on the multi-head attention mechanism and perform feature aggregation to generate spatial correlation features; input the long-range temporal features and spatial correlation features into a 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.

[0089] According to a third aspect of the embodiments of the present invention, An electronic device is provided, comprising: processor; a memory for storing processor-executable instructions; The processor is configured to call the instructions stored in the memory to execute the aforementioned method.

[0090] A fourth aspect of the embodiments of the present invention is: A computer-readable storage medium is provided, on which computer program instructions are stored. When the computer program instructions are executed by a processor, the aforementioned method is implemented.

[0091] 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 carrying computer-readable program instructions for executing various aspects of the present invention.

[0092] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or replace some or all of the technical features therein with equivalents. 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: Adopting adaptive grid division algorithm to divide urban area into multiple grid units, dynamically adjusting the size of grid units based on regional density threshold, and generating gridded urban area; obtaining multi-source spatiotemporal data of gridded urban area, performing spatiotemporal alignment according to the preset number of time segments, and generating multi-dimensional time series feature set; performing feature extraction and fusion operation on multi-dimensional time series feature set through multi-modal feature fusion algorithm, and generating multi-dimensional feature vector of 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 2, characterized in that Based on the regional density threshold sequence, grid splitting and merging are performed to generate a gridded urban area, including: When the grid optimization function value is greater than the corresponding regional density threshold and the spatial dependency is less than a first preset threshold, the initial grid unit is non-uniformly split into four, and when the grid optimization function value of an adjacent grid unit is less than the corresponding regional density threshold and the spatial dependency is greater than a second preset threshold, the adjacent grid units are adaptively merged; The non-uniform quadripartitioning includes: calculating the density gradient field in the initial grid unit, determining the main direction and the secondary direction of the density change based on the density gradient field; calculating the relative importance measure of each position in the initial grid unit, and constructing a weight matrix; In the main direction, based on the cumulative distribution function of the weight matrix, the position of the first splitting line is determined so that the sum of the weighted density values ​​on both sides of the first splitting line is equal; in the secondary direction, the positions of the second splitting line are calculated for the two sub-regions divided by the first splitting line respectively so that the weighted density values ​​of the four sub-grids satisfy the balance constraint; Dividing the initial grid unit into four sub-grids, wherein 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; The adaptive merging includes: calculating the shape similarity and density difference of adjacent grid cells, and merging adjacent grid cells whose shape similarity is greater than a third preset threshold and whose density difference is less than a fourth preset threshold into a new grid cell, wherein the boundary of the new grid cell is obtained by fitting the boundary points of the adjacent grid cells using the least squares method.

4. 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.

5. 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.

6. 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.

7. 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.

8. 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 7, 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 a 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.

9. 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 described in any one of claims 1 to 7.

10. 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 7 is implemented.

Images