Identification of Urban Functional Zones Based on POI Density and Marginalized Graph Autoencoder

Abstract

With rapid urbanization, urban functional zones have become important for rational government and resource allocation. Points of interest (POIs), as informative and open-access data, have been widely used in studies of urban functions. However, most existing studies have failed to address unevenly or sparsely distributed POIs. In addition, the spatial adjacency of analysis units has been ignored. Therefore, we propose a new method for identifying urban functional zones based on POI density and marginalized graph autoencoder (MGAE). First, kernel density analysis was utilized to obtain the POI density and spread the effects of POIs to the surroundings, which enhanced the data from unevenly or sparsely distributed POIs considering the barrier effects of main roads and rivers. Second, MGAE performed feature extraction in view of the spatial adjacency to integrate features from the POIs of the surrounding units. Finally, the k-means algorithm was used to cluster units into zones, and semantic recognition was applied to identify the function category of each zone. A case study of Changzhou indicates that this method achieved an overall accuracy of 90.33% with a kappa coefficient of 0.88, which constitutes considerable improvement over that of conventional methods and can improve the performance of urban function identification.

1. Introduction

Urban areas are important spaces for humans to live in and produce [1]. With accelerating urbanization, sites of social activities gradually converge to form functional zones with distinct characteristics, such as residential, industrial, and commercial areas. Urban functional zones reflect the complex spatial structures and social functions [2]. These zones are important for understanding the spatial distributions of urban economic and social activities [3]. Determining the changes in urban functional zones and land use [4] is not only beneficial for urban decision making and policy formulation [5], but is also crucial for the economy, society, and ecology of urban areas. Urban functional zones are widely utilized in many studies of the spatial distribution of urban carbon emissions [6,7,8,9], the quality of urban development [10,11,12,13,14], and the urban heat island effect [15,16,17,18,19], among others. Thus, understanding the distribution of urban functional zones is essential for optimizing urban spatial structures and future city planning [20].

Urban functional zones are commonly identified using multisource geographic data, such as remote sensing data [21,22,23,24,25,26,27], points of interest (POIs) [28,29,30], vehicle trajectory data [31,32,33,34,35], and mobile phone data [36,37,38,39,40,41]. Urban remote sensing data are mainly used to extract features in terms of their natural attributes (e.g., spectral, shape, and textural features) [42]. However, the formation of urban functions depends on the spatial agglomeration of living and production activities. Information obtained from urban remote sensing data alone cannot effectively reflect social information. In the last decade, the emergence of geospatial big data has led to novel methods for urban function research. Cai et al. [41] used GPS trajectory and check-in data to construct a travel pattern graph, to which the Node2vec method was applied to discern the features of urban areas for identifying the functional zones. Zhou et al. [34] used taxi trajectory data and proposed an hour-day-spectrum (HDS) approach to generate distribution waveforms of taxi pick-up and drop-off points to identify different social functions. Song et al. [39] used aggregated mobile phone data containing work and residence attributes to automatically identify urban functional zones based on a self-organizing map neural network model.

POIs, as a type of geospatial big data, describe the main locations and information of urban human activity, record the detailed geographic location and functional labels of urban places, and have advantages of easy access and processing. POIs are mainly obtained via the external interfaces of commercial map providers, such as Baidu and Amap. POIs are increasingly used to identify urban functions and monitor urban dynamics. In particular, applying POIs for urban functional identification requires the extraction of the social features of the analysis units. Two types of methods are commonly used:

• Statistics-based methods: These methods count the number of POIs located within each unit to calculate indicators. The units are often divided by road network data or grids and are large in size in order to ensure that they contain sufficient POIs for subsequent analysis. The indicators include the frequency density, category factor, and term frequency-inverse document frequency (TF-IDF) to quantify the socioeconomic characteristics. Chen et al. [30] calculated the POI density of each urban functional zone using the POI quadrat density method to identify single functional, mixed functional, and no-data areas within urban areas. Xie et al. [1] used the category factors of POIs as feature vectors and fused remote sensing and trajectory data to identify urban functional zones. Chen et al. [27] calculated the TF-IDF index and frequency density of POIs within units to quantify the social features of urban functional zones. However, POIs are unevenly distributed in urban areas, whereas some emerging urban districts contain sparse or no POIs. Therefore, statistical methods are insufficient for obtaining regional social attributes and judging zone types.
• Spatial-context-based methods: These methods combine the spatial contextual relationships of POIs and use machine learning methods to mine potential semantic features for functional identification. Commonly used models include Word2Vec [43] and Doc2Vec. Zhang et al. [44] proposed the GeoSemantic2vec algorithm for extracting urban functional zones from POI semantic and location information. Sun et al. [45] proposed the Block2vec model that uses the skip-gram framework to map the spatial correlations between POIs as well as mapping the study units into high-dimensional vectors to identify urban functions. Zhai et al. [28] combined POIs and a simplified Place2vec model to construct a POI-based spatial contextual relationship to detect urban functional zones at a neighborhood scale. However, in addition to the spatial adjacency between POIs, such a relationship also exists between the study units within urban areas, which is neglected in the feature extraction process.

Thus, existing methods using POIs to identify urban functional zones have two shortcomings. First, they fail to address unevenly or sparsely distributed POIs, rendering it impossible to make reasonable inferences about emerging urban areas in the urban functional identification process. Previous studies have often identified these zones as undefined zones, which cannot be used to guide urban spatial structure optimization. Second, previous functional identification studies have ignored the spatial adjacency between the analysis units. As the units are not isolated from each other, there are correlations between neighboring units (e.g., commercial and residential units are often adjacent, whereas industrial and residential units are well separated) [46]. Ignoring these spatial adjacency results in deficiency in the spatial context.

To address these problems, we developed an urban functional identification model based on POI density and marginalized graph autoencoder (MGAE). The study had dual objectives: (1) to mitigate the problem of unevenly or sparsely distributed POIs, the model obtains the POIs density using kernel density analysis that considers barriers (e.g., main roads and rivers) and spreads the effects of POIs to the surrounding areas, and (2) to integrate the spatial context of POIs and units to generate the representation vectors, the model introduces MGAE, a machine learning model that can apply graph autoencoder to learn feature representation, for feature learning. The spatial adjacencies of unit and POI features were integrated using MGAE to improve the accuracy of urban function identification. Our study is finally able to identify the type of urban functions in the study area. It can better characterize the distribution of urban areas and provide references for urban planning.

2. Study Area and Dataset

2.1. Study Area

The downtown region of Changzhou was selected as the study area (Figure 1), which is surrounded by three expressways (i.e., Jiangyi, Huwu, and Hulong). Changzhou is located in southern Jiangsu Province in eastern China, and is a hub city in the Yangtze River Delta region. The city has a high level of urban development, an urbanization rate of over 70%, and diverse urban functions; thus, it is highly suitable for urban function identification. The study area covers 391.78 km² and has a population of over 10 million.

2.2. Data Sources and Preprocessing

The data used in this study include POIs and land use data for Changzhou. The interface provided by Baidu maps (https://lbsyun.baidu.com/, accessed on 17 May 2022) was utilized to obtain the POIs in 2021. The land use data were obtained from the Natural Resources and Planning Bureau of Changzhou.

2.2.1. POIs

A total of 6,137,875 POIs in Changzhou were acquired using the application programming interface (API), and coordinate transformation was performed. Each POI contained attributes including name, latitude, longitude, and label. However, the raw POIs included duplicated pickups and some types exhibited low relevance to urban function. Therefore, a partial dataset was extracted to obtain 185,018 POIs, which were reclassified into 14 categories (

Table 1. POI reclassification standards and numbers.

Reclassified Label	Description	Number	Proportion
Finance	ATMs; banks; pawn shops; credit unions	1871	1.01%
Car service	Car detailing; car sales; car accessories; car inspection centers	5092	2.75%
Corporate and factory	Plants and mines, parks; agriculture, forestry, and horticulture; office buildings, companies	25,058	13.54%
Culture and media	Exhibition galleries; cultural palaces; radio and television	590	0.32%
Food	Chinese restaurants; snack fast food restaurants; cake and dessert shops; foreign restaurants; bars	22,313	12.06%
Governmental and public organizations	Government agencies	5032	2.72%
Hotel	Star hotels; express hotels; hostels; guesthouses	1549	0.84%
Public facility	Medical; living services	21,692	11.72%
Residence	Interior buildings; dormitories; residential areas	40,906	22.11%
Science and education	Training institutions; primary schools; kindergartens; secondary schools; universities	4875	2.63%
Shopping mall	Department stores; shopping centers; home appliances; digital; stores; shopping areas	41,452	22.40%
Sports and recreation	Sports fitness; entertainment	5026	2.72%
Tourism attraction	Tourist attractions; water systems; natural features	1465	0.79%
Transportation facility	Entrances/exits; subway stations; bus stops; bus lines; subway lines; transportation facilities	8097	4.38%

), including finance, car service, culture and media, food, governmental and public organizations, hotel, public facility, residence, science and education, shopping mall, sports and recreation, tourism attraction, and transportation facility. The fields of each POI included the name, latitude, longitude, and reclassified label.

2.2.2. Land Use Data

In our study, land use data were used to divide the study area, which constituted vector data produced by the Natural Resource Administration of the Chinese government that is updated annually and comprises patches. Land use data were produced based on remote sensing image interpretations and field verification. These data can reflect the use of land at the time of investigation. The land use data had six categories, including farmland, woodland, grassland, water, residential land, and unused land. The patches in land use data were selected as the analysis units, and were assigned a function, such as “commercial”, following the evaluation of each patch. After removing roads and water, 36,710 units were used in the experiment.

3. Methodology

In this study, urban function was regarded as a comprehensive result impacted by POIs within the surrounding area, as well as the local structural information of the study area. Therefore, a method based on POI density and MGAE was proposed to identify urban functional zones. The general procedure of this method is illustrated in Figure 2. POIs were first preprocessed, including data de-duplication and reclassification, after which the effects of the POIs were spread to the surrounding area via kernel density analysis. Spatial adjacency between units was converted into a graph to represent local spatial structural features. The MGAE was used to integrate the spatial structure and content of units. Specifically, it could combine content between central and neighboring units in the graph and generate representation vectors.

Finally, the representation vectors were clustered using the k-means algorithm to divide the entire area into zones with relatively homogeneous urban functions, and semantic recognition was performed to determine the urban functional zone categories. Kernel density analysis, feature extraction, and clustering were performed using Python 3.0 as the programming language.

3.1. Kernel Density of POIs Considering Barrier Effects

Owing to the application of electronic map and location-based services (LBS), the number of POIs in urban areas has increased substantially; however, uneven POI distribution still exists. To address these problems, the kernel density analysis method was applied to fit the spatial distributions of POIs and spread their effects from a point to their surrounding areas. Then, the densities of different types of POIs were used to generate the features of each unit for subsequent identification [47].

Because the main roads and rivers within a city will hinder the flow of people and economic activity on both sides, the spread of POIs should stop when encountering a road or river. In this study, expressways, primary roads, railroads, and rivers were set as barriers in the kernel density calculations. The roads and water of the study area were based on data from Map World in the National Platform of Common Geospatial Information Services (https://www.tianditu.gov.cn/, accessed on 13 April 2023). When a barrier was present between POIs and the grid to be calculated, these POIs were considered as belonging to other grids and were not incorporated in the accumulated kernel density calculation (Figure 3). The kernel function used in this study was the quartic kernel; the kernel density of grid $\mathrm{P}(\mathrm{x},\mathrm{y})$ in space was calculated as follows:

$$\mathrm{P}(\mathrm{x},\mathrm{y}) =\frac{1}{{\mathrm{n}\mathsf{\pi }\mathrm{R}}^2}\sum _{\mathrm{i}=1}^{\mathrm{n}}\frac{3}{\mathsf{\pi }}{\left(1-{\left(\frac{{\mathrm{D}}_{\mathrm{i}}(\mathrm{x},\mathrm{y})}{\mathrm{R}}\right)}^2\right)}^2$$

(1)

where $\mathrm{R}$ is the bandwidth, $\mathrm{n}$ is the total number of POIs within the bandwidth range $\mathrm{R}$ and that have no barriers between the points and grid, and ${\mathrm{D}}_{\mathrm{i}}(\mathrm{x},\mathrm{y})$ is the distance between the grid center $\mathrm{P}(\mathrm{x},\mathrm{y})$ and the i-th POI within the bandwidth range $\mathrm{R}$ and without barriers between the point and grid (${\mathrm{D}}_{\mathrm{i}}(\mathrm{x},\mathrm{y}) <\mathrm{R}$).

The choice of the bandwidth plays an important role in the results [48,49]. An excessively large bandwidth will produce results that are more biased toward showing global-scale spatial variations, while local details are lost. In contrast, a very small bandwidth is more suitable for expressing local spatial variations, but misses the overall spatial variations. After considering road and water barriers and spatial effects in the study area, the final bandwidth used in this study was taken as 100 m.

3.2. Extracting Features Using MGAE

Each unit is not isolated in an urban area, but is interconnected with and influenced by the surrounding units. Therefore, the urban function depends not only on the social attributes of a single unit, but also on the influences of the surrounding units. In this study, the units and their adjacent areas were converted into graphs to reflect the spatial structural features, and the MGAE was applied to perform unsupervised feature extraction of the local structure and generate the representation vectors [50].

A graph construction procedure is illustrated in Figure 4. The graph consists of nodes and edges; binary adjacency matrix A is often used to describe graph structure. The graph is generally defined as $\mathrm{G}=(\mathrm{V},\mathrm{E},\mathrm{A})$. The node set $\mathrm{V}=\{{\mathsf{\nu }}_0,{\mathsf{\nu }}_1,\dots ,{\mathsf{\nu }}_{\mathrm{n}-1}\}$ represents patches in the land use data, and each node ${\mathsf{\nu }}_{\mathrm{i}}$ contains a d-dimensional feature vector $\{{\mathrm{x}}_1,{\mathrm{x}}_2,\dots ,{\mathrm{x}}_{\mathrm{d}}\}$ that represents the density values of various types of POIs in patches. $\mathrm{E}=\{{\mathrm{e}}_0,{\mathrm{e}}_1,\dots ,{\mathrm{e}}_{\mathrm{m}-1}\}$ is the edge set. Each ${\mathrm{e}}_{\mathrm{i}}$ denotes an edge connecting two nodes. If a common edge exists between two patches, the two corresponding nodes are connected by one edge.

The MGAE model is a graph-based learning method with the advantage of not requiring labeled data for training. The MGAE model has the advantage of good performance and a fast training speed. It can also use the convolutional operation on the graph to extract structural features (Figure 4). The MGAE model replaces $\mathrm{f}\left(\mathrm{X}\right)$ in the autoencoder loss function (Equation (2)) with the convolution formula ${\stackrel{\mathrm{~}}{\mathrm{D}}}^{-\frac{1}{2}}\stackrel{\mathrm{~}}{\mathrm{A}}{\stackrel{\mathrm{~}}{\mathrm{D}}}^{-\frac{1}{2}}\stackrel{\mathrm{~}}{\mathrm{X}}\mathrm{W}$ in the graph convolutional network, and adds a regularization term to construct the loss function as follows:

$$\mathrm{L}={‖\mathrm{X}-{\stackrel{\mathrm{~}}{\mathrm{D}}}^{-\frac{1}{2}}\stackrel{\mathrm{~}}{\mathrm{A}}{\stackrel{\mathrm{~}}{\mathrm{D}}}^{-\frac{1}{2}}\stackrel{\mathrm{~}}{\mathrm{X}}\mathrm{W}‖}^2+\mathsf{\lambda }{‖\mathrm{W}‖}_{\mathrm{F}}^2$$

(2)

$$\stackrel{\mathrm{~}}{\mathrm{A}}=\mathrm{A}+\mathrm{I}$$

(3)

where $\mathrm{A}$ denotes the adjacency matrix of the graph data, $\stackrel{\mathrm{~}}{\mathrm{D}}$ denotes the degree matrix of $\stackrel{\mathrm{~}}{\mathrm{A}}$, $\mathrm{W}\in {\mathrm{R}}^{\mathrm{d}\times \mathrm{d}}$ is the parameter matrix, ${‖\mathrm{W}‖}_{\mathrm{F}}^2$ is the regularization term, $\mathsf{\lambda }$ is a trade-off, $\mathrm{X}$ is the original vector matrix, and $\stackrel{\mathrm{~}}{\mathrm{X}}$ is the vector matrix after random corruption.

3.3. Urban Function Clustering

Based on the vectors representation generated by the MGAE, further cluster analysis was conducted to divide the entire area into zones with relatively homogeneous urban functions. Cluster analysis is an unsupervised classification process that groups similar patches into the same cluster, that is, patches of similar urban function are identified in one cluster. In urban function studies, various clustering algorithms have been widely used, including the k-means algorithm and the OPTICS clustering algorithm. In this study, the k-means algorithm was used for cluster analysis.

The clustering results were evaluated using external and internal test methods. The external test method was used to compare the clustering results with the correct classification results of the dataset to determine the merits of the clustering method. The internal test method relied only on the dataset to evaluate the accuracy of the clustering results. As no correct labels existed for the dataset, the internal test method was used as the criterion for setting the training parameters of the model. We used two clustering internal test methods to select the optimal value of cluster number M: the sum of the squared error (SSE), as calculated using Equation (4), and the silhouette coefficient (SIL), as calculated using Equation (5) [51]:

$$\mathrm{SSE}=\sum _{\mathrm{i}=1}^{\mathrm{m}}\sum _{{\mathrm{x}}_{\mathrm{i}}\in {\mathsf{\Omega }}_{\mathrm{i}}}{\left|{\mathrm{x}}_{\mathrm{i}} - {\overline{\mathrm{x}}}_{\mathrm{i}}\right|}^2$$

(4)

where $\mathrm{m}$ is the number of functional zones, ${\mathsf{\Omega }}_{\mathrm{i}}$ denotes all units classified into the i-th functional zone, ${\mathrm{x}}_{\mathrm{i}}$ denotes the representation vector of the unit, and ${\overline{\mathrm{x}}}_{\mathrm{i}}$ denotes the mean value of all units in the i-th functional zone.

$$\mathrm{SIL}=\frac{1}{\mathrm{n}}\sum _{\mathrm{j}=1}^{\mathrm{n}}\frac{a-b}{\underset{}{\max }\left(a,b\right)}$$

(5)

where n is the number of units, a denotes the average distance of the j-th unit from the units of interzones, and b denotes the average distance of the j-th unit from the other units in intrazones.

3.4. Urban Functional Semantic Recognition

Although the geometry of the urban functional zones was extracted after clustering, the functional types of the zones remained unclear. Therefore, the semantic identification of each zone type was needed to determine the urban function of each zone. Thus, two common metrics for partitioning were used: frequency density (FD) and category factor (CF) [52,53]. The metrics were calculated using the following equations:

$${\mathrm{FD}}_{\mathrm{ki}}=\frac{{\mathrm{V}}_{\mathrm{ki}}}{\sum _{\mathrm{j}\in {\mathsf{\Omega }}_{\mathrm{i}}}{\mathrm{S}}_{\mathrm{j}}}$$

(6)

where ${\mathrm{V}}_{\mathrm{ki}}$ denotes the number of the k-th POI type in the i-th functional zone, ${\mathsf{\Omega }}_{\mathrm{i}}$ denotes the set of all units assigned to the i-th functional zone, and ${\mathrm{S}}_{\mathrm{j}}$ denotes the area of the j-th unit. A higher ${\mathrm{FD}}_{\mathrm{ki}}$ indicates that the k-th type of POI is more densely distributed in the i-th functional zone.

$${\mathrm{CF}}_{\mathrm{ki}}=\left({\mathrm{V}}_{\mathrm{ki}}/{\mathrm{V}}_{\mathrm{i}}\right)/\left({\mathrm{V}}_{\mathrm{k}}/\mathrm{V}\right)$$

(7)

$${\mathrm{V}}_{\mathrm{k}}=\sum _{\mathrm{i}=0}^{\mathrm{M}}{\mathrm{V}}_{\mathrm{ki}}$$

(8)

$${\mathrm{V}}_{\mathrm{i}}=\sum _{\mathrm{k}=0}^{\mathrm{T}}{\mathrm{V}}_{\mathrm{ki}}$$

(9)

$$\mathrm{V}=\sum _{\mathrm{k}=0}^{\mathrm{T}}\sum _{\mathrm{i}=0}^{\mathrm{M}}{\mathrm{V}}_{\mathrm{ki}}$$

(10)

where, ${\mathrm{V}}_{\mathrm{i}}$ denotes the total number of the various POIs in the i-th functional zone, and ${\mathrm{V}}_{\mathrm{k}}$ and $\mathrm{V}$ denote the numbers of the k-th type of POI and the total value of the POIs, respectively. ${\mathrm{CF}}_{\mathrm{ki}}$ denotes the level of aggregation of the k-th POI type in the i-th functional zone and ${\mathrm{CF}}_{\mathrm{ki}}$ compares the value of the k-th type of POIs in the i-th functional zone with the average value of all zones. A larger ${\mathrm{CF}}_{\mathrm{ki}}$ indicates that the k-th type of POIs is more concentrated in the i-th functional zone. In particular, if the ${\mathrm{CF}}_{\mathrm{ki}}$ value is greater than 1, it suggests that the average of the k-th POI type in the i-th functional zone is greater than the average of all k-th POI type in all zones. This means that the level of aggregation of the k-th POI type in the i-th functional zone is relatively high in all zones.

4. Results and Analysis

4.1. Determination of Model Parameter

To test the clustering performance and determine the optimal number of clusters (M), we varied the number of clusters from 3 to 19 and repeated the analysis 10 times with each set of parameters to ensure the accuracy of the results. Figure 5 shows that the SSE and SIL gradually decreased as the number of clusters increased. The slope of the SSE curve varied widely when K ranged from 5 to 7, indicating that the optimal number of clusters should be within that range based on the elbow rule. When M = 5, the SIL curve reached a local maximum. Therefore, we selected M = 5 as the number of clusters for the experiment and subsequent analysis.

4.2. Accuracy Evaluation and Performance Comparison

A confusion matrix was used to evaluate the accuracy of the urban functional zone identification. Sixty patches in each of the five types of urban functional zones were selected randomly. Then, high-resolution remote sensing images and Map World data were utilized to artificially judge the urban functions. The evaluation results are listed in

Table 2. Confusion matrix results of urban functional zones.

	Cluster 1	Cluster 2	Cluster 3	Cluster 4	Cluster 5	Total
Cluster 1	58	1	0	1	0	60
Cluster 2	0	53	0	1	6	60
Cluster 3	0	2	58	0	0	60
Cluster 4	0	6	0	51	3	60
Cluster 5	8	1	0	0	51	60
Total	66	63	58	53	60	300

. The overall accuracy (OA) of the urban functional zone obtained using the method proposed herein reached 90.33%, and the kappa coefficient was 0.88.

To evaluate the urban function recognition performance, the proposed model was compared with commonly used models, including the Latent Dirichlet Allocation (LDA), TF-IDF, and marginalized Stacked Denoising Autoencoder (mSDA) [54]. mSDA is an unsupervised learning method for attribute data mining that uses linear denoisers as the basis for learning feature representation, simplifying the parameter estimation method, and reducing the training time. In total, 300 manually identified units were selected for each comparison (

Table 3. Overall accuracies and kappa coefficients of the functional area classification models.

Index	Model	OA	Kappa
1	LDA	62.68%	0.52
2	TF-IDF	57.35%	0.54
3	mSDA	34.58%	0.26
4	Proposed model	90.33%	0.88

). Because TF-IDF considers only the distribution frequency of the POIs and ignores their semantic information, this method often has poor identification performance. Although LDA and mSDA extract semantic information of POIs, they ignore the spatial relationships between units in feature extraction; thus, they have lower accuracies than the proposed method. These experiments show that incorporating spatial relationships into attribute feature extraction can mine the features of districts more effectively, which in turn allows better performance in extracting urban functional zones.

4.3. Urban Functional Zone Results

In Changzhou, five functional zones were identified using the proposed method (Figure 6 and

Table 4. Overall POI frequency densities and category factor.

Frequency Density
	C1	C2	C3	C4	C5
Science and education	1.953	40.677	439.445	79.841	19.102
Corporate and factory	42.986	125.100	639.795	306.240	92.937
Food	5.615	206.407	2101.777	410.648	70.115
Hotel	0.249	14.094	120.966	35.756	4.673
Finance	0.219	18.934	118.130	53.171	3.252
Tourism attraction	1.554	7.767	38.747	10.348	8.696
Public facility	7.822	187.733	2076.261	374.303	74.047
Governmental and public organizations	2.880	43.085	60.483	66.548	24.395
Sports and recreation	1.330	50.711	357.227	96.667	14.791
Car service	3.328	30.691	27.406	111.475	25.363
Shopping mall	13.945	359.059	4127.951	776.116	132.032
Transportation facility	2.362	47.689	238.151	105.754	21.385
Residence	17.253	399.193	639.795	538.780	200.408
Culture and media	0.199	4.840	31.186	10.937	2.586
Category Factor
	C1	C2	C3	C4	C5
Science and education	0.715	0.986	1.485	0.998	1.025
Corporate and factory	3.092	0.596	0.425	0.752	0.980
Food	0.456	1.109	1.575	1.139	0.834
Hotel	0.290	1.085	1.298	1.420	0.796
Finance	0.212	1.210	1.052	1.753	0.460
Tourism attraction	1.976	0.654	0.455	0.449	1.620
Public facility	0.651	1.035	1.596	1.065	0.904
Governmental and public organizations	1.037	1.027	0.201	0.819	1.288
Sports and recreation	0.474	1.197	1.176	1.177	0.773
Car service	1.168	0.713	0.089	1.336	1.305
Shopping mall	0.607	1.034	1.657	1.153	0.842
Transportation facilities	0.765	1.022	0.712	1.170	1.015
Residence	0.744	1.139	0.254	0.793	1.266
Culture and media	0.604	0.971	0.872	1.132	1.148

). The regional clustering results are as follows:

Industrial zone (Cluster 1). This zone was mainly located at the edges of the study area and covered a large area. The overall FD value was low (Table 4), with a maximum value of only 42.986, and the distribution of POIs was sparse compared with those in other zones. The FD and CF values of the companies and industries were high. Thus, Cluster 1 is an industrial zone with large numbers of companies and industries and sparsely distributed infrastructure, and has the main production functions of the urban area.

Residential zone (Cluster 2). Cluster 2 had a high FD value (Table 4) of residences, residential buildings were densely distributed in the zone, and commercial and the number of community service POIs were also relatively high. Based on the CF value, all POIs in this zone were more relatively distributed. Therefore, Cluster 2 is a residential zone with densely distributed residential neighborhoods and rich functions.

Central business zone (Cluster 3). Cluster 3 was mainly distributed in the center of the study area, where commercial buildings such as South Street and Wanda Plaza are located. The POI frequency densities of businesses and life services were high in this zone (Table 4). Therefore, Cluster 3 is a central business zone with a high density of facilities and densely distributed business facilities, which take on the main commercial functions of the urban area.

Entertainment zone (Cluster 4). In this zone, the proportions of the hotel, finance, and entertainment POIs were higher than those in other districts. Some residential and living facilities were distributed in this zone, including People’s Park and Xintiandi Plaza. Therefore, this zone is an entertainment zone that undertakes urban leisure and entertainment functions.

Public services zone (Cluster 5). The tourism attraction and governmental organizations POIs within this zone were densely distributed and differed considerably from those in other zones (Table 4). Government buildings such as Changzhou Olympic Sports Center, train stations in Changzhou, and Changzhou library are located in this zone, as well as many restaurants, convenience stores, hotels, and other service establishments. Therefore, Cluster 5 is a public service zone that undertakes urban public service functions and meets the daily needs of people in the city.

The urban function identification results of some typical zones were compared with online maps and high-resolution remote sensing images for verification (

Table 5. Comparisons of the identification results with remote sensing images and Map World.

Identification District	Identification Results	Map World	Remote Sensing Images
South Street
Changzhou Olympic Sports Center
Pubei Community
Xingrong High-tech Company

). The results indicate that the proposed method can effectively identify urban functions in central Changzhou City with high recognition accuracy.

5. Discussion

5.1. Effect of Barriers on the Spread of POIs

To investigate the effect of traffic on the spread of POIs and on urban functional zones, we conducted a comparative experiment without considering the traffic barrier phenomenon. The results with and without road and river barriers are shown in Figure 7. Without roads as barriers, the influence of some public service POIs near Changzhou Railway Station increased the values of the surrounding transportation facilities and factories, thereby causing the community to be incorrectly identified as a public service zone (Figure 8a). Figure 8b shows the Lanling home community. Without rivers as barriers, the interaction between the factory POIs on the northern side of the river and the residential POIs on the southern side caused the industries on the northern side to be incorrectly identified as residential zones, whereas the community on the southern side was identified as a public service zone. Figure 8c shows the Pengxin community. The Love Square area located in the northwest contained a large number of POIs related to entertainment and commerce. Without roads as barriers, the excessive commercial POIs reduced the role of residential POIs in the community, leading to incorrect identification. These examples indicate that the addition of barriers to kernel density analysis can effectively control the overspread of POIs and prevent the POIs of different zones from interacting with each other.

5.2. Stability of the Proposed Method against POI Dilution

To examine the stability of the proposed method, various types of POI points were removed from the total dataset using a random selection method at 20%, 30%, 50%, 60%, 80%, and 90% in equal proportions. The experimental result labels of the deleted data were statistically analyzed to determine their OA and kappa coefficients with the experimental result labels of the source data and the manual judgment of the true labels (

Table 6. Overall accuracy (OA) and kappa coefficient analysis after data dilution.

Percentage Dilution	OA (True Label)	Kappa (True Label)	OA (Experiment Label)	Kappa (Experiment Label)
20%	79.33%	0.74	93.13%	0.79
30%	78.33%	0.73	92.62%	0.77
50%	72.00%	0.68	91.73%	0.75
60%	71.67%	0.65	91.23%	0.73
80%	69.00%	0.61	89.65%	0.69
90%	53.33%	0.41	80.26%	0.40

and Figure 9).

The results indicate that when the data dilution was less than 60%, the urban functional partitioning results did not differ considerably from the experimental results based on the original undiluted POIs, the label matching remained at 90%, and the overall accuracy with the manually identified labels remained above 70%. This indicates that the proposed method is stable when considering POI dilution and does not require a large number of POIs. However, when the amount of POI data fell below a certain threshold, the model could not obtain sufficient information from the POIs to correctly identify urban functions. In this scenario, it would be necessary to supplement the data sources to identify urban functional areas more comprehensively and accurately.

6. Conclusions

In this study, a new method for the identification of urban functional zones based on POI density and MGAE was proposed. The proposed method was applied in a case study of Changzhou City and the results were compared with those obtained using the LDA, TF-IDF, and mSDA methods for validation. The main conclusions are as follows:

• Kernel density analysis is effective for enhancing sparsely and unevenly distributed POIs to identify urban functional zones by obtaining the value of POI density and spreading the influences of POIs from points of origin to their surroundings. Thus, patches that do not contain POIs can also acquire the social features of the surrounding POIs, thereby solving the problem of uneven POI distribution. Moreover, kernel density analysis with barriers performed better in terms of representing the social features of the study area.
• Combining the spatial adjacency of the analysis units and social features can improve the performance of urban functional zone identification. This method addresses the problem of ignoring the functional correlations between units and more comprehensively explores the potential social information of POIs within units. The experiments demonstrate that our model performed better than the LDA, TF-IDF, and mSDA models, and the recognition accuracy was approximately 20% better.

In future work, we will integrate multiple source data of cities into the model. For example, social data such as vehicle trajectory data and mobile phone data reflect human activities to some extent, which have effects on urban functional zone. Moreover, optimizing the methods for quantifying the spillover effect of different social data on urban functional zoning or other applications remains to be explored. In addition, different applications of the results of urban functional zoning should be evaluated. A clear understanding the different types of urban functional zones, the structural patterns among these, and the evolution and current status of urban function will serve to provide guidance for urban governance and planning.