Cross-modal data fusion-based urban functional area identification method

Abstract

The invention discloses a cross-modal data fusion-based urban functional area identification method, and relates to the technical field of urban space-time big data analysis and geographic information processing, and the method comprises the steps: 1, constructing an analysis unit set according to a region boundary; 2, multi-source city data are obtained with the analysis units in the analysis unit set as units, the multi-source city data comprise built-up environment data and human activity data, and the multi-source city data are preprocessed; 3, constructing a population activity time image PATI according to the preprocessed human activity data, and performing dynamic time sequence feature vector extraction on the population activity time image PATI; and step 4, calculating a static built-up environment feature vector according to the preprocessed built-up environment data, fusing the static built-up environment feature vector with the dynamic time sequence feature vector to obtain a comprehensive feature vector, and clustering the comprehensive feature vector to obtain a function partition. The method gives consideration to robust recognition, mixing degree quantification and online updating capability, and is suitable for scenes of urban physical examination, planning evaluation, operation scheduling and the like.

Description

Technical Field

[0001] This invention relates to the field of urban spatiotemporal big data analysis and geographic information processing technology, and in particular to a method for identifying urban functional areas through cross-modal data fusion. Background Technology

[0002] As highly complex systems, cities' internal functional zones are often determined by a combination of activities and the physical environment. While traditional urban planning methods define functional zones based on initial zoning designs, actual functions often exhibit hybrid and dynamically evolving characteristics as cities expand and their socio-economic conditions change rapidly. Accurate classification and identification of urban functions are crucial for urban planning, resource allocation, and smart city development. Traditional research relying on remote sensing, GIS, and statistical yearbooks can depict macroscopic forms but struggles to reflect dynamic activity patterns such as intraday / weekly rhythms and seasonal migrations. With the increased availability of data such as mobile positioning, travel trajectories, and street view images, cross-modal data offers new opportunities for identifying hybrid functions and revealing evolutionary processes.

[0003] However, existing methods still face three key challenges: (1) single-modal or offline one-time segmentation is difficult to simultaneously characterize "static morphology + dynamic activity" and lacks the ability to quickly update to changes; (2) multi-source data differ significantly in sampling frequency, spatial granularity and quality, making cross-modal alignment and robust fusion difficult; (3) most outputs are hard classifications, unable to quantify mixed land use and uncertainty expression. Existing urban functional area identification methods have technical problems such as difficulty in simultaneously characterizing static physical space and dynamic human activities, difficulty in aligning and fusing multi-source heterogeneous data, rigid classification results and lack of uncertainty quantification, and inability to update incrementally with new data. Therefore, there is an urgent need for a technical approach that uses a unified spatial unit as a carrier, can integrate multiple modalities, output soft classifications and support continuous updates. Summary of the Invention

[0004] The technical problem to be solved by the present invention is to overcome the shortcomings of the prior art and provide a method for identifying urban functional areas by cross-modal data fusion. It integrates multi-source urban data on a unified spatial unit, outputs the membership degree of primary and secondary functions and the mixed entropy, and forms a serviceable vector result and a visual interpretation view, which is applicable to scenarios such as urban physical examination, planning evaluation and operation scheduling.

[0005] To solve the above-mentioned technical problems, the present invention adopts the following technical solution:

[0006] A method for identifying urban functional zones through cross-modal data fusion proposed according to the present invention includes:

[0007] Step 1: Construct a set of analysis units based on the region boundaries;

[0008] Step 2: Obtain multi-source city data by taking the analysis units in the analysis unit set as the unit. The multi-source city data includes built environment data and human activity data. Preprocess the multi-source city data.

[0009] Step 3: Construct a population activity time image (PATI) based on the preprocessed human activity data, and extract its dynamic temporal feature vector.

[0010] Step 4: Calculate the static built environment feature vector based on the preprocessed built environment data, and fuse it with the dynamic time series feature vector to obtain the comprehensive feature vector. Then, cluster the comprehensive feature vector to obtain functional partitions.

[0011] As a further optimization of the urban functional area identification method based on cross-modal data fusion described in this invention, step 1 involves constructing an analysis unit set, as detailed below:

[0012] The regional boundary is the city's administrative boundary or planning management boundary;

[0013] Based on the input region boundaries, with the region's top-left corner (x min , y max ) and the bottom right corner (x) max , y min Using Δx and Δy as the boundaries, generate the mesh by setting the mesh resolutions Δx and Δy.

[0014] ;

[0015] We obtain the unit set G={g i,j |i=1..N x j=1..N y};

[0016] Where, x min x max y min y max These are the left, right, lower, and upper boundaries of the rectangle enclosing the region boundary; Δx and Δy are the grid resolutions in the east-west and north-south directions, respectively; N x N is the number of grid columns in the horizontal x-direction. y It is the number of grid rows in the vertical y-direction; g i,j It is the grid in the i-th row and j-th column.

[0017] As a further optimization of the urban functional area identification method based on cross-modal data fusion described in this invention, in step 2, the built environment data includes buildings, road networks, green spaces, water systems, street view images, points of interest (POIs), and areas of interest (AOIs), while the human activity data is hourly population distribution data.

[0018] As a further optimization of the urban functional area identification method based on cross-modal data fusion described in this invention, step 2 includes preprocessed built environment data comprising preprocessed building and road network data, preprocessed green space and water system data, preprocessed street view image data, and preprocessed POI and AOI data; specifically as follows:

[0019] Step 2.1: Obtain Tianditu tiles from the Tianditu platform according to the regional boundaries. Based on the obtained Tianditu tiles, extract buildings and road networks according to color and vectorize them to form building vector data and road network vector data. The building vector data and road network vector data are the building and road network parts in the preprocessed built environment data.

[0020] Based on building vector data, calculate and analyze the building density, floor area ratio, and average height within the unit.

[0021] Based on the road network vector data, calculate the road network density, intersection density, and compactness centrality within the analysis unit;

[0022] Step 2.2: Based on the obtained Tianditu tiles, extract green areas and water systems according to color and vectorize them to form green area vector data and water system vector data; the green area vector data and water system vector data are the green area and water system parts in the preprocessed built environment data;

[0023] Calculate the green coverage rate within the analysis unit based on green space vector data;

[0024] Calculate the river coverage rate within the analysis unit based on the water system vector data;

[0025] Step 2.3: Based on the road network vector data obtained in Step 2.1, generate sampling points at preset intervals, and then call the Baidu Map API to obtain street view images based on the latitude and longitude coordinates of the sampling points, thus obtaining the street view image portion of the preprocessed built environment data.

[0026] The semantic segmentation network DeepLabV3+ was used to obtain the proportion of nine types of elements in each street view image, including roads, sidewalks, buildings, walls, plants, sky, people, cyclists, and cars. Based on this, four indicators corresponding to each sampling point were calculated: road aspect ratio, green view rate, sky visibility, and enclosure degree.

[0027] The four indicators of the analysis unit are calculated, and the average value of the sampling points falling into the analysis unit is used to represent the road height-to-width ratio, green view rate, sky visibility, and enclosure degree of the analysis unit.

[0028] Step 2.4: Call the Baidu Maps API to obtain POI and AOI data within the area boundary;

[0029] Based on the acquired POI data, a mapping is established from Baidu Maps POI industry classification to city function category set C, forming reclassified POI data; the reclassified POI data is the POI data portion in the preprocessed built environment data.

[0030] For each analysis unit g and each functional category c∈C, count the number of points of interest (POIs) n for each analysis unit g. g,c Calculate the POI density (POID) of the functional category set C in the analysis unit. g,c ;

[0031] The TF-IDF algorithm is used to calculate the characteristics of each city functional category in the analysis unit. :

[0032] ;

[0033] ;

[0034] ;

[0035] in, This represents the characteristics of class c functions within the analysis unit g; This indicates the relative proportion of class c functions within analysis unit g; n g,c c' represents the POI count of class c functions within analysis unit g; c' represents traversing any category among all function categories. IDF represents the total number of POIs across all functional categories within analysis unit g. c |G| represents the global rarity of class c functionality; |G| represents the total number of analysis units containing at least one POI within the region boundary; g' represents traversing all analysis units. This indicates the number of units containing class c functionality across all analysis units;

[0036] The preprocessed human activity data is as follows:

[0037] The Baidu Maps API is called to obtain population distribution data within the regional boundary and aggregated into various analysis units to form hourly population density data; the hourly population density data is the preprocessed human activity data.

[0038] As a further optimization of the urban functional area identification method based on cross-modal data fusion described in this invention, step 3 is as follows:

[0039] Based on preprocessed one-week hourly human activity data, a population activity time image PATI with a size of 7 days × 24 hours × 2 channels is constructed for each analysis unit g; where channel 1 represents the absolute population intensity channel and channel 2 represents the relative population change channel.

[0040] The PATI vector is input into a convolutional neural network (CNN) to extract the temporal pattern features of each analysis unit, thus obtaining the dynamic temporal feature vector z corresponding to analysis unit g. g .

[0041] As a further optimization of the urban functional area identification method based on cross-modal data fusion described in this invention, step 3 includes:

[0042] Step 3.1: Construct the PATI temporal image of population activity;

[0043] Based on the obtained population density data, the population density at time t on day d within each analysis unit g is calculated. d=1,…,7, where d is the weekday index, 1=Monday,…,7=Sunday; t=0,…,23, where t is the hour index; ≥0, according to P g (d,t) Construct a population activity time map with dimensions of 7 days × 24 hours × 2 channels:

[0044] ;

[0045] in, Represents the PATI tensor. Represents the real number field. This indicates channel 1 within analysis unit g at time t on day d. This represents the pixel value of channel 1 at time t on day d within analysis unit g. This represents channel 2 at time t on day d within analysis unit g. This represents the pixel value of channel 2 at time t on day d within analysis unit g;

[0046] The PATI temporal population activity image includes channel 1 and channel 2, where,

[0047] Channel 1 represents the absolute intensity channel, which characterizes the absolute population density intensity that is comparable across analysis units and time periods;

[0048] First, logarithmically compress the population density of analysis unit g, then normalize all analysis units within the region boundary to 0–1:

[0049] ;

[0050] in, This represents the population density value of analysis unit g at time t on day d after logarithmic compression. Represents the global The minimum value, Represents the global The maximum value, , This is a numerical stability term used to prevent the denominator from being zero; a larger pixel value in this channel indicates a higher absolute population density.

[0051] Channel 2 represents the relative change channel, which depicts the increase or decrease in population density of each analysis unit in adjacent time periods;

[0052] Log-smooth the population density of analysis unit g:

[0053] ;

[0054] Define the index prev(d,t) of the previous time period relative to the current time period (d,t), and connect them in a circular fashion hourly. That is, the time period before 0:00 is 23:00 of the previous day; and the time period before the beginning of the week is connected to the end of the week.

[0055] ;

[0056] in, This represents the index of the previous time period (d,t).

[0057] Calculate the adjacent differences between the analysis unit g and the previous time period:

[0058] ;

[0059] in, This represents the logarithmic difference between adjacent hours of the analysis unit g, i.e., the magnitude and direction of population density change;

[0060] To obtain the range [0,1] and be robust to extreme values, a robust normalization is performed on the one-week difference sequence of the analysis unit g, using 1% and 99% as quantiles:

[0061] ;

[0062] Among them, a g b represents the 1st percentile of the difference sequence of analysis unit g within one week. g This represents the 99th percentile of the difference sequence of analysis unit g over a week;

[0063] Step 3.2, for analysis unit g, As input, a convolutional neural network (CNN) is used to extract the pattern semantics from the PATI structure of the population activity time image, thus obtaining the dynamic temporal feature vector of that unit:

[0064] ;

[0065] in, This represents the dynamic temporal characteristic vector of the analysis unit g. This indicates that a Convolutional Neural Network (CNN) is used for convolutional feature extraction mapping. Represents dynamic time series feature vectors Dimensions Indicates length is The set of real number vectors.

[0066] As a further optimization of the urban functional area identification method based on cross-modal data fusion described in this invention, the method for obtaining the comprehensive feature vector in step 4 is as follows:

[0067] Based on the preprocessed images of buildings, road networks, green spaces, water systems, street views, and POI data, the static built environment feature vector s is calculated for analysis unit g. g ;

[0068] For each analysis unit g, the static physical environment feature vector s will be... g With dynamic temporal feature vector z g The combined feature vector y is obtained by fusing and dimensionality reduction. g .

[0069] As a further optimization of the urban functional area identification method based on cross-modal data fusion described in this invention, step 4 includes:

[0070] Step 4.1: Calculate the static built environment feature vector;

[0071] For analysis unit g, static physical environment features are constructed and standardized into vectors based on preprocessed images of buildings, road networks, green spaces, water systems, street views, and points of interest (POIs). ;

[0072] S g = [BD g FAR g H g , RD g ID g CC g GR g , RR g HWR g GVI g SVF g , ENG g ,POID g,c , w g,c ];

[0073] in, This represents the static built environment feature vector of the analysis unit g. Represents dynamic time series feature vectors Dimensions, BD gIndicates building density, FAR g Indicates floor area ratio, H g Indicates average height, RD g Indicates road network density, ID g Indicates intersection density, CC g Indicates compactness centrality, GR g Indicates green space coverage rate, RR g Indicates river coverage, HWR g Indicates the road's height-to-width ratio, GVI g Indicates green visibility ratio, SVF g Indicates sky visibility, ENG g Indicates enclosure degree, POID g Indicates POI density, w g,c Indicates the characteristics of urban functional categories;

[0074] Step 4.2: Integrate static built environment feature vectors and dynamic temporal feature vectors;

[0075] For the analysis unit g, its static built environment feature vector s g and dynamic temporal feature vector z g Weighted concatenation and PAC dimensionality reduction are performed to obtain the comprehensive feature vector. ;

[0076] ;

[0077] in, This represents the comprehensive eigenvector of the analysis unit g. Represents the comprehensive feature vector Dimensions The weights represent the feature vectors of the static built environment. The weights represent the dynamic temporal feature vectors. This indicates that PAC was used for principal component analysis dimensionality reduction.

[0078] Step 4.3: The comprehensive feature vector obtained in Step 4.2 { Input the K-Means clustering algorithm to divide the set of analysis units into clusters of equal number of categories, each cluster representing a type of urban functional area, thus completing the functional area identification.

[0079] As a further optimization of the urban functional area identification method based on cross-modal data fusion described in this invention, step 5 is included after step 4, and step 5 includes:

[0080] Step 5: Using the AOI data from the built environment data obtained in Step 2, construct a functional category prototype. Calculate the main function, secondary function, and functional mixing index of each analysis unit by comparing the comprehensive feature vector with the functional category prototype, and output the vector result diagram and interpretation view.

[0081] As a further optimization of the urban functional area identification method based on cross-modal data fusion described in this invention, step 5 is as follows:

[0082] Step S5.1: Based on the functional category set C, using the AOI data, select high-confidence regions. A high-confidence region refers to AOI data in which the functional category label exists and is unique. Use the high-confidence AOI as the prototype sample of the corresponding functional category c, and calculate the arithmetic mean of the comprehensive feature vectors of the analysis units covered by the prototype samples of each functional category to obtain the comprehensive feature mean vector of each functional category, which is used as the prototype of each functional category.

[0083] Step S5.2: Combine the mean vector of the comprehensive features corresponding to the functional category prototype with the comprehensive feature vector y of each analysis unit g. g Similarity measures are performed and functional membership is calculated to obtain the primary function, secondary function and membership of each analysis unit, and the functional mixture entropy and its normalized value are calculated.

[0084] Step S5.3: Output the global functional partitions and write the primary and secondary functions of each analysis unit and their membership degree into the feature attribute table.

[0085] Compared with the prior art, the present invention, employing the above technical solution, has the following technical effects:

[0086] (1) Improved spatial alignment and measurement consistency: Using regular grids as a unified analysis unit, heterogeneous data such as buildings, road networks, green spaces, water systems, streetscapes, POIs, and population density are spatially aligned at the same scale. This avoids the boundary offset and measurement error caused by repeated resampling and reprojection between multiple data sources in traditional methods, and allows indicators within the same unit to be directly compared and superimposed, thus improving the spatial consistency and repeatability of functional identification.

[0087] (2) Enhanced ability to depict static and dynamic elements: Simultaneously utilizing the temporal characteristics of static built environment and hourly population density within the same unit, it makes up for the shortcomings of previous methods based solely on land cover or POI that cannot reflect the intensity and rhythm of use during specific time periods. At the same time, in adjacent or overlapping areas of residential, commercial, and public service uses, it can effectively distinguish plots with similar forms but different time rhythms of use, thereby improving the accuracy of functional area boundary and type identification.

[0088] (3) Mixed functions can be quantitatively expressed to support refined management: output the membership degree and mixing degree of primary and secondary functions, replace the traditional hard classification results, intuitively quantify the degree and uncertainty of mixed land use, and carry out differentiated management based on the membership degree threshold and mixing degree level in applications such as urban renewal unit control planning, public service facility layout assessment, and use conflict investigation;

[0089] (4) Explainable and verifiable: Both static and dynamic features are traceable. The category prototype is obtained from the high-confidence area of ​​AOI. The category meaning is clear and it is easy to carry out consistency verification with the planned land use and facility layout.

[0090] (5) Universality: It relies on publicly available or widely accessible data sources, has a small number of parameters, a clear implementation path, and can be quickly deployed in different cities and scales. Attached Figure Description

[0091] Figure 1 This is a flowchart illustrating the method.

[0092] Figure 2 This is a schematic diagram of the Population Activity Time Image (PATI); where (a) is the absolute population intensity channel 1 and (b) is the relative population change channel 2.

[0093] Figure 3 The functional area identification results are based on Wuxi City as an example.

[0094] Figure 4 This is a schematic diagram of the random verification area of ​​the embodiment; where (a) is a residential area, (b) is a commercial service area, and (c) is a green space and plaza area. Detailed Implementation

[0095] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be described in detail below with reference to the accompanying drawings and specific embodiments.

[0096] The method of this invention uses a unified regular grid within the regional boundary as the analysis unit, and combines multi-source heterogeneous data such as buildings, road networks, green spaces, water systems, street view images, points of interest (POI), areas of interest (AOI), and hourly population distribution data. It forms a functional identification framework that takes into account both physical space and human activity attributes, can output functional labels and mixed land intensity, and supports rapid incremental updates.

[0097] Specifically, this invention takes the region boundary as input and automatically generates a regular grid as the analysis unit; it extracts multimodal features at the unit scale:

[0098] (1) Static built environment characteristics, including building form and density, road network structure indicators, street view semantics, green space and water system coverage, and POI facility structure;

[0099] (2) Dynamic temporal features: The population activity time image PATI dual-channel encoding is proposed to generate dynamic temporal images of population density (absolute population intensity + relative population change), and CNN is used for convolutional feature extraction and mapping.

[0100] The two types of features are concatenated and dimensionality reduced in the fusion layer to form a comprehensive feature vector. Then, a functional prototype is constructed based on the AOI and the similarity is calculated according to the comprehensive features. The main function, secondary function and membership degree of the analysis unit are output, and the functional mixed entropy and its normalized value are calculated. Functional partitions are obtained according to the main function.

[0101] The following effects can be achieved:

[0102] (1) Cross-modal unification: Using unit mesh as a common carrier, the problem of inconsistent granularity and coordinate inconsistency between different modes is solved;

[0103] (2) Temporal semantics can be learned: PATI transforms population density sequences into visual semantics, and CNN can capture weekly-daily rhythms and morphological differences;

[0104] (3) Soft classification and mixing degree: The main and secondary functions, membership degree and entropy value are calculated and output by comprehensive feature similarity calculation to explicitly quantify mixed land use and uncertainty;

[0105] (4) Explainable and portable: Static / dynamic features are visualized and traced; robust to missing data and noise, applicable to different city and regional scales.

[0106] To achieve the above objectives, the overall steps are as follows: Figure 1 As shown in the flowchart, the present invention specifically adopts the following technical solution:

[0107] 1. Region boundary input and unified analysis unit construction

[0108] Based on the input region boundaries, with the region's top-left corner (x min , y max ) and the bottom right corner (x) max , y min Using Δx and Δy as the boundaries, generate the mesh by setting the mesh resolutions Δx and Δy.

[0109] ;

[0110] We obtain the unit set G={g i,j |i=1..N x j=1..N y};

[0111] Where, x min x max y min y max These are the left, right, lower, and upper boundaries of the rectangle enclosing the region boundary; Δx and Δy are the grid resolutions in the east-west and north-south directions, respectively; N x N is the number of grid columns in the horizontal x-direction. y It is the number of grid rows in the vertical y-direction; g i,j It is the grid in the i-th row and j-th column.

[0112] Specifically, for object allocation, point features are assigned grids according to grid rules; polygon features are weighted and summarized based on area ratio; and line features are summarized based on length ratio or buffer area ratio. A unified projection to the metric coordinate system is used to ensure consistency in density and area calculations.

[0113] 2. Multi-source data acquisition and preprocessing

[0114] Built environment data includes buildings, road networks, green spaces, water systems, street view images, points of interest (POIs), and areas of interest (AOIs). Human activity data consists of hourly population distribution data.

[0115] 2.1 Obtain Tianditu tiles from the Tianditu platform according to the regional boundaries. Based on the obtained Tianditu tiles, extract buildings and road networks by color and vectorize them to form building vector data and road network vector data, thus obtaining the building and road network components in the preprocessed built environment data; calculate the building density, floor area ratio, and average height within the analysis unit based on the building vector data; calculate the road network density, intersection density, and compactness centrality within the analysis unit based on the road network vector data.

[0116] 2.2 Based on the acquired Tianditu (sky map) tiles, green areas and water systems are extracted by color and vectorized to form green area vector data and water system vector data, thus obtaining the green area and water system components in the preprocessed built environment data; based on the green area vector data, the green area coverage rate within the analysis unit is calculated; based on the water system vector data, the river coverage rate within the analysis unit is calculated.

[0117] 2.3 Based on the road network vector data, sampling points are generated at preset intervals. Then, the street view images are obtained by calling the Baidu Maps API based on the latitude and longitude coordinates of the sampling points, resulting in the street view image portion of the preprocessed built environment data. The semantic segmentation network DeepLabV3+ is used to obtain the proportion of nine types of elements in each street view image, including roads, sidewalks, buildings, walls, plants, sky, people, cyclists, and cars. Based on this, four micro-perception indicators, namely road aspect ratio, green view rate, sky visibility, and enclosure degree, are calculated for each sampling point. The micro-perception indicators of the analysis unit are calculated, and the average value of the sampling points falling into the analysis unit is used to represent the micro-perception indicators corresponding to the analysis unit.

[0118] 2.4 Use the Baidu Maps API to obtain POI and AOI data within the area boundary;

[0119] Based on the acquired POI data, a mapping is established from Baidu Maps POI industry classifications to a set of city function categories C={Residential, Commercial, Green, Public, LifeService, Finance, Company}, forming reclassified POI data and obtaining the POI data portion of the preprocessed built environment data. Among them, the city function categories are divided with reference to the major categories in the urban construction land classification system: Residential represents residential areas, Commercial represents commercial service areas, Green represents green space and square areas, Public represents public management and public utility areas, LifeService represents public service areas, Finance represents roads and transportation areas, and Company represents industrial and logistics warehousing areas.

[0120] For each analysis unit g and each functional category c∈C, count the number of points of interest (POIs) n for each analysis unit g. g,c Calculate the POI density (POID) of each functional category in the analysis unit. g,c ;

[0121] The TF-IDF algorithm is used to calculate the characteristics of each city functional category in the analysis unit. :

[0122] ;

[0123] ;

[0124] ;

[0125] in, This represents the characteristics of class c functions within the analysis unit g; This indicates the relative proportion of class c functions within analysis unit g; n g,c c' represents the POI count of class c functions within analysis unit g; c' represents traversing any category among all function categories. IDF represents the total number of POIs across all functional categories within analysis unit g. c |G| represents the global rarity of class c functionality; |G| represents the total number of analysis units containing at least one POI within the region boundary; g' represents traversing all analysis units. This indicates the number of units containing class c functionality across all analysis units;

[0126] Step 2.5: Call the Baidu Maps API to obtain population distribution data within the regional boundary, and aggregate it into each analysis unit to form hourly population density data, thus obtaining preprocessed human activity data.

[0127] 3. Construction of Temporal Population Activity Image (PATI) and Extraction of Dynamic Temporal Feature Vectors

[0128] 3.1 Constructing a temporal image of population activity using PATI

[0129] Based on the obtained population density data, the population density P at time t on day d within each analysis unit g is calculated. g (d,t), d=1,…,7, d is the weekday index, 1=Monday,…,7=Sunday; t=0,…,23, t is the hour index; P g (d,t)≥0, P g (d,t) represents the population density at time t on day d within the analysis unit g; according to P g (d,t) Construct a population activity time map with dimensions of 7 days × 24 hours × 2 channels:

[0130] ;

[0131] in, Represents the PATI tensor. Represents the real number field. This indicates channel 1 within analysis unit g at time t on day d. This represents the pixel value of channel 1 at time t on day d within analysis unit g. This represents channel 2 at time t on day d within analysis unit g. This represents the pixel value of channel 2 at time t on day d within analysis unit g;

[0132] like Figure 2 As shown, the PATI temporal population activity image includes channel 1 and channel 2, where,

[0133] like Figure 2In (a), channel 1 represents the absolute intensity channel, which characterizes the absolute population density intensity that is comparable across analysis units and time periods;

[0134] First, logarithmically compress the population density of analysis unit g, then normalize all analysis units within the region boundary to 0–1:

[0135] ;

[0136] in, This represents the population density value of analysis unit g at time t on day d after logarithmic compression. P' represents the global P' g The minimum value of (d,t) P' represents the global P' g The maximum value of (d,t) , This is a numerical stability term used to prevent the denominator from being zero; a larger pixel value in this channel indicates a higher absolute population density.

[0137] like Figure 2 In (b), channel 2 represents the relative change channel, which depicts the increase or decrease in population density of each analysis unit in adjacent time periods;

[0138] Log-smooth the population density of analysis unit g:

[0139] ;

[0140] Define the index prev(d,t) of the previous time period relative to the current time period (d,t), and connect them in a circular fashion hourly. That is, the time period before 0:00 is 23:00 of the previous day; and the time period before the beginning of the week is connected to the end of the week.

[0141] ;

[0142] in, This represents the index of the previous time period (d,t).

[0143] Calculate the adjacent differences between the analysis unit g and the previous time period:

[0144] ;

[0145] in, This represents the logarithmic difference between adjacent hours of the analysis unit g, i.e., the magnitude and direction of population density change;

[0146] To obtain the range [0,1] and be robust to extreme values, a robust normalization is performed on the one-week difference sequence of the analysis unit g, using 1% and 99% as quantiles:

[0147] ;

[0148] Among them, a g b represents the 1st percentile of the difference sequence of analysis unit g within one week. g This represents the 99th percentile of the difference sequence of analysis unit g over a week;

[0149] 3.2 Convolutional Feature Extraction using CNN

[0150] For analysis unit g, As input, a convolutional neural network (CNN) is used to extract the pattern semantics from the PATI structure of the population activity time image, thus obtaining the dynamic temporal feature vector of that unit:

[0151] ;

[0152] in, This represents the dynamic temporal characteristic vector of the analysis unit g. This indicates that a Convolutional Neural Network (CNN) is used for convolutional feature extraction mapping. Represents dynamic time series feature vectors Dimensions Indicates length is The set of real number vectors.

[0153] 4. Integration of static built environment characteristics and dynamic temporal characteristics

[0154] 4.1 Calculate the feature vector of the static built environment

[0155] For analysis unit g, static physical environment features are constructed and standardized into vectors based on preprocessed images of buildings, road networks, green spaces, water systems, street views, and points of interest (POIs). ;

[0156] S g = [BD g FAR g H g , RD g ID g CC g GR g , RR g HWR g GVI g SVF g , ENG g ,POID g,c , w g,c ];

[0157] in, This represents the static built environment feature vector of the analysis unit g. Represents dynamic time series feature vectors Dimensions, BD g Indicates building density, FAR g Indicates floor area ratio, H g Indicates average height, RD g Indicates road network density, ID g Indicates intersection density, CC g Indicates compactness centrality, GR g Indicates green space coverage rate, RR g Indicates river coverage, HWR g Indicates the road's height-to-width ratio, GVI g Indicates green visibility ratio, SVF g Indicates sky visibility, ENG g Indicates enclosure degree, POID g Indicates POI density, w g,c Indicates the characteristics of urban functional categories;

[0158] 4.2 Integrate static built environment feature vectors and dynamic temporal feature vectors;

[0159] For the analysis unit g, its static built environment feature vector s g and dynamic temporal feature vector z g Weighted concatenation and PAC dimensionality reduction are performed to obtain the comprehensive feature vector. ;

[0160] ;

[0161] in, This represents the comprehensive eigenvector of the analysis unit g. Represents the comprehensive feature vector Dimensions The weights represent the feature vectors of the static built environment. The weights represent the dynamic temporal feature vectors. This indicates that PAC was used for principal component analysis dimensionality reduction.

[0162] Step 4.3: The comprehensive feature vector obtained in Step 4.2 { The K-Means clustering algorithm is used to divide the set of analysis units into 7 clusters according to the number of categories in set C. Each cluster represents a type of urban functional area, thus completing the functional area identification. Figure 3 The image shows the functional area identification results using Wuxi City as an example.

[0163] 5. Comprehensive Feature Similarity Measurement and Functional Partition Output

[0164] 5.1 Screening of high-confidence samples for AOI

[0165] Based on the POI and AOI data obtained in step 2.4, for each functional category c∈C, a set of high-confidence areas A is selected based on the Baidu Map AOI classification labels. c , as a prototype of the functional category.

[0166] For all analysis units g∈G, the comprehensive features are used. Calculate the prototype mean and covariance for each functional category:

[0167]

[0168]

[0169]

[0170] Where c represents the function category index, A c This represents the set of high-confidence samples for functional category c. This represents the prototype mean vector of function category c. w represents the estimated weight of functional category c within analysis unit g. g,c This represents the characteristics of class c functions within the analysis unit g. This represents the diagonal covariance matrix of the functional category c, and diag(*) represents the operator that places the components of the input vector on the diagonal to form a diagonal matrix.

[0171] 5.2 Calculate the membership degree of functions

[0172] For all analysis units g∈G, the comprehensive eigenvectors are used. The Mahalanobis distance to the prototype mean vector of functional category c is used to obtain the similarity between the analysis unit g and each functional category, and then normalized to the functional membership degree.

[0173]

[0174]

[0175] in, The comprehensive eigenvector of the analysis unit g. Mahalanobis distance to the prototype mean vector of function category c τ represents the functional membership degree of analysis unit g to functional category c, and τ represents the numerical stability constant. This represents the matrix transpose operation.

[0176] Calculate the primary function, secondary function, function mix, and normalized function mix of the analysis unit g:

[0177]

[0178]

[0179]

[0180]

[0181] in, This indicates the main function of the analysis unit g. This indicates the secondary function of the analysis unit g. This indicates the functional mixing degree of the analysis unit g. The normalized value representing the functional mixing degree of analysis unit g.

[0182] 5.3 Presentation of Results

[0183] Write the primary and secondary functions, their membership properties, functional mix H, and normalized functional mix H_norm of each analysis unit into an appropriate feature attribute table. Output the membership degree and membership radar chart for each functional category corresponding to the analysis unit. (Example...) Figure 4 The diagram shown is a schematic of the random verification region output in the embodiment. Figure 4 In this context, (a)residential refers to a residential area. Figure 4 (b) commercial refers to the commercial services area. Figure 4 (c)green in the text refers to green spaces and plaza areas.

[0184] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for urban functional area identification through cross-modal data fusion, characterized in that, The application comprises the following steps: Step 1, constructing an analysis unit set according to a regional boundary; Step 2, acquiring multi-source urban data in the analysis unit set, the multi-source urban data comprising built environment data and human activity data, and preprocessing the multi-source urban data; Step 3, constructing a population activity time image (PATI) according to the preprocessed human activity data, and extracting a dynamic time sequence feature vector therefrom; Step 4, calculating a static built environment feature vector according to the preprocessed built environment data, fusing the dynamic time sequence feature vector and the static built environment feature vector to obtain a comprehensive feature vector, and clustering the comprehensive feature vector to obtain a functional district. 2.The urban functional area identification method of cross-modal data fusion according to claim 1, characterized in that, In step 1, the analysis unit set is constructed as follows: The regional boundary is a city administrative boundary or a planning management boundary; According to the input region boundary, with the region upper left corner (x min , y max ) and right bottom corner (x max , y min ) as boundary, set grid resolution Δx, Δy, generate grid: ； get the set of units G = {g i,j ∣i=1..N x ,  j=1..N y} wherein x min , x max , y min , y max are the left border, the right border, the lower border, and the upper border of the outer wrapping rectangle of the region border; Δx, Δy are the resolution of the grid in the east-west direction and the south-north direction; N x is the number of horizontal x-direction grid columns, N y is the number of vertical y-direction grid rows; g i,j is the grid in the i-th row and the j-th column. 3.The urban functional area identification method of cross-modal data fusion according to claim 1, characterized in that, In step 2, the built environment data comprises buildings, road networks, green spaces, water systems, street view images, points of interest (POIs), and areas of interest (AOIs), and the human activity data is hourly population distribution data.

4. The urban functional area identification method of cross-modal data fusion according to claim 3, characterized in that, In step 2, the preprocessed built environment data comprises preprocessed buildings, road network parts, preprocessed green space and water system parts, preprocessed street view image parts, and preprocessed POI and AOI parts, and is as follows: Step 2.1, acquiring a map slice according to the regional boundary from a map platform, extracting buildings and road networks according to the acquired map slice, and vectorizing the buildings and road networks to form building vector data and road network vector data, which are the buildings and road network parts of the preprocessed built environment data; According to the building vector data, the building density, volume rate, and average height in the analysis unit are calculated; According to the road network vector data, the road network density, intersection density, and closeness centrality in the analysis unit are calculated; Step 2.2, extracting green spaces and water systems according to the acquired map slice, and vectorizing the green spaces and water systems to form green space vector data and water system vector data, which are the green space and water system parts of the preprocessed built environment data; According to the green space vector data, the green space coverage rate in the analysis unit is calculated; According to the water system vector data, the river coverage rate in the analysis unit is calculated; Step 2.3, generating sampling points according to a preset interval based on the road network vector data obtained in step 2.1, and then calling a Baidu map API to obtain street view images according to the latitude and longitude coordinates of the sampling points, thereby obtaining the street view image part of the preprocessed built environment data; The semantic segmentation network DeepLabV3+ is used to obtain the proportions of a total of 9 elements, i.e., roads, sidewalks, buildings, walls, plants, skies, people, cyclists, and cars, in each street view image, and the road width-height ratio, green view rate, sky visibility, and enclosure degree of each sampling point are calculated based on the proportions; The four indexes of the analysis unit are calculated, and the average values of the sampling points falling in the analysis unit represent the road width-height ratio, green view rate, sky visibility, and enclosure degree of the analysis unit; Step 2.4, calling a Baidu map API to obtain POI and AOI data within the regional boundary; According to the obtained POI data, a mapping of the POI industry classification of Baidu Map to the city function category set C is established, and reclassified POI data is formed; the reclassified POI data is the POI data part in the preprocessed built environment data; For each analysis unit g and each functional category c e C, count the number of POIs n in each analysis unit g g,c ; calculate the POI density POID of the set of functional categories C in the analysis unit g g,c ; The TF-IDF algorithm is used to calculate and analyze the characteristics of each city function category in the unit : ； ； ； wherein, represents the characteristic of the c-th category of functions within the analysis unit g; represents the relative proportion of the c-th category of functions within the analysis unit g; n g,c represents the POI count of the c-th category of functions within the analysis unit g; c’ represents any of the categories of functions; represents the total number of POIs of all categories of functions within the analysis unit g; IDF c represents the global rarity of the c-th category of functions; |G| represents the total number of analysis units within the region boundary that contain at least 1 POI point; g’ represents all analysis units; represents the number of units containing the c-th category of functions among all analysis units; The preprocessed human activity data is as follows: The population distribution data within the region boundary is obtained by calling the Baidu Map API, and is aggregated to each analysis unit to form the hourly population density data; the hourly population density data is the preprocessed human activity data. 5.The urban functional area identification method of cross-modal data fusion according to claim 1, characterized in that, Step 3 is as follows: According to the preprocessed one-week hourly human activity data, a population activity time image PATI with a size of 7 days x 24 hours x 2 channels is constructed for each analysis unit g; wherein, channel 1 represents a population absolute intensity channel, and channel 2 represents a population relative change channel; The PATI is input into a convolutional neural network (CNN) to extract pattern features of each analysis unit in a time structure, and a dynamic time sequence feature vector z corresponding to the analysis unit g is obtained g .

6. The method of claim 3, wherein, Step 3 includes: Step 3.1, constructing a population activity time image PATI; According to the acquired population density data, the population density in each analysis unit g on the dth day at the tth hour , d = 1, …, 7, d is a week index, 1 = Monday, …, 7 = Sunday; t = 0, …, 23, t is an hour index; ≥ 0, according to P g (d, t) Construct a population activity time image with a size of 7 days x 24 hours x 2 channels: ； wherein, denotes the PATI tensor, denotes the real field, denotes the pixel value of channel 1 at day d and time t within analysis unit g, denotes the pixel value of channel 1 at day d and time t within analysis unit g, denotes the pixel value of channel 2 at day d and time t within analysis unit g, denotes the pixel value of channel 2 at day d and time t within analysis unit g, The population activity time image PATI includes channel 1 and channel 2, wherein, Channel 1 represents an absolute intensity channel, which depicts the absolute population density intensity across analysis units and time periods; First, the population density of the analysis unit g is logarithmically compressed, and then all analysis units within the region boundary are normalized to 0-1: ； wherein, represents the density value of the population density after logarithmic compression at the dth day and the tth hour of the analysis unit g, represents the global minimum value of represents the global maximum value of , is a numerical stability term for preventing the denominator from being 0; the greater the pixel value of this channel indicates the higher the absolute population density; Channel 2 represents a relative change channel, which depicts the increase and decrease amplitude of the population density of each analysis unit between adjacent time periods; The population density of the analysis unit g is logarithmically smoothed: ； Define the index prev(d,t) of the previous time period relative to the current time period (d,t), and perform ring connection hour by hour, i.e. the previous time period of 0 o'clock is 23 o'clock of the previous day; the first week is connected to the weekend again: ； wherein, denotes the index of the previous period of (d, t); Calculate the adjacent difference of the analysis unit g with the previous time period: ； wherein, represents the log-difference of adjacent hours of the analysis unit g, i.e. the magnitude and direction of the population density change; In order to obtain the range of [0,1] and be robust to extreme values, the one-week difference sequence of the analysis unit g is normalized with 1% and 99% as quantiles: ； where a g represents the 1% quantile of the difference series of the analysis unit g within a week, b g represents the 99% quantile of the difference series of the analysis unit g within a week; Step 3.2, for analysis unit g, the following is calculated: As input, the pattern semantics in the population activity time image, PATI, structure are extracted by a convolutional neural network, CNN, resulting in a dynamic time-series feature vector for this unit: ； wherein, denotes a dynamic timing feature vector of the analysis unit g, denotes a convolution feature extraction mapping with a convolutional neural network CNN, denotes a dynamic timing feature vector of dimension, denotes a set of real-valued vectors of length .

7. The method of claim 1, wherein, The method for obtaining the comprehensive feature vector in step 4 is as follows: According to the pre-processed building, road network, green land, water system, street view image, POI data, a static built environment feature vector s is calculated for the analysis unit g g ; For each analysis unit g a static physical environment feature vector s g with dynamic timing feature vector z g fusing and reducing dimensionality to a comprehensive feature vector y g . 8.The urban functional area identification method of cross-modal data fusion according to claim 1, characterized in that, Step 4 includes: Step 4.1, calculating a static built environment feature vector; For the analysis unit g, according to the pre-processed building, road network, green land, water system, street view image, point of interest POI data, the static physical environment features are constructed and standardized into vectors ; S g = [BD g , FAR g , H g , RD g , ID g , CC g , GR g , RR g , HWR g , GVI g , SVF g , ENG g , POID g,c ,w g,c ] wherein, represents a static built environment feature vector of the analysis unit g, represents a dynamic temporal feature vector dimension of BD g represents building density, FAR g represents volume rate, H g represents average height, RD g represents road network density, ID g represents intersection density, CC g represents closeness centrality, GR g represents green coverage, RR g represents river coverage, HWR g represents road height-width ratio, GVI g represents sky view factor, SVF g represents sky view factor, SVF g represents enclosure, POID g represents POI density, w g,c represents urban function category features; Step 4.2, fusing the static built environment feature vector and the dynamic time series feature vector; For the analysis unit g, the static built environment feature vector s g and the dynamic timing feature vector z g are weighted spliced and PAC dimension reduction is performed to obtain the comprehensive feature vector ; ； wherein, represents a comprehensive feature vector of the analysis unit g, represents a comprehensive feature vector dimension of the comprehensive feature vector, represents a weight of the static built environment feature vector, represents a weight of the dynamic timing feature vector, represents a principal component analysis dimension reduction using PAC; Step 4.3: The comprehensive feature vector obtained in Step 4.2 { Input the K-Means clustering algorithm to divide the set of analysis units into clusters of equal number of categories, each cluster representing a type of urban functional area, thus completing the functional area identification. 9.The urban functional area identification method of cross-modal data fusion according to claim 1, characterized in that, After step 4, step 5 is included, and step 5 includes: Step 5, using the AOI data in the built environment data obtained in step 2 to construct a function category prototype, calculating the main function, secondary function and function mixing degree index of each analysis unit by comparing the comprehensive feature vector with the function category prototype, and outputting a vector result map and an interpretation view.

10. The urban functional area identification method of cross-modal data fusion according to claim 9, characterized in that, Step 5 is as follows: Step S5.1, based on the function category set C, using the AOI data, selecting a high confidence area, which refers to the AOI data in which the function category label exists and is unique; taking the high confidence AOI as the prototype sample of the corresponding function category c, and performing arithmetic averaging on the analysis unit comprehensive feature vector covered by each function category prototype sample to calculate the comprehensive feature mean vector of each function category as the prototype of each function category; Step S5.2, the comprehensive feature mean vector corresponding to the functional category prototype is compared with the comprehensive feature vector y of each analysis unit g g The similarity measure is performed and the functional membership is calculated, the main function, the secondary function and the membership of each analysis unit are obtained, and the functional mixed entropy and its normalized value are calculated; Step S5.3, outputting a global function partition, and writing the main and secondary functions and their membership of each analysis unit into an attribute table.

Images