CN113806419A - Urban area function identification model and method based on space-time big data - Google Patents

Abstract

The invention discloses an urban area function identification model and an urban area function identification method based on space-time big data, which are characterized in that: dividing urban plots to be researched according to an urban road network and a grid method to serve as target areas; constructing a dynamic time sequence diagram, and generating node characteristics and a graph adjacency matrix of the dynamic time sequence diagram; constructing a dynamic graph representation learning model GAT-GRU; calculating each land block B output by the model through the trained GAT-GRU model_iHigh level semantic feature representation of f_TAnd performing k-means clustering, and performing semantic recognition on a clustering result by combining POI distribution. The method is oriented to the development characteristics of large city density and high functional complexity in China, and realizes clustering and identification of functional areas with fine granularity on time by constructing a dynamic time sequence diagram, so that the method is more in line with the functional area identification thought of the development characteristics of cities in China; according to the invention, from the perspective of space-time interaction, a deep learning model is constructed to extract space-time characteristics, so that the automatic identification of the urban functional area is realized, and a new thought is provided for the traditional urban functional area research.

Description

Urban area function identification model and method based on space-time big data

Technical Field

The invention belongs to the technical field of data mining and urban calculation, and relates to an urban area function identification model and an urban area function identification method based on space-time big data.

Background

In recent years, along with popularization of relevant applications of smart cities, a large amount of real-time, high-precision and full-coverage mobile space-time big data generated by resident travel provides effective data support for relevant researches on urban area function identification. The analysis method driven by the space-time big data overcomes the defects of high economic cost, easiness in subjective judgment and the like in the traditional research based on questionnaire survey and remote sensing image data; in addition, the mobile space-time big data can describe the close relation between human activities and the urban spatial structure in a finer granularity, and provides powerful data support for effectively analyzing the human activities and researching the urban spatial structure. On the other hand, the successful application of the Deep Learning (Deep Learning) technology in the fields of image processing, auxiliary diagnosis and the like shows the capability of Deep Learning in modeling and solving complex problems, and provides a good technical basis for smart city analysis.

At present, much work is focused on clustering by combining with the movement space-time big data of residents, so that the functions of urban areas are effectively and timely identified. Generally, a resident travel mode is described Based on travel data and interest point data, and a machine learning algorithm is further used for identifying regional functions of a specific city, for example, Liu and Gao are Based on Urban Taxi GPS track data and interest point data (POI) data, and division and identification of Urban land functions are respectively realized by adopting Time Series data formed by combining a K-nearest neighbor algorithm (KNN) and a Gaussian Mixture Model (GMM) with track data (Liu X, Y tie, Zhang X, et al. The existing research shows that the effectiveness of city functional area identification by using the space-time data of resident movement is achieved, but in China, most cities have the development characteristics of high density and high land utilization mixing degree, and few single pure functional areas exist, which is also discouraged in the modern city planning concept of China. On the one hand, urban land areas present the characteristic of a combination of multiple land functions. On the other hand, there is a phenomenon that the function presented by the region of the land changes with time. However, the current urban area land use identification research still needs to be improved in functional concentration, and particularly, accurate identification which is more in accordance with the actual situation cannot be achieved for land areas with mixed functions (such as land areas where commercial and residential mixed buildings are located).

In addition, spatial correlation between urban plots is also important for region clustering and functional discovery. The urban spatial structure and the travel activities of residents mutually influence each other, and can be used as a research foundation for urban function discovery. Some existing studies considering spatial correlation of urban plots generally model a city into a graph structure, and analyze the graph structure based on a graph theory method. For example, Xu et al proposes a method for calculating affinity based on complex network node similarity to measure density distribution of various functions in an urban plot (yaying Xu, Gongfu L, Chao X, et al. affinity-based human mobility pattern for improved regional function discovery), shivers et al proposes a region pattern graph construction algorithm, builds a frequent trip pattern mining algorithm to discover frequently-occurring movement patterns in the region pattern graph, and finally realizes discovery of urban plot functions (e.g., wangyue, meiyi man, etc. the city function region discovery method based on trip pattern). Although the above existing method based on graph theory describes the relationship between the urban plots, it is not sufficiently effective to describe the time-varying characteristics of the interaction between the urban plots, that is, the relationship between the plots is not fixed and does not change in different time periods of the day, for example, the number of commutes between the working area and the residential area is higher than at other times during the early peak period, and therefore the corresponding relationship is more relevant than at other times.

Disclosure of Invention

Aiming at the problem of low urban plot function identification aggregation degree in the prior art, particularly the urban development phenomenon with high land utilization mixing degree, an urban area function identification model and an identification method based on space-time big data are provided, and the method is more suitable for functional area identification of urban development characteristics in China.

In order to achieve the purpose, the invention adopts the following technical scheme to solve the problem:

a city region function identification model and an identification method based on space-time big data specifically comprise the following steps:

step 1: dividing urban plots to be researched according to an urban road network and a grid method to serve as target areas;

step 2: constructing a dynamic time sequence diagram, and generating node characteristics and a graph adjacency matrix of the dynamic time sequence diagram;

the method comprises the following substeps:

step 21: each land block B in the target area B is taken as a node of the dynamic time sequence chart to generate each land block B_iAnd using the characteristic as a node characteristic of the dynamic time sequence chart;

step 22: calculating v time of each node in the target area B at integer time period t_i，v_zThe number of connected orders is recorded as a directed weighted graph G^(t)The adjacency matrix of (a);

defining a directed weighted graph G corresponding to an integer time period t^(t)Is A ═ w_iz)_N×NWherein w is_izRepresenting a node B_iAnd node B_zThe weights between, are calculated by:

wherein (B)_{j_up}，B_{j_down}) Indicating order j from starting point B_{j_up}To the end point B_{j_down}The OD (Origin-Destination) stream of (2), card represents the function of the number of elements of the solution, i.e. the starting point B_{j_up}To the end point B_{j_down}The number of connected orders is recorded;

and step 3: constructing a dynamic graph representation learning model GAT-GRU to obtain a trained GAT-GRU model;

and 4, step 4: calculating the four time periods of the working day peak in the morning and evening and the weekend peak in the morning and evening in a fixed time range respectively through the GAT-GRU model trained in the step 3, and outputting each land parcel B by the model_iIs represented by high-level semantic feature of f'_TAnd performing k-means clustering, and performing semantic recognition on a clustering result by combining POI distribution.

Further, the step 1 comprises the following sub-steps:

step 11, extracting road network data of a city to be researched, and determining an administrative region E according to the road network data; wherein E { (lon)₁，lat₁)，...，(lon_i，lat_i)，...，(lon_n，lat_n)}，(lon_i，lat_i) Representing longitude and latitude coordinates of a position point i in geographic space;

step 12, processing the administrative region E by applying a raster data model to obtain a land parcel set Z; wherein Z ═ { polygon ═₁，...polygon_j，...polygon_s}，polygon_jRepresents a rectangular ground plane composed of four points a, b, c and d in geographic space, and is represented as polygon_j＝{(lon_a，lat_a，lon_b，lat_b，lon_c，lat_c，lon_d，lat_d)}；

Step 13, taking an intersection of the administrative area E and the land parcel set Z to obtain a land parcel set B in the administrative area E, and using the land parcel set B as a target area B; wherein B ═ B₁，...，B_i，...，B_NN is the number of blocks in the administrative region E, B_iIs the ith parcel in the target area B.

Further, the step 21 includes the following sub-steps:

step 211, taking all taxi orders in the target area B within a fixed time range T to obtain an order set J;

wherein J ═ { J ═ J₁，j₂，...，j_i，...，j_n}. n is the number of all taxi orders in the fixed time range T of the target area B, j_iThe ith order in the order set J;

for a valid taxi order: j ═ lon_up，lat_up，lon_down，lat_down，t_up，t_downIn which lon_up，lat_upIndicates the boarding location, lon, of the order j_dowm，lat_downIndicating the drop-off location, t, of the order j_upFor getting-on corresponding time stamp, t_downA timestamp corresponding to the departure of the vehicle;

step 212, time-slicing each order using the formula, i.e., t for the order_up，t_downRespectively dividing the time periods into corresponding integer time periods, wherein the integer time periods refer to each hour in a fixed time range T;

t_up＝Floor(t_up)，t_down＝Floor(t_down)；

step 213, traverse the set of orders J, for each order, calculate the parcel B using the equation_iAt each integer time period t (t ∈ [0, 23 ]]) The number of getting-on and getting-off in the vehicle is as follows:

if (lon)_up，lat_up)∈B_iThen B is_iAt t_upNumber of getting on vehicles at time

Adding 1 in an accumulated way;

if (lon)_down，lat_dowm)∈B_iThen B is_iAt t_downNumber of alights at time

Adding 1 in an accumulated way;

step 214, traverse all plots B_iCalculating land parcel B_iNet input flow at integer time period t:

each land block B_iIs recorded as the ith node v in the dynamic time sequence chart G_iWill be

As node v_iThe characteristic of integer time period T is that the dynamic time sequence G ═ G in the fixed time range T under study⁽¹⁾，G⁽²⁾，...，G^(t)，...，G^(T)Is expressed as (f) corresponding to the node feature matrix₁，f₂，...，f_T) (ii) a Wherein G is^(t)Representing a directed weighted graph corresponding to an integer time period t.

Further, the step 3 comprises the following sub-steps:

step 31, designing a spatial layer representation module based on the weighted graph attention network through a GAT network model, and extracting spatial features of each node in the weighted graph;

the specific design process of the spatial layer representation module is to use a directed weighted graph G corresponding to an integer time period t^(t)In each node B_iAs a target node v_iCalculating a target node v by using the GAT network model shown in the formulas 2 and 3_iAttention weight between the target node and the neighbor node is obtained, namely the target node v_iRepresenting the feature vector after aggregating the features of the neighbor nodes;

in particular, w_ijRepresenting a target node v_iWith its neighbor node v_jAttention weight of between, wherein W^lIs a directed weighted graph G^(t)Target node v of l-th layer needing learning in GAT network_iShared weight transformation parameter of, N_iRepresenting a target node v_iSet of neighbor nodes of (A)_ij，A_ikRespectively represent target nodes v_iWith its neighbor node v_j，v_kWeight in between;

representing the feature vector representation of the target node vi at the l-th layer;

neighbor node v representing target node vi_j，v_kFeature vector representation at layer i; a is a parameterized vector to be learned by a network, | | represents series operation, and σ (·) is a nonlinear activation function;

for the target node v_iFeature vector representation obtained on the basis of the l-th layer

Next, the eigenvector representation of the L +1 layer network is calculated by formula 3

Step 32, designing a time layer representation module based on a recurrent neural network through a GRU network model;

for each node v_iExtracting T directed weighted graphs { G }⁽¹⁾，...，G^(t)，...，G^(T)The corresponding eigenvector representation in the GRU network model is used as the input of the GRU network model corresponding to T integer time periods, and then the corresponding directed weighted graph G of the GRU network model corresponding to the last integer time period T is output^(T)The feature vector representations of all the nodes at this time are combined to form a feature matrix f'_TAnd representing f 'as a high-level semantic feature of each node after spatial information and historical time information of neighbor nodes of each node are fused in the last integer time period T'_T；

Step 33, selecting a network model loss function, repeating the iteration steps 31-32, and performing optimization training on the dynamic graph representation learning model GAT-GRU to obtain a trained GAT-GRU model;

further, in step 33, a loss function is selected as shown in the following formula:

wherein the content of the first and second substances,

representing the probability of neighbor nodes around node v at time t,

representing the set of nodes that pass at a fixed random walk step of node v at time t,

represents G^(t)U denotes a certain node in the negative sample distribution, w_nIs an adjustable hyperparameter, represents a negative sampling rate, and sigma is a sigmoid activation function.

Further, the step 4 comprises the following sub-steps:

step 41, dividing the fixed time range T into four time periods of working day early peak, working day late peak, weekend early peak and weekend late peak; calculating the four time periods respectively through the dynamic graph representation learning model GAT-GRU obtained in the step 3, and outputting each land block B_iIs represented by high-level semantic feature of f'_T；

Step 42, training the four time periods to obtain a directional weighted graph G^(T)The feature vector of each node v represents f'_TAs input to the k-means algorithmThen get the directed weighted graph G^(T)Clustering results of all nodes in the cluster;

step 43, combining the clustering results of all the nodes obtained in step 42 with the node v_iCorresponding land parcel B_iAnd carrying out functional semantic annotation on the POI distributed on the POI.

Further, in the step 43, a WTF-IDF algorithm is adopted for functional semantic annotation.

Compared with the prior art, the invention has the following technical effects:

1) the method is oriented to the development characteristics of large city density and high functional complexity in China, and realizes clustering and identification of functional areas with fine granularity in time by constructing a dynamic time sequence diagram, so that the method is more in line with the functional area identification thought of the development characteristics of cities in China;

2) from the perspective of space-time interaction, a deep learning model is constructed to extract space-time characteristics, so that automatic identification of the urban functional area is realized, and a new idea is provided for traditional urban functional area research.

Drawings

FIG. 1 is a schematic diagram of a structure of a dynamic graph representation learning model according to the present invention.

FIG. 2 is a schematic diagram illustrating a step of dividing a target area in an embodiment; the three partial graphs from left to right respectively represent the original road network E, the grated region Z and the study region B.

FIG. 3 is a visualization of the results of the study in the examples; wherein, the four sub-graphs are respectively visualization graphs under four working modes.

FIG. 4 is a graph comparing the entropy of functional purity of the model of the present invention and related work (Liu, Gao) in cities.

FIG. 5 is a comparison verification chart of the actual geographic image data of the distribution of the Goodpasture map at the junior sports center according to the method of the present invention in the embodiment. Wherein, (a) is real geographical image data of the distribution of the high-grade map; (b) the flow trend graph of the original traffic flow data in four time modes is used.

The invention is further explained below with reference to the drawings and the detailed description.

Detailed Description

The invention provides an urban area function identification model and an urban area function identification method based on space-time big data, which specifically comprise the following steps:

step 1: dividing urban plots to be researched according to an urban road network and a grid method to serve as target areas;

step 11, extracting road network data of a city to be researched, and determining an administrative region E (as shown in fig. 2 (a)) according to the road network data; wherein E { (lon)₁，lat₁)，...，(lon_i，lat_i)，...，(lon_n，lat_n)}，(lon_i，lat_i) Representing the latitude and longitude coordinates of a location point i in geospatial space.

Step 12, processing the administrative region E by applying a raster data model to obtain a land parcel set Z (as shown in fig. 2 (b)); wherein Z ═ { polygon ═₁，...polygon_j，...polygon_s}，polygon_jRepresents a rectangular ground plane composed of four points a, b, c and d in geographic space, and is represented as polygon_j＝{(lon_a，lat_a，lon_b，lat_b，lon_c，lat_c，lon_d，lat_d)}；

Step 13, taking an intersection of the administrative area E and the land parcel set Z to obtain a land parcel set B in the administrative area E as a target area B (as shown in fig. 2 (c)); wherein B ═ B₁，...，B_i，...，B_NN is the number of blocks in the administrative region E, B_iIs the ith plot in the target area B;

step 2: constructing a dynamic time sequence diagram, and generating a node characteristic and a graph adjacency matrix of the dynamic time sequence diagram for fusing spatial and temporal relation characteristics between target areas B; the method specifically comprises the following substeps:

step 21: each land block B in the target area B is taken as a node of the dynamic time sequence chart to generate each land block B_iAnd as a node characteristic of the dynamic timing graph. The method specifically comprises the following steps:

step 211, taking all taxi orders in the target area B within a fixed time range T to obtain an order set J;

wherein J ═ { J ═ J₁，j₂，...，j_i，...，j_n}. n is the number of all taxi orders in the fixed time range T of the target area B, j_iThe ith order in the order set j;

specifically, for a valid taxi order: j ═ lon_up，lat_up，lon_dowm，lat_down，t_up，t_downIn which lon_up，lat_upIndicates the boarding location, lon, of the order j_down，lat_downIndicating the drop-off location, t, of the order j_upFor getting-on corresponding time stamp, t_downAnd the corresponding time stamp of the departure is shown.

Step 212, time-slicing each order using the formula, i.e., t for the order_up，t_downRespectively dividing the time periods into corresponding integer time periods, wherein the integer time periods refer to each hour in a fixed time range T;

t_up＝Floor(t_up)，t_down＝Floor(t_down)；

step 213, traverse the set of orders J, for each order, calculate the parcel B using the equation_iAt each integer time period t (t ∈ [0, 23 ]]) The number of getting-on and getting-off in the vehicle is as follows:

if (lon)_up，lat_up)∈B_iThen B is_iAt t_upNumber of getting on vehicles at time

Adding 1 in an accumulated way;

if (lon)_down，lat_down)∈B_iThen B is_iAt t_downNumber of alights at time

Adding 1 in an accumulated way;

step 214, traverse all plots B_iCalculating land parcel B_iAt an integer time period tNet input flow of (2):

each land block B_iIs recorded as the ith node v in the dynamic time sequence chart G_iWill be

As node v_iThe characteristics of the integral time period T, the dynamic time sequence diagram in the fixed time range T of the research

G＝{G⁽¹⁾，G⁽²⁾，...，G^(t)，...，G^(T)Is expressed as (f) corresponding to the node feature matrix₁，f₂，...，f_T). Wherein G is^(t)Representing a directed weighted graph corresponding to an integer time period t;

step 22: calculating v time of each node in the target area B at integer time period t_i，v_zThe number of connected orders is recorded as a directed weighted graph G^(t)The adjacency matrix of (a);

specifically, a directed weighted graph G corresponding to the integer time period t is defined^(t)Is A ═ w_iz)_N×NWherein w is_izRepresenting a node B_iAnd node B_zThe weights between, are calculated by:

wherein (B)_{j_up}，B_{j_down}) Indicating order j from starting point B_{j_up}To the end point B_{j_down}The OD (Origin-Destination) stream of (2), card represents the function of the number of elements of the solution, i.e. the starting point B_{j_up}To the end point B_{j_down}The number of orders linked between them.

Dynamic timing diagrams, also known as dynamic networks, where the structure of the graph evolves over time, as dynamic additions, subtractions, or changes in attributes of nodes and edges occur at a certain timestamp. Considering that the positions of taxis in each block of a city for getting on and off the taxi and the traffic are also dynamically changed, if the crowds in a certain block area have different destination demands under different timestamps, the correspondingly reflected traffic can also show different changing trends. Based on the time sequence diagram principle, the invention models the time and space characteristics of taxies in each block of the city.

And step 3: dynamic graph representation is designed to represent learning model GAT-GRU for further discovering each node B_i，B_zThe potential characteristics of the spatiotemporal relationship of the distributions; the method specifically comprises the following substeps:

step 31, designing a spatial layer representation module based on the weighted graph attention network through a GAT network model (namely the weighted graph attention network) to extract spatial features of each node in the weighted graph;

the specific design process of the spatial layer representation module is to use a directed weighted graph G corresponding to an integer time period t^(t)In each node B_iAs a target node v_iCalculating a target node v by using the GAT network model shown in the formulas 2 and 3_iAttention weight between the target node and the neighbor node is obtained, namely the target node v_iAnd representing the feature vector after aggregating the features of the neighbor nodes.

In particular, w_ijRepresenting a target node v_iWith its neighbor node v_jAttention weight of between, wherein W^lIs a directed weighted graph G^(t)Target node v of l-th layer needing learning in GAT network_iNi represents the target node v_iSet of neighbor nodes of (A)_ij，A_ikRespectively represent target nodes v_iWith its neighbor node v_j，v_kWeight in between;

representing the feature vector representation of the target node vi at the l-th layer;

representing a target node v_iV of a neighbor node_j，v_kFeature vector representation at layer i; a is a parameterized vector to be learned by the network, | | represents concatenation operation, σ (·) is a nonlinear activation function, and a leak relu function is adopted in the embodiment. Therefore, the target node v can be obtained by equation 2_iWith its neighbor node v_jAttention weight in between.

Further, for the target node v_iFeature vector representation obtained on the basis of the l-th layer

Next, the eigenvector representation of the L +1 layer network is calculated by formula 3

In step 31, contribution degrees of different neighbor nodes to the target node are automatically learned through the GAT network model, so that feature vector representation of the target node after neighbor node features are aggregated is obtained.

Step 32, designing a time layer representation module based on a recurrent neural network through a GRU network model;

after step 31, each integer time period t outputs a corresponding directed weighted graph G^(t)Representing the feature vector of each node v after fusing the features of surrounding neighbor nodes; then extracting T directed weighted graphs { G ] for each node vi⁽¹⁾，...，G^(t)，...，G^(T)The corresponding eigenvector representation in the GRU network model is used as the input of the corresponding T integer time periods in the GRU network model so as to output the corresponding T integer time periods in the last GRU network modelInteger time period T corresponds to a directed weighted graph G^(T)The feature vector representations of all the nodes at this time are combined to form a feature matrix f'_TAnd representing f 'as a high-level semantic feature of each node after spatial information and historical time information of neighbor nodes of each node are fused in the last integer time period T'_T。

Step 33, taking the formula 5 as a network model loss function, repeating the iteration steps 31-32, and performing optimization training on the dynamic graph representation learning model GAT-GRU to obtain a trained GAT-GRU model;

in the training process, the goal is for each node v_i(i.e., each plot B of a city_i) In the case that the node is expressed as a vector, the probability of the occurrence of the points adjacent to the node is expected to be the maximum, so that the original graph structure is reconstructed more accurately, and meanwhile, the embedded vector of each node can effectively represent high-level semantic information in the graph structure, so that the loss function shown in the formula 5 is selected.

Wherein the content of the first and second substances,

representing the probability of neighbor nodes around node v at time t,

representing the set of nodes that pass at a fixed random walk step of node v at time t,

represents G^(t)U denotes a certain node in the negative sample distribution, w_nIs an adjustable hyperparameter, represents a negative sampling rate, and sigma is a sigmoid activation function.

In step 3, aiming at the differentiated attributes of the urban functional areas in different time periods, the traditional static map method is difficult to effectively describe the dynamic evolution of the urban functional areas. The embodiment provides a dynamic Graph representation model based on a Graph Attention Network (GAT) and a Gated circulation Unit Network (GRU), which is used for representing the dynamic evolution rules of city functional areas in different time periods. Specifically, space structure information among all the land parcels at a certain moment is obtained through a GAT network, then historical flow change information of a certain land parcel is obtained through a GRU network, and space and time characteristic information of each land parcel of a city is extracted through a dynamic graph representation learning algorithm GAT-GRU.

And 4, step 4: for each land block B obtained in the step 3_iIs represented by high-level semantic feature of f'_TPerforming k-means clustering, and performing semantic recognition on a clustering result by combining POI distribution; the method specifically comprises the following substeps:

step 41, dividing the fixed time range T into four time periods of working day early peak, working day late peak, weekend early peak and weekend late peak; in the embodiment, the early peak is 7-9 points, and the late peak is 17-19 points; respectively representing the learning model GAT-GRU training by the dynamic graph obtained in the step 3 in the four time periods, and outputting each land parcel B_iIs represented by high-level semantic feature of f'_T；

Step 42, training the four time periods to obtain a directional weighted graph G^(T)The feature vector of each node v represents f'_TAs input to the k-means algorithm, a directed weighted graph G is then obtained^(T)Clustering results of all nodes in the cluster;

since the functional areas in the urban plots lack real labels, the embodiment applies a typical unsupervised algorithm k-means to cluster the features extracted by the learning model based on the dynamic graph representation in step 3. The basic idea of the k-means clustering algorithm is to divide samples into different categories according to the similarity between the samples, and the similarity calculation method adopted in this embodiment is an euclidean distance method. In addition, the clustering effect is influenced by the k value of the clustering number, and the optimal k value is selected according to two indexes of cluster inertia index (SSE) and contour coefficient (SH).

Step 43, combining the clustering results of all the nodes obtained in step 42 with the node v_iCorresponding land parcel B_iAnd carrying out functional semantic annotation on the POI distributed on the POI.

Preferably, a WTF-IDF algorithm is adopted for functional semantic annotation.

Example (b):

the embodiment is based on that 7065936 taxi order data of a target analysis city (city, Sichuan province) from No. 11/month No. 1 to No. 30/month in 2016 are obtained from a data open platform, and each record field comprises an order ID, order starting and ending time, and starting point longitude and ending point longitude and latitude corresponding to the order ID. And matching the getting-on and getting-off positions in the order record to the corresponding city plot according to the longitude and latitude to construct a net input stream (order ending quantity-order starting quantity) as the regional traffic data of the city plot at a certain time point.

The urban plot POI data of the embodiment is derived from a high-grade developer platform, about 550000 POI record points in 13 industry categories in a 2019 research area are crawled, and the classification is shown in table 1.

TABLE 1 POI data Classification

The embodiment respectively performs step 41 on the corresponding directed weighted graphs G in four time modes (working day morning peak T1, working day evening peak T2, weekend morning peak T3 and weekend evening peak T4)^(T)The feature vector of each node v in (a) represents f'_TK-means clustering is carried out, the two indexes of SSE and SH are comprehensively considered, the optimal k value is selected to be 6 (respectively expressed as Cl, C2, C3, C4, C5 and C6), and the WTF-IDF values corresponding to the 13 POI types are shown in Table 2. The first three POI type values with highest proportion of each category are added with bold font in the table and are used as the valuesAnd marking a standard for the functional semantics of the parcel. It can be seen from the table that after the time and space characteristics of the land parcel are extracted by the method provided by the invention, the function difference of each category is obvious, namely the clustering effect is good. Meanwhile, for the same type of land, corresponding functionality under different time modes (T1, T2, T3 and T4) is also different, the development characteristic of high functional complexity of the selected object in the embodiment is reflected, and the phenomenon can be well excavated.

TABLE 2 WTF-IDF values for each type of area POI under different time patterns

In the embodiment, a metropolis is taken as a target city to perform identification and analysis of functional areas, and in consideration of differences of resident travel on weekdays and weekends and different travel demands of early and late peaks every day, in this embodiment, four time pattern analysis corresponding to 22 working days of early peaks (7-9 points), late peaks (17-19 points) and 8 holidays within one month of 11 months are selected, and a visualization result is shown in fig. 3. The urban function is correspondingly changed along with the flowing dynamics of residents, so that the urban area function is not static and unchangeable, if certain specific areas are subjected to static analysis, the obtained result is a composite function, the method provided by the invention considers that the specific function presented by each area in each time period is captured better in a short time, and further the area division function with finer granularity on the time dimension can be realized. To illustrate the effectiveness of the method, the entropy of functional purity of urban areas is introduced to analyze the degree of functional aggregation of the areas, and the calculation formula 7 of the entropy of functional purity of urban areas is as follows:

wherein, Ent (C)_i) Represents the ith type land parcel B_iEntropy of functional purity of (1), p_kIs shown in the i-th type land block B_iAnd j represents the total number of the POI types.

The function aggregation degree of each type of region can be reflected by calculating the entropy value, and the smaller the entropy value is, the higher the POI distribution purity of the corresponding type of region is, the stronger the function aggregation degree is. Fig. 4 shows the results of the method of the present invention and the results provided by the work (Liu, Gao) in the same city, compared with the functional purity entropy in the city region, where max represents the record with the highest purity entropy in all the category regions, and min represents the record with the lowest purity entropy in all the category regions. As can be seen from fig. 4, for the regions where the actual distribution is multifunctional mixture, the method of the present invention exhibits the lowest entropy value, i.e. the highest degree of regional aggregation. This is because the method fully considers the dynamic changes of the region functions under different time modes, so that the identified function aggregation is higher. For the areas with functions biased to be unified, the areas where the "Chunxi roads" are located are divided into a single category in the result of the Gao recognition, so the entropy value is the lowest, the areas are clustered by the flow characteristics in the Gao, the flow of the areas where the "Chunxi roads" are always in a high value, and the areas are finally divided into a category.

In addition, in order to test the effect of the method on the dynamic identification of the urban functional area, the identification result obtained by the research of the embodiment is compared and analyzed with POI data and an area flow trend graph which are truly distributed on a Goodpasture map, and the experimental result shows that most of functional partitions which evolve along with time can be effectively identified through the method. Fig. 5 is a distribution diagram of real data points of a high-grade map of an area where a sports center of an adult city is located and a recording trend diagram of orders of getting on and off the bus in the area in one month on working days and on rest days. From a high-rise map distribution perspective, a number of government offices (shown as dark circles) are located near the metropolitan sports center, while a famous business center building wealth center is also located in this area. In addition, large scientific and educational service facilities such as sports centers, art theaters and art museums are distributed in the area, and the open time of the facilities is combined, so that the area can be preliminarily shown as a typical function mixed area, and groups can be in different needs to move in the area under different time modes. In combination with the flow trend chart, the getting-off quantity of the area is always larger than the getting-on quantity before 17 o' clock, namely, the situation of group inflow always indicates that the residents arrive at the area in the daytime because of certain functional requirements of the area. The traffic trends of the working days and the rest days before 19 o' clock are obviously different, the working days fluctuate more frequently, a plurality of extreme values are generated, and the commuting time characteristics are more obvious, so that the traffic flow distribution system is more consistent with the attributes of business offices. The whole fluctuation range is small during the rest day, the inflow in the morning is low, the inflow is gradually increased to about 10 points, the whole day is stable, the travel purpose of entertainment and leisure of residents is better met, and the probability that the region shows the functions of the cultural entertainment area is higher and the functions of the corresponding business area are weakened in the time mode by combining the distribution of POI. Through the combined analysis of the POI distribution and the flow trend in one day, the obtained result is basically consistent with the function evolution situation of the GAT-GRU method which only considers the flow characteristics and the space interaction characteristics of the early peak and the late peak (the interval between the first two vertical lines and the interval between the second two vertical lines in the graph 5).

Based on the result comparison and verification analysis, the effectiveness of the model provided by the invention on fine-grained city functional area identification is shown. By effectively identifying the complex and time-varying characteristics of urban land functions, urban regional structure attributes are comprehensively understood by urban planning researchers, effective information is provided for realizing urban planning with finer granularity, support is provided for tasks such as commercial site selection, traffic flow prediction and the like, and finally, the urban space is utilized more reasonably and efficiently to serve.

Claims (7)

1. A city region function recognition model and recognition method based on space-time big data are characterized by comprising the following steps:

step 1: dividing urban plots to be researched according to an urban road network and a grid method to serve as target areas;

step 2: constructing a dynamic time sequence diagram, and generating node characteristics and a graph adjacency matrix of the dynamic time sequence diagram;

the method comprises the following substeps:

step 21: each land block B in the target area B is taken as a node of the dynamic time sequence chart to generate each land block B_iAnd using the characteristic as a node characteristic of the dynamic time sequence chart;

step 22: calculating v time of each node in the target area B at integer time period t_i，v_zThe number of connected orders is recorded as a directed weighted graph G^(t)The adjacency matrix of (a);

defining a directed weighted graph G corresponding to an integer time period t^(t)Is A ═ w_iz)_N×NWherein w is_izRepresenting a node B_iAnd node B_zThe weights between, are calculated by:

wherein (B)_{j_up}，B_{j_down}) Indicating that order j is from starting point Bj _ up to ending point B_{j_down}The OD (Origin-Destination) stream of (2), card represents the function of the number of elements of the solution, i.e. the starting point B_{j_up}To the end point B_{j_down}The number of connected orders is recorded;

and step 3: constructing a dynamic graph representation learning model GAT-GRU to obtain a trained GAT-GRU model;

and 4, step 4: calculating the four time periods of the working day peak in the morning and evening and the weekend peak in the morning and evening in a fixed time range respectively through the GAT-GRU model trained in the step 3, and outputting each land parcel B by the model_iHigh level semantic feature representation of f_TAnd performing k-means clustering, and performing semantic recognition on a clustering result by combining POI distribution.

2. The spatio-temporal big data-based urban area function recognition model and recognition method according to claim 1, wherein the step 1 comprises the following sub-steps:

step 11, extracting road network data of a city to be researched, and determining an administrative region E according to the road network data; wherein E { (lon)₁，lat₁)，...，(lon_i，lat_i)，...，(lon_n，lat_n)}，(lon_i，lat_i) Representing longitude and latitude coordinates of a position point i in geographic space;

step 12, processing the administrative region E by applying a raster data model to obtain a land parcel set Z; wherein Z ═ { polygon ═₁，...polygon_j，...polygon_s}，polygon_jRepresents a rectangular ground plane composed of four points a, b, c and d in geographic space, and is represented as polygon_j＝{(lon_a，lat_a，lon_b，lat_b，lon_c，lat_c，lon_d，lat_d)}；

Step 13, taking an intersection of the administrative area E and the land parcel set Z to obtain a land parcel set B in the administrative area E, and using the land parcel set B as a target area B; wherein B ═ B₁，...，B_i，...，B_NN is the number of blocks in the administrative region E, B_iIs the ith parcel in the target area B.

3. The spatio-temporal big data-based urban area function recognition model and recognition method according to claim 1, wherein the step 21 comprises the following sub-steps:

step 211, taking all taxi orders in the target area B within a fixed time range T to obtain an order set J;

wherein J ═ { J ═ J₁，j₂，...，j_i，…，j_n}. n is all taxis in the fixed time range T of the target area BNumber of orders, j_iThe ith order in the order set J;

for a valid taxi order: j ═ lon_up，lat_up，lon_down，lat_down，t_up，t_downIn which lon_up，lat_upIndicates the boarding location, lon, of the order j_down，lat_downIndicating the drop-off location, t, of the order j_upFor getting-on corresponding time stamp, t_downA timestamp corresponding to the departure of the vehicle;

step 212, time-slicing each order using the formula, i.e., t for the order_up，t_downRespectively dividing the time periods into corresponding integer time periods, wherein the integer time periods refer to each hour in a fixed time range T;

t_up＝Floor(t_up)，t_down＝Floor(t_down)；

step 213, traverse the set of orders J, for each order, calculate the parcel B using the equation_iAt each integer time period t (t ∈ [0, 23 ]]) The number of getting-on and getting-off vehicles in the vehicle;

if (lon)_up，lat_up)∈B_iThen B is_iAt t_upNumber of getting on vehicles at time

Adding 1 in an accumulated way;

if (lon)_down，lat_down)∈B_iThen B is_iAt t_downNumber of alights at time

Adding 1 in an accumulated way;

step 214, traverse all plots B_iCalculating land parcel B_iNet input flow at integer time period t:

each land block B_iIs recorded as the ith node v in the dynamic time sequence chart G_iWill be

As node v_iThe characteristic of integer time period T is that the dynamic time sequence G ═ G in the fixed time range T under study⁽¹⁾，G⁽²⁾，...，G^(t)，...，G^(T)Is expressed as (f) corresponding to the node feature matrix₁，f₂，…，f_T) (ii) a Wherein G is^(t)Representing a directed weighted graph corresponding to an integer time period t.

4. The spatio-temporal big data-based urban area function recognition model and recognition method according to claim 1, wherein the step 3 comprises the following sub-steps:

step 31, designing a spatial layer representation module based on the weighted graph attention network through a GAT network model, and extracting spatial features of each node in the weighted graph;

the spatial layer represents a module specific design process: the directed weighted graph G corresponding to the integer time period t^(t)In each node B_iAs a target node v_iCalculating a target node v by using the GAT network model shown in the formulas 2 and 3_iAttention weight between the target node and the neighbor node is obtained, namely the target node v_iRepresenting the feature vector after aggregating the features of the neighbor nodes;

in particular, w_ijRepresenting a target node v_iWith its neighbor node v_jAttention weight of between, wherein W¹Is a directed weighted graph G^(t)Target node v of l-th layer needing learning in GAT network_iShared weight transformation parameter of, N_iRepresenting a target node v_iSet of neighbor nodes of (A)_ij，A_ikRespectively represent target nodes v_iWith its neighbor node v_j，v_kWeight in between;

representing a target node v_iFeature vector representation at layer i;

respectively represent target nodes v_iV of a neighbor node_j，v_kFeature vector representation at layer 1; a is a parameterized vector to be learned by a network, | | represents series operation, and σ (·) is a nonlinear activation function;

for the target node v_iThe feature vector obtained on the basis of the l-th layer represents f_i ^lNext, the eigenvector representation of the L +1 layer network is calculated by formula 3

Step 32, designing a time layer representation module based on a recurrent neural network through a GRU network model;

for each node v_iExtracting T directed weighted graphs { G }⁽¹⁾，...，G^(t)，...，G^(T)The corresponding eigenvector representation in the GRU network model is used as the input of the GRU network model corresponding to T integer time periods, and then the corresponding directed weighted graph G of the GRU network model corresponding to the last integer time period T is output^(T)The feature vector representations of all the nodes at this time are combined to form a feature matrix f'_TAnd representing f 'as a high-level semantic feature of each node after spatial information and historical time information of neighbor nodes of each node are fused in the last integer time period T'_T；

And 33, selecting a network model loss function, repeating the iteration steps 31-32, and performing optimization training on the dynamic graph representation learning model GAT-GRU to obtain a trained GAT-GRU model.

5. The spatio-temporal big data-based urban area function recognition model and recognition method according to claim 4, wherein in the step 33, a loss function is selected as shown in the following formula:

wherein the content of the first and second substances,

representing the probability of neighbor nodes around node v at time t,

representing the set of nodes that pass at a fixed random walk step of node v at time t,

represents G^(t)U denotes a certain node in the negative sample distribution, w_nIs an adjustable hyperparameter, represents a negative sampling rate, and sigma is a sigmoid activation function.

6. The spatio-temporal big data-based urban area function recognition model and recognition method according to claim 4, wherein the step 4 comprises the following sub-steps:

step 41, dividing the fixed time range T into four time periods of working day early peak, working day late peak, weekend early peak and weekend late peak; calculating the four time periods respectively through the dynamic graph representation learning model GAT-GRU obtained in the step 3, and outputting each land block B_iHigh level semantic feature representation of f_T；

Step 42, training the four time periods to obtain a directional weighted graph G^(T)The feature vector representation f of each node v in_TAs input to the k-means algorithm, a directed weighted graph G is then obtained^(T)Clustering results of all nodes in the cluster;

step 43, combining the clustering results of all the nodes obtained in step 42 with the node v_iCorresponding land parcel B_iAnd carrying out functional semantic annotation on the POI distributed on the POI.

7. The spatio-temporal big data-based urban area function recognition model and recognition method according to claim 6, wherein in the step 43, a WTF-IDF algorithm is adopted for functional semantic annotation.

Images