Automatic extraction of built-up areas in Chinese urban agglomerations based on the deep learning method using NTL data

Abstract

Rapidly and accurately extracting built-up areas is an essential prerequisite of urbanization research. There have been many studies on the extraction of built-up areas using remote sensing technologies. So far, few studies have been conducted to evaluate the applicability of the deep learning method to extract built-up areas under the condition that only nighttime light (NTL) data are used. This study proposed a deep learning method to extract the built-up areas using NTL data, and applied the method to analyze the spatial and temporal changes of the built-up areas in Chinese two urban agglomerations from 2000 to 2020. The results show that the U-Net deep learning method can be used to extract built-up areas efficiently under the condition that only NTL data are used. The proposed method was able to improve the accuracy of built-up area extraction significantly compared to the existing method. For the extraction of built-up areas in large regions with long time series, the proposed method can facilitate the work and improve the processing efficiency. The gravity centre of the built-up areas in the Central Plains Urban Agglomeration migrated south-eastward, and the gravity centre of the built-up areas in the Shandong Peninsula Urban Agglomeration migrated south-westward, with these gravity centres gradually approaching the geometric centres of the corresponding urban agglomerations. The built-up areas in the Central Plains and Shandong Peninsula Urban Agglomerations grew rapidly, increasing by 4.14 times and 3.73 times from 2000 to 2020, respectively. The built-up areas in the Central Plains Urban Agglomeration expanded faster, while the urban development degree of the Shandong Peninsula Urban Agglomeration was higher. The urban distributions and development modes of these two urban agglomerations were quite different. The Central Plains Urban Agglomeration tended to further agglomerate, while the Shandong Peninsula Urban Agglomeration tended to disperse.

Keywords:

• Built-up area
• nighttime light data
• deep learning
• urban agglomeration

1. Introduction

The built-up area refers to the artificial surfaces dominated by impervious structures such as buildings and roads which are made up of concretes, asphalt, metals, etc (Sharma et al. 2016). The change of built-up areas indicates the speed of urban expansion, and the built-up area distribution demonstrates the urban structure and direction of urban development (Wang et al. 2021). Extracting built-up areas in a timely and accurate manner is critical for optimizing land use patterns and supporting sustainable urban development, management and planning (Feng et al. 2019; Sun et al. 2020; Zheng et al. 2021).

There have been many studies on the extraction of built-up areas using remote sensing technologies (Bhatti and Tripathi 2014; Feng et al. 2019; He et al. 2020; Sun et al. 2020; Li et al. 2021; Wang et al. 2021; Shojaei et al. 2022). According to the data sources used in those studies, there are two main kinds of methods: traditional optical remote sensing data and nighttime light (NTL) data. Bhatti et al. used Landsat 8 Operational Land Imager (OLI) data to improve the extraction accuracy of built-up areas through a modified normalized difference built-up index (NDBI) approach (Bhatti and Tripathi 2014). Liu et al. used ZY3 multiview optical remote sensing imagery to extract 45 urban built-up areas across the world by fusing a series of architectural features (Liu et al. 2019). Wang et al. (2021) used multisource high-spatial-resolution optical remote sensing images to extract the built-up areas in Zhengzhou by a deep learning method and analysed the changes of built-up areas in the city from 2016 to 2020. The extraction of built-up areas by traditional optical remote sensing data is essentially a classification process. A single data point covers a small area, resulting in high extraction costs for large built-up areas. In addition, due to the complex structure and diverse ground objects inside cities, traditional optical remote sensing data are prone to the phenomena of mixed pixels, same objects of different spectra, and different objects of same spectrum, all of which cause some interference to the built-up-area extraction results.

NTL data record nocturnal artificial light on the Earth’s surface and provide unique observations of human activities (Elvidge et al. 1997b; Elvidge et al. 2001). At present, the major sensors used for the remote sensing of NTL include the Defense Meteorological Satellite Program-Operational Linescan System (DMSP-OLS) (Elvidge et al. 1997a), National Polar-orbiting Partnership-Visible Infrared Imaging Radiometer Suite (NPP-VIIRS) (Elvidge et al. 2013a; Elvidge et al. 2013b; Elvidge et al. 2013c), LuoJia1-01 (Li et al. 2018; Li et al. 2019), Jilin-1 (JL1-3B) (Zheng et al. 2018), JL1-07/08 (Zhao et al. 2019), and Earth Remote Observation System-B (EROS-B) sensors (Levin et al. 2014; Katz and Levin 2016). NTL data are applied to estimate gross domestic product (GDP) (Zhao et al. 2017), model carbon dioxide emissions (Shi et al. 2016; Shi et al. 2018), monitor disasters (Zhao et al. 2018), examine ecological light pollution (Horton et al. 2019), and map built-up areas (Zheng et al. 2021). Many studies have shown that NTL data are the best choice for built-up area extraction (Small et al. 2005; Shi et al. 2014; Zhou et al. 2014; Zhou et al. 2015; Li and Zhou 2017; Li et al. 2018; Zhou et al. 2018; Ouyang et al. 2019; Hu et al. 2020; Xu et al. 2020). The main method is the threshold method, in which the best threshold is determined by experience, statistics or clustering, and a threshold is used to divide the built-up area and nonbuilt-up area (Imhoff et al. 1997; Small et al. 2005; He et al. 2006b; Liu et al. 2012; Zhou et al. 2015; Xie and Weng 2016; Zhou et al. 2018; Zhao et al. 2020; Li et al. 2021). The determination of the threshold value is susceptible to image noise and other factors, and the NTL intensity differs among different cities with different economic development levels. The threshold obtained for one city cannot be applied to other cities. Therefore, the threshold method is not a global method and is not applicable for studies of large regions.

Deep learning is a data-driven technology. Through multilayer neural network learning, this method can automatically acquire the deep, nonlinear features of the research object (Yuan et al. 2020). This technology has been widely used in computer vision, speech recognition, natural language processing, object recognition and image classification research (Lecun et al. 2015). Xu et al. (2020) tried to apply the artificial neural network method to extract built-up areas using NTL data, and the research has shown that artificial neural networks can extract built-up areas quickly and effectively.

An urban agglomeration is defined as a group of cities consisting of one or more mega-cities at the core of at least three other large cities within a specific geographical area, all of which rely on a developed transport and communication infrastructure network, have a compact spatial organization, and close economic links, ultimately leading to the formation of conurbations and a high degree of integration (Fang 2015). China has invested both financially and academically to facilitate the formation and growth of urban agglomerations (Fang and Yu 2017). The urban agglomeration is the ultimate urban spatial form for China’s New Urbanization.

In this study, the Central Plains Urban Agglomeration and Shandong Peninsula Urban Agglomeration in the middle and lower reaches of the Yellow River in China were taken as the study areas. The objectives of this study comprise: (1) developing a method that can be applied to NTL data for automatic extraction of built-up areas; and (2) evaluating the applicability of the proposed method to extract built-up areas for large regions with long time series; and (3) analyzing the spatial and temporal changes of the built-up areas in Chinese two urban agglomerations from 2000 to 2020.

2. Related work

2.1. Extraction of built-up areas based on traditional optical remote sensing data

Traditional optical remote sensing data generally refer to remote sensing images obtained by sensors that utilize solar radiation such as Landsat TM, SPOT HRV, and Sentinel-2A. Built-up areas based on traditional optical remote sensing data are generally considered to be composed of impervious surfaces, such asrooftops, roads, sidewalks, parking lots, and other manmade concrete surfaces. The commonly used methods for extracting built-up areas include manual interpretation, statistical index-based methods, classification methods, and deep learning neural networks.

Conventionally, manual interpretation through remote sensing images has been used to extract built-up areas. However, manual interpretation is labor-intensive, prohibitively expensive for large area, and difficult to keep the interpretation results consistent (Zhou and Wang 2008).

The statistical index-based methods establish the index models based on the spectral differences of different land use types. The commonly statistical index-based methods include ISA (Carlson and Arthur 2000), IBI (Xu 2008), NDISI (Xu 2010), BCI (Deng and Wu 2012), MNDIS (Liu et al. 2013), BAEM (Bhatti and Tripathi 2014), NDII (Wang et al. 2015), and CBI (Sun et al. 2016). The statistical index-based methods are usually used to extract the built-up areas on low and medium spatial resolution images. The problems of mixed pixels and spectral similarity affect the applicability of the index methods to some extent. Another challenge for all statistical indices is the determination of the optimal threshold value to classify each pixel as either impervious or pervious (Parekh et al. 2021).

Extensive research has been conducted on the application of classification techniques to extract built-up areas. The maximum likelihood classifier was used to perform land cover classification in medium-sized Indian cities (Chettry 2022). Object-based threshold classification technique was used to identify LULC categories in Ankara City (Cengiz et al. 2022). Support vector machine (SVM) was used to obtain an annual map of impervious surfaces in Wuhan city with Landsat Thematic Mapper (TM), Enhanced Thematic Mapper (ETM+), and Operational Land Imager (OLI) images from 1987 to 2016 (Shi et al. 2017). The classification and regression tree (CART) was used to predict impervious surfaces in Chicago, Venice, and Guangzhou (Wang et al. 2017). Classification methods need large volumes of labelled data for training and evaluation. This requires significant effort for manual labelling of data and visual interpretation of results. As a result, these approaches typically focus on predicting impervious surfaces in specific geographic regions (Parekh et al. 2021). Another challenge for these approaches is deciphering mixed pixels representing a combination of impervious surfaces and other land cover types (Zhang and Du 2015).

Recent years deep learning neural networks have been used successfully to extract built-up areas. The method for automatic extraction of impervious surfaces from high-resolution remote sensing images based on deep learning (AEISHIDL) had higher accuracy and automation level compared with other four representative methods (Huang et al. 2019). The different deep learning neural networks were trained by using TensorFlow to predict impervious surfaces from Landsat 8 images (Parekh et al. 2021). The automatic framework of mapping impervious surface growth with long-term Landsat imagery based on temporal deep learning network was proposed and efficiently reflected the impervious surface growth (Yin et al. 2022). The Small Attention Hybrid Unet (SAH-Unet) network was proved to have superior and accuracy for the extraction of impervious surface information (Chang et al. 2023). Although deep learning neural networks have shown certain advantages in extracting built-up areas, large volumes of training samples and the networks’ generalization have been currently challenging issues in this area of research.

Generally, traditional optical remote sensing data are limited by their geographic coverage; they require a large amount of human and computational resources to extract one-time built-up area information for the whole nation (Liu et al. 2012). Therefore, it is difficult to obtain built-up area information for multiple years at a large spatial scale using these types of data.

2.2. Extraction of built-up areas based on NTL data

Given the advances in NTL satellite sensors and technologies, satellite-observed NTLs have emerged as unique geospatial data products that provide a measure of the lighting brightness observed at night from space (Zhao et al. 2019). NTL satellite sensors are suitable for detecting the dynamic process of built-up areas at a large spatial scale due to its appropriate spatial and temporal resolution. NTL data have been demonstrated to be valuable for extracting built-up areas from global to local scales (Zhao et al. 2020). The methods for extracting built-up areas based on NTL data can be grouped into two major categories.

The first category is the threshold-based method. The threshold is defined as the NTL data DN value above which the pixel is classified as built-up area. A single threshold technique for using DMSP-OLS NTL data was originally proposed to map urban areas in the United States (Imhoff et al. 1997). In a subsequent analysis, many authors recognized the limitations of using a single threshold for a global analysis (Small et al. 2005). Therefore, researchers adopted multiple thresholds among regions (Sutton et al. 2001; Henderson et al. 2003; Small et al. 2005; He et al. 2006b; Liu et al. 2012; Li et al. 2021) and used segmentation-based methods to identity local-optimized thresholds (Zhou et al. 2015; Xie and Weng 2016; Zhou et al. 2018; Zhao et al. 2020). The threshold-based methods have been widely used for extracting built-up areas, but the accuracy can decrease because of the empirical and subjective selection of threshold values.

The second category is based on the supervised classification. (Cao et al. 2009) developed a support vector machine (SVM)-based region-growing algorithm to semi-automatically extract built-up areas from DMSP-OLS NTL data. (Jing et al. 2015) proposed four machine learning methods, including classification and regression tree (CART), k-nearest neighbor (k-NN), random forest (RF), and support vector machine (SVM), to extract built-up areas by using DMSP-OLS NTL data and MODIS data. (He et al. 2019) proposed an FCN-based method to perform detection of global built-up areas from 1992–2016, and proved that the deep learning method based on the supervised classification had great potential to effectively detect global built-up areas. (Xu et al. 2020) developed the artificial neural network method to extract built-up areas using NTL data, which further demonstrated that the deep learning method based on the supervised classification can extract built-up areas quickly and effectively.

3. Study area and data

3.1. Study area

The study area includes two urban agglomerations in the middle and lower reaches of the Yellow River in China: the Central Plains Urban Agglomeration (inland urban agglomeration) and the Shandong Peninsula Urban Agglomeration (coastal urban agglomeration). According to the Central Plains Urban Agglomeration Development Plan approved by the State Council on December 28, 2016, the Central Plains Urban Agglomeration is located between 31.38° N and 37.82° N and between 110.21° E and 118.19° E and includes the whole province of Henan and some cities in Shanxi, Hebei, Shandong and Anhui provinces, covering an area of 287,000 square kilometres (as shown in Figure 1(a)). According to the Development Plan of Shandong Peninsula Urban Agglomeration (2016-2030) published by the Shandong Provincial Government in 2017, Shandong Peninsula Urban Agglomeration includes 17 cities in Shandong Province and has an overall layout of ‘two circles and four districts’ (as shown in Figure 1(b)). This agglomeration is located in the eastern coastal region and thus has obvious regional advantages in building a new pattern of internal and external openings.

Figure 1. The study area: (a) Central Plains Urban Agglomeration and (b) Shandong Peninsula Urban Agglomeration.

3.2. NTL datasets

Commonly used NTL datasets include DMSP-OLS stable NTL data, NPP-VIIRS NTL data, and LuoJia-1 data. DMSP-OLS, which acquires NTL images at a 2.7-km resolution without on-board radiometric calibration, is the most commonly used sensor for NTL research, with composites released annually between 1992 and 2013. These data have disadvantages, including the lack of on-orbit radiance calibration, saturation issues, and blooming issues (Cao et al. 2014; Cao et al. 2019; Levin et al. 2020), that limit their potential applications. NPP-VIIRS NTL data are of higher quality (e.g. have a higher spatial resolution of ∼500 m) and have a superior detection ability, but the short available time span causes problems when long-term analysis is required (Chen et al. 2021). Luojia-1 NTL data have a spatial resolution of 130 m at the nadir point and provide abundant spatial detail information (Li et al. 2019). However, these data have a shorter available time span.

Chen et al. built an extended time series (2000–2018) of NPP-VIIRS-like NTL data through a new cross-sensor calibration from DMSP-OLS NTL data (2000–2012) and a composition of monthly NPP-VIIRS NTL data (2013–2018) (Chen et al. 2021). The extended NPP-VIIRS-like NTL data (2000–2018) have an excellent spatial pattern and temporal consistency. In addition, the resulting product could be easily updated and provide a useful proxy to monitor the dynamics of demographic and socioeconomic activities for longer time periods than those of existing products. The product has been updated until 2020. This study chose the extended NPP-VIIRS-like NTL data (2000-2020) as the research NTL dataset.

3.3. Other data

To obtain the reference data of the built-up areas, the statistical data of the built-up areas for each county were derived from the China Urban Construction Statistical Yearbook and China County Seat Construction Statistical Yearbook. In this study, we used 1:1 million vector data of administrative divisions from the national Geographic Information Resources Directory Service system. To facilitate a comparison between the methods, Landsat 8 OLI imagery on the Central Plains Urban Agglomeration in 2020 was downloaded.

4. Methodology

4.1. Overall process of the proposed method

The overall process of the proposed method is illustrated in Figure 2. As shown in Figure 2, the proposed method involves three steps:

Figure 2. Overall process of the proposed method.

1. Sample set production: The extended NPP-VIIRS-like NTL data of the Central Plains Urban Agglomeration from 2015 to 2018 are selected for sample set production. First, the extended NPP-VIIRS-like NTL data are preprocessed by reprojection and transformation; then, the NTL data of the Central Plains Urban Agglomeration and the NTL data base on the county administration units are acquired by using the vector data of the Central Plains Urban Agglomeration and the vector data based on the county administration units respectively; using the ISDC method (shown in Section 4.2), the NTL data base on the county administration units are classified into built-up areas and nonbuilt-up areas according to the statistical data of the built-up areas; based on the classification image, the patch-based sample sets are acquired with the operation of cropping, enhancement and equalization.

2. Deep neural network architecture: Based on U-Net network, the initial neural network is constructed; then, the network is trained and adjusted by using the training sample sets and validation sample sets; the final neural network for built-up areas is completed through multiple iterations of the above operations.

3. Change analysis of the built-up areas in the two urban agglomerations: The extended NPP-VIIRS-like NTL data for the study area are input into the neural network, and the built-up areas of the two urban agglomerations are automatically extracted; then, the change analysis of the two urban agglomerations’ built-up areas is completed from the perspectives of gravity centre migration and landscape pattern.

4.2. Sample set production

The extended NPP-VIIRS-like NTL data of the Central Plains Urban Agglomeration from 2015 to 2018 were selected for sample set production. Tang (2017) showed that the logarithmic transformation of NPP-VIIRS NTL data can effectively suppress the drastic fluctuations of pixel values in the built-up area centre and enhance the homogeneity of the overall pixel value distribution. Therefore, formula (1) was used to transform the extended NPP-VIIRS-like NTL data: (1) $${\mathrm{ }DN}^{\mathrm{‘}}=\ln \hspace{0.17em}\left(DN+1\right)$$(1) where ${DN}^{‘}$ represents the pixel’s grey value after transformation and $DN$ represents the pixel’s grey value before transformation.

On the basis of He et al. (2006a), the improved statistical data comparison (ISDC) method was used to determine the built-up area and nonbuilt-up area labels of each county. The major steps of the method were as follows:

• Step 1. According to the statistical data of the built-up areas, the number of pixels belonging to the built-up areas in the extended NPP-VIIRS-like NTL data was determined and marked as N.

• Step 2. The pixel grey values of the extended NPP-VIIRS-like NTL data were rearranged into one-dimensional array A from large to small.

• Step 3. The grey value corresponding to the N^th element in array A was taken as threshold B.

• Step 4. The pixel grey values of the extended NPP-VIIRS-like NTL data were compared with threshold B, and then the labels of built-up areas and nonbuilt-up areas were determined.

The image patch size was set to 128 × 128. Since the number of pixels in the built-up area was less than that in the nonbuilt-up area, it is ensured that the central pixels of the image patch belonged to the built-up area, thereby improving the imbalance of the sample set. Considering the small volumes of the sample set, the sample set was processed for data enhancement. A total of 3200 image patches were acquired and divided into training samples and validation samples according to an 8:2 ratio.

4.3. Deep neural network architecture

Compared to traditional optical remote sensing images, the semantic meaning of NTL data is relatively simple, and the semantic information is relatively singular. These features are similar to those of medical images, and the extraction of built-up areas is also a two-classification, semantic image segmentation task. Therefore, U-Net (Ronneberger et al. 2015), which is excellent in the field of medical image segmentation, was selected as the deep neural network architecture for the extraction of built-up areas (as shown in Figure 3). U-Net is a symmetric encoding-decoding network with a U-shaped structure. The contracting path (the left half) is an encoding structure that is composed of four downsampling modules, each of which contains two 3 × 3 convolution layers and a 2 × 2 max pooling layer. The expansive path (the right half) is a decoding structure that is composed of four upsampling modules and one output layer. Each upsampling module contains an upsampling of the feature map followed by a 2 × 2 convolution that halves the number of feature channels, a concatenation with the correspondingly cropped feature map from the contracting path, and two 3 × 3 convolution layers. At the output layer, a 1 × 1 convolution is used to map each 64-component feature vector to the desired number of classes.

Figure 3. The deep neural network architecture applied in this study.

4.4. Indicators for the change analysis of the built-up areas in the two urban agglomerations

1. Gravity centre migration indicators

   The gravity centre of the built-up areas in the urban agglomeration was calculated as the average value of the gravity centres of all built-up areas within the urban agglomeration. By calculating the gravity centre for different periods, the spatial distribution changes of the built-up areas in the urban agglomeration could be reflected. The gravity centre coordinates were calculated using formula (2): (2) $$\begin{array}{c} X_t=\frac{\sum _{i=1}^n{w}_i{x}_i}{\sum _{i=1}^n{w}_i},Y_t=\frac{\sum _{i=1}^n{w}_i{y}_i}{\sum _{i=1}^n{w}_i} \end{array}$$(2) where $X_t$ and $Y_t$ are the gravity centre coordinates of the built-up areas in the urban agglomeration in year t; $x_i$ and $y_i$ are the gravity centre coordinates of the i-th built-up area patch in the urban agglomeration; n indicates the total number of built-up area patches in the urban agglomeration; and $w_i$ indicates the weight, i.e., the area of the i-th built-up area patch divided by the total area of all built-up area patches.

   The gravity centre migration distance and direction can reflect the movement trajectory of the gravity centre and were obtained in this study as shown in formulas (3) and (4), respectively: (3) $$\begin{array}{c}{D}_{t,t+k}=\sqrt{{\left(X_{t+k}-X_t\right)}^2+{\left(Y_{t+k}-Y_t\right)}^2} \end{array}$$(3) (4) $$\begin{array}{c}{\theta }_{t,t+k}=\{\begin{array}{c}\arctan (\frac{Y_{t+k}-Y_t}{X_{t+k}-X_t}), X_{t+k}-X_t>0 \\ \arctan (\frac{Y_{t+k}-Y_t}{X_{t+k}-X_t})+\pi , X_{t+k}-X_t<0,Y_{t+k}-Y_t\ge 0\\ \arctan (\frac{Y_{t+k}-Y_t}{X_{t+k}-X_t})-\pi , X_{t+k}-X_t<0,{\mathrm{Y}}_{t+k}-Y_t<0\\ \frac{\pi }{2}, X_{t+k}-X_t=0,Y_{t+k}-Y_t>0\\ -\frac{\pi }{2}, X_{t+k}-X_t=0,Y_{t+k}-Y_t<0\\ no \mathit{migration}, X_{t+k}-X_t=0,Y_{t+k}-Y_t=0\end{array}\end{array}$$(4)

   where $D_{t,t+k}$ is the migration distance of the gravity centre of the built-up areas in the urban agglomeration from year t to year t + k; ${\theta }_{t,t+k}$ is the migration direction of the gravity centre of the built-up areas in the urban agglomeration from year t to year t + k; and $X_{t+k}$ and $Y_{t+k}$ are the gravity centre coordinates of the built-up areas in the urban agglomeration in year t + k.

2. Landscape pattern indicators

   Landscape pattern indicators can be used to express the characteristics and spatial distribution of the built-up area patches in the urban agglomeration, as shown in Table 1.

Table 1. Descriptions of landscape pattern indicators.

Type	Calculation formula	Value range	Unit	Description
Total patch area	$\mathit{TPA}={\sum }_{i=1}^n{a}_i$	TPA > 0	km^2	The total area of built-up patches
Number of patches	$NP=n$	NP > 0		The number of patches can reflect the fragmentation of the built-up patches and has a good positive correlation with the fragmentation degree of the built-up patches.
Average patch area	$\mathit{APA}=\frac{{\sum }_{i=1}^n{a}_i}{n}$	$\mathit{APA}$>0	km^2	The average patch area can reflect the average size of built-up patches.
Patch density	$PD=\frac{n}{A}\times {10}^8$	PD > 0	per 100 km^2	The patch density refers to the number of patches per 100 km2 and can reflect the density of built-up patches.
Largest patch index	$\mathit{LPI}=\frac{\mathrm{Max}(a_1,a_2,\dots ,a_n)}{A}\times 100$	0 < LPI $\le$100	%	The largest patch index is the percentage of the largest built-up patch area in the total area of the study area, indicating the impact of the largest built-up patch on the study area.
Total edge length	$\mathit{TEL}={\sum }_{i=1}^n{e}_i$	TEL > 0	km	The total edge length of all built-up patches in the study area
Average edge density	$\mathit{AED}=\frac{{\sum }_{i=1}^n{e}_i}{A}\times {10}^8$	AED $\ge$0	km per 100 km^2	The average edge density can reflect the degree of dividing the study area by the built-up patches’ boundary.
Landscape shape index	$\mathit{LSI}=\frac{0.25{\sum }_{i=1}^n{e}_i}{\sqrt{A}}$	LSI ≥ 1		The landscape shape index can reflect the shape change of the built-up patches. The larger the value is, the more complex the patch shape is.
Aggregation index	$AI=\left[\frac{g_{ii}}{\max \to g_{ii}}\right]\times 100$	0 ≤ AI ≤ 100	%	The aggregation index can reflect the connectivity between the built-up patches. The smaller the value is, the more discrete the patches are.

Notes: n represents the total numbers of the built-up patches, $a_i$ represents the area of the i-th built-up patch (${\mathrm{m}}^2),$ A represents the total area of the study area (${\mathrm{m}}^2),$ $\mathrm{Max}(a_1,a_2,\dots ,a_n)$ represents the area of the largest built-up patch in the study area (${\mathrm{m}}^2),$ $e_i$ represents the edge length of the i-th built-up patch ($km$), $g_{ii}$ represents the connectivity number between the built-up patches, and $\max \to g_{ii}$ represents the maximum connectivity number between the built-up patches.

5. Results and discussion

5.1. Evaluation of the deep neural network

Through repeated experiments, we set the batch size to 8 and the learning rate to 1e-4 for deep neural network training. The trained network was evaluated with the validation data. Precision, recall, F1 score and intersection over union (IoU) were used as the evaluation indicators, as shown in formulas (5), (6), (7) and (8), respectively. The definitions of the parameters in the formulas are shown in Table 2. TP (True Positive) refers to the number of pixels for which the real categorization is a built-up area and the category prediction by the deep neural network was also the built-up area. FP (False Positive) refers to the number of pixels for which the real category is a nonbuilt-up area and the category predicted by the deep neural network was a built-up area. FN (False Negative) refers to the number of pixels for which the real category was the built-up area and the category predicted by the deep neural network was a nonbuilt-up area. TN (True Negative) refers to the number of pixels for which the real category was the nonbuilt-up area and the category predicted by the deep neural network was also a nonbuilt-up area. (5) $$\begin{array}{c} \mathit{Precision}=\frac{TP}{TP+FP} \end{array}$$(5) (6) $$\mathit{Recall}=\frac{TP}{TP+FN}\mathrm{ }$$(6) (7) $$F1=\frac{\mathit{Precision}\times \mathit{Recall}}{\mathit{Precision}+\mathit{Recall}}\times 2$$(7) (8) $$\begin{array}{c}IoU=\frac{TP}{TP+FN+FP} \end{array}$$(8)

Table 2. The confusion matrix.

	Built-up area (ground truth)	Nonbuilt-up area (ground truth)
Built-up area (prediction)	TP	FP
Nonbuilt-up area (Prediction)	FN	TN

The precision, recall, F1 score and IoU of the trained network were found to be 0.9692, 0.9777, 0.9734 and 0.9483, respectively, on the validation data; these values all meet the accuracy requirements.

5.2. Comparison analysis

In order to prove the advantages of our method, the built-up area extraction method (BAEM) (Bhatti and Tripathi 2014), which has been cited 146 times according to the statistics of Web of Science until now, was selected as a comparison method. The BAEM method using Landsat-8 OLI imagery comprised four major steps: preprocessing and examination of satellite data, image enhancement through resolution merging, producing the BAEM_OLI imagery, and finding the optimal threshold value to segregate built-up from non-built-up areas in the BAEM_OLI imagery. Our method and BAEM method were applied simultaneously to extract the built-up areas in the Central Plains Urban Agglomeration in 2020, and compared in terms of accuracy and operational efficiency. The reference data of the built-up areas were acquired by the ISDC method (shown in Section 4.2).

5.2.1. Accuracy

As shown in Figure 4, the extraction results of our method and BAEM method were both in good agreement with the reference data in most counties. This demonstrated that both methods had high accuracy of extracting the built-up areas.

Figure 4. Extraction results of the built-up areas in the Central Plains Urban Agglomeration in 2020: (a) the built-up areas extracted by our method overlaid with the reference data, (b) the reference data overlaid with the built-up areas extracted by our method, (c) the built-up areas extracted by BAEM method overlaid with the reference data, and (d) the reference data overlaid with the built-up areas extracted by BAEM method.

Comparing Figure 4(a) and (b), the mismatches on the ranges of the built-up areas between our method and the reference data were mainly reflected in County 0 to 12. In County 0 to 3, our method had a smaller range of built-up areas. Through field investigations, we found that in County 0 to 3, the pixels that were not identified as built-up areas by our method while belonged to built-up areas in the reference data located in the urban suburbs. Although these urban suburbs were composed of buildings, concretes, and asphalt, the social and economic activities of human beings were infrequent, and people were accustomed to calling these regions ‘ghost towns’. We preferred not to classify these regions as built-up areas. However, (Zhao et al. 2020) emphasized that the pixels that belonged to built-up areas in the urban suburbs with low intensity of human activity were easy to be ignored because of the effect of shrink of NTL. In order to further understand the reasons, the multisource data could be integrated to improve the proposed method in future research. In County 4 to 12, our method had a bigger range of built-up areas. Through field investigations, we found that in County 4 to 12, the pixels that were identified as built-up areas by our method while did not belong to built-up areas in the reference data located in the towns. With the implementation of China’s rural revitalization policy, the economic development of these towns had been rapid, and the scale of the built-up areas had been expanding. However, due to the lack of statistical data and lagging management, the built-up areas of these towns were generally not included in the statistics. This kind of built-up areas can be well recognized by our method.

Comparing Figure 4(c) and (d), the mismatches on the ranges of the built-up areas between BAEM method and the reference data were reflected in County 0 to 2, County 5 to 7, County 9, and County 13 to 54. In County 0 to 2 and County 13 to 44, BAEM method had a smaller range of built-up areas. The built-up areas of County 13 had not even been extracted at all. Through field investigations and the results from the methods of NDII (Wang et al. 2015) and CBI (Sun et al. 2016), we verified that the omission built-up areas by BAEM method in the above counties should be identified as built-up areas really. In County 5 to 7, County 9, and County 45 to 54, BAEM method had a bigger range of built-up areas. Through field investigations, we found that in these counties, the pixels that were identified as built-up areas by BAEM method while did not belong to built-up areas in the reference data mainly located in the urban suburbs. As mentioned in the previous paragraph, these urban suburbs are known as ‘ghost towns’ due to their infrequent human socio-economic activities. From a socio-economic perspective, it is reasonable to classify these urban suburbs as nonbuilt-up areas. Therefore, BAEM method has certain commission error in these counties.

As shown in Table 3, compared to BAEM method, the precision, recall, F1 score and IoU of our method increased by 7.19%, 0.30%, 3.68% and 5.78%, respectively. This shows that our method has a higher accuracy.

Table 3. Comparison of different methods.

Method	Precision	Recall	F1 score	IoU
Our method	0.9160	0.8642	0.8894	0.8008
BAEM method	0.8441	0.8612	0.8526	0.7430

5.2.2. Operational efficiency

The main hardware configuration of our operating computer was as follows:

• ① Central Processing Unit (CPU): Intel i7-4790 3.60 GHz

• ② Video adapter: NVIDIA GeForce GT 720

• ③ Random Access Memory (RAM): Hynix/Hyundai 8GB

We calculated the running time of our method and BAEM method for extracting the built-up areas in the Central Plains Urban Agglomeration in 2020 respectively.

1. Running time of our method

   Our method had a running time of 59 seconds, including: 24 seconds for model prediction; 35 seconds for prediction results post-processing. For the work of extracting the built-up areas in the Central Plains Urban Agglomeration covering an area of 287,000 square kilometres, a 59-second running time is pretty short.

2. Running time of BAEM method

   BAEM method had a cumulative running time of 8.3 hours, including: 6.4 hours for downloading, preprocessing and examination of Landsat-8 OLI imagery; 0.3 hour for image enhancement; 0.5 hour for producing the BAEM_OLI imagery; 1.1 hours for finding the optimal threshold value to segregate built-up from non-built-up areas in the BAEM_OLI imagery and producing the extraction results. Due to the influence of weather and imaging conditions, high-quality Landsat-8 OLI imageries covering the large regions were difficult to be obtained. This resulted in spending more time in the first stage. For the large region like the Central Plains Urban Agglomeration, finding the optimal threshold value is quite time-consuming. As a result, a considerable amount of time was spent in the final stage.

Comparing the running time of our method and BAEM method, we found that the operational efficiency of our method was significantly higher than that of BAEM method. Moreover, our method requires no human intervention, and the entire operation process is fully automated.

5.3. Change analysis of the built-up areas in the two urban agglomerations

The built-up areas in the Central Plains Urban Agglomeration and Shandong Peninsula Urban Agglomeration were extracted from 2000 to 2020 based on our method. For convenience, the change results of the built-up areas were displayed in a 5-year cycle, as shown in Figure 5. The built-up areas in the two urban agglomerations grew rapidly in the past 20 years, and the growth rates of individual cities were relatively large.

Figure 5. Changes in the built-up areas in the study area from 2000 to 2020: (a) Central Plains Urban Agglomeration and (b) Shandong Peninsula Urban Agglomeration.

Tables 4 and 5 show the statistics of the gravity centre migrations of the built-up areas in the two urban agglomerations. Before 2015, the gravity centre migration distance in the Central Plains Urban Agglomeration was large, indicating that the expansion of the built-up areas in this region was unbalanced during this period, and the unidirectional expansion speed was fast. After 2015, the gravity centre migration distance in the Central Plains Urban Agglomeration was small, indicating that the expansion of built-up areas in this region was relatively balanced during this period. The gravity centre migration distance in the Shandong Peninsula Urban Agglomeration was disorderly and fluctuated greatly.

Table 4. The statistics of the gravity centre migration of the built-up areas in the Central Plains Urban Agglomeration.

Year	X coordinate (m)	Y coordinate (m)	Migration distance (km)	Migration direction (degree)	Remarks
2000	813287.41	3750585.42			The gravity centre of the built-up areas in the previous year was set as the origin of the coordinate axis, and the following coordinate system of the migration direction was established. Notes: The Y axis is the north direction, and the X axis is the east direction.
2001	809663.72	3765045.85	14.91	104.07
2002	805639.48	3748628.10	16.90	−103.77
2003	812316.54	3740823.85	10.27	−49.45
2004	809719.51	3742417.02	3.05	148.47
2005	803127.77	3738681.01	7.58	−150.46
2006	811778.70	3735797.50	9.12	−18.43
2007	821749.30	3726795.12	13.43	−42.08
2008	815813.28	3723487.25	6.80	−150.87
2009	821691.83	3737095.59	14.82	66.64
2010	826950.46	3725605.76	12.64	−65.41
2011	829890.45	3721627.84	4.95	−53.53
2012	829678.92	3727594.77	5.97	92.03
2013	829521.63	3716064.03	11.53	−90.78
2014	837365.63	3716242.52	7.85	1.30
2015	841705.35	3714279.95	4.76	−24.33
2016	842490.85	3713648.03	1.01	−38.82
2017	841667.76	3712200.30	1.67	−119.62
2018	841869.40	3711005.38	1.21	−80.42
2019	840856.85	3712657.74	1.94	121.50
2020	839869.74	3713138.12	1.10	154.05

Table 5. The statistics of the gravity centre migration of the built-up areas in the Shandong Peninsula Urban Agglomeration.

Year	X coordinate (m)	Y coordinate (m)	Migration distance (km)	Migration direction (degree)	Remarks
2000	1184962.75	3966467.83			The gravity centre of the built-up areas in the previous year was set as the origin of the coordinate axis, and the following coordinate system of the migration direction was established. Notes: The Y axis is the north direction, and the X axis is the east direction.
2001	1187833.72	3970144.15	4.66	52.01
2002	1180253.30	3964134.97	9.67	−141.60
2003	1187336.91	3966743.51	7.55	20.22
2004	1182284.83	3961566.14	7.23	−134.30
2005	1183046.87	3962406.31	1.13	47.79
2006	1189501.68	3961308.70	6.55	−9.65
2007	1192738.75	3960306.47	3.39	−17.20
2008	1193096.84	3961502.66	1.25	73.33
2009	1184525.51	3957356.73	9.52	−154.19
2010	1183208.81	3958405.37	1.68	141.47
2011	1182604.38	3957872.33	0.81	−138.59
2012	1174984.17	3955409.34	8.01	−162.09
2013	1177924.01	3956889.60	3.29	26.73
2014	1176144.51	3955113.99	2.51	−135.06
2015	1177952.45	3955068.48	1.81	−1.44
2016	1177911.58	3955161.81	0.10	113.65
2017	1174734.31	3953310.67	3.68	−149.77
2018	1176158.91	3953562.81	1.45	10.04
2019	1172289.24	3953400.81	3.87	−177.60
2020	1172814.07	3952513.07	1.03	−59.41

The gravity centre migration trajectories of the built-up areas in the Central Plains Urban Agglomeration and Shandong Peninsula Urban Agglomeration are shown in Figure 6. From 2000 to 2020, the gravity centre of the built-up areas in the Central Plains Urban Agglomeration was located in the metropolitan area of Zhengzhou, close to the geometric centre of the agglomeration. In 2000, the gravity centre was located in Yuanyang County, Xinxiang City. In 2003, it moved to Zhongmou County, Zhengzhou City. In 2010, it moved to the Kaifeng municipal district. In terms of the migration directions, it migrated 8 times south-eastward, 5 times south-westward, 2 times north-eastward, and 5 times north-westward. Compared to the gravity centre in 2000, it ultimately migrated south-eastward. It can be judged that the northern cities with relatively high development degrees had slower expansion speeds, while the southern cities with lower development degrees had faster expansion speeds. The cities in the eastern region expanded faster than those in the western region. From 2000 to 2020, the gravity centre of the built-up areas in the Shandong Peninsula Urban Agglomeration was located in the Qingdao metropolitan area. Excluding 2001 (located in Qingzhou City) and 2002 (located at the junction of Qingzhou City and Linqu County), the gravity centre was located in Weifang City Linqu County. In terms of the migration directions, it migrated 4 times south-eastward, 8 times south-westward, 6 times north-eastward, and 2 times north-westward. From 2000 to 2007, the gravity centre shifted south-eastward, and from 2008 to 2020, it shifted south-westward. Compared to the gravity centre in 2000, it finally migrated south-westward. It can be judged that the northern cities with relatively high development degrees had slower expansion speeds, while the southern cities with lower development degrees had faster expansion speeds. Before 2007, the eastern cities expanded faster than the western cities, and after 2007, the western cities expanded faster than the eastern cities.

Figure 6. The gravity centre migration trajectories of the built-up areas in the study area from 2000 to 2020: (a) Central Plains Urban Agglomeration and (b) Shandong Peninsula Urban Agglomeration.

Table 6 shows the calculation results of the landscape pattern indicators of the built-up areas in the Central Plains Urban Agglomeration. The TPA increased rapidly from 2000 to 2020, with an average annual growth rate of 8.53% and a total increase of 4.14 times. This shows that the scope of the built-up area expanded rapidly over the past 21 years. The change trends of NP and PD were the same and increased at average annual growth rates of 7.14% and 7.42%, respectively. The NP and PD increased significantly in 2012 and 2013, which indicates that a large number of emerging towns appeared during this period, and the agglomeration entered the stage of high-speed development. After 2013, these two indicators fell back, which indicates that emerging cities and towns began to integrate and connect with large urban patches.

Table 6. The calculation results of the landscape pattern indicators of the built-up areas in the Central Plains Urban Agglomeration.

Year	TPA (km^2)	NP	PD (per 100 km^2)	LPI (%)	TEL (km)	AED (km per 100 km^2)	LSI	AI (%)
2000	1665.50	328	0.11	0.0314	3236.0	1.132	19.7378	76.6489
2001	1847.50	336	0.12	0.0490	3373.0	1.180	19.6163	78.0805
2002	2599.00	414	0.14	0.0680	4576.0	1.601	22.4657	78.7303
2003	2712.00	396	0.14	0.0781	4679.0	1.637	22.3971	79.1874
2004	3558.00	464	0.16	0.0767	6154.0	2.153	25.7741	79.0221
2005	3508.00	442	0.15	0.0838	5852.0	2.047	24.7089	79.8074
2006	3638.75	433	0.15	0.0984	5743.5	2.009	23.7479	80.9304
2007	3577.75	474	0.17	0.0868	6005.0	2.101	25.0333	79.6773
2008	3787.75	471	0.16	0.1125	6174.0	2.160	24.9960	80.2795
2009	3736.00	509	0.18	0.0836	6451.0	2.257	26.3306	79.0642
2010	4203.50	617	0.22	0.1014	7639.5	2.672	29.3923	77.8770
2011	4976.50	654	0.23	0.1438	7997.0	2.798	28.2933	80.4599
2012	4944.50	905	0.32	0.0937	8890.5	3.110	31.5496	78.0644
2013	6387.25	1240	0.43	0.1374	11448.0	4.005	35.8563	78.0338
2014	6730.00	1125	0.39	0.1439	11164.5	3.906	34.0213	79.6976
2015	6771.75	1092	0.38	0.1445	11146.5	3.899	33.8424	79.8715
2016	6827.00	1095	0.38	0.1474	11417.0	3.994	34.5559	79.5395
2017	7403.25	1137	0.40	0.1795	11955.5	4.182	34.7217	80.2415
2018	7771.00	1166	0.41	0.1957	12519.5	4.380	35.5411	80.2750
2019	8198.00	1211	0.42	0.1867	13049.5	4.565	36.0165	80.5109
2020	8568.25	1303	0.46	0.2043	13836.0	4.840	37.3450	80.2215
Annual growth rate	8.53%	7.14%	7.42%	9.82%	7.54%	7.54%	3.24%	8.53%

The average annual growth rate of the LPI was 9.82%, which was higher than that of the TPA. This shows that the Zhengzhou metropolitan area expands faster than all built-up areas in the Central Plains Urban Agglomeration. Moreover, the LPI showed an increasing trend year by year, which indicates that the influence of the Zhengzhou metropolitan area in the Central Plains Urban Agglomeration is constantly increasing. As the development core of the networked spatial pattern in the Central Plains Urban Agglomeration, the Zhengzhou metropolitan area has played a good radiation role, fulfilled the responsibility of ‘one core’, achieved deep integration with surrounding cities, and further promoted the development of peripheral cities.

The change trends of the TEL, AED and LSI increased each year, with average annual growth rates of 7.54%, 7.54% and 3.24%, respectively. This shows that the shape complexity of the built-up areas in the Central Plains Urban Agglomeration was gradually increasing, and the scope of the built-up areas was gradually expanding. Taking 2014 as the dividing point, the average growth rate of the LSI from 2000 to 2013 was 4.70%, and the average growth rate from 2014 to 2020 was 1.57%. These results indicate that with the passage of time, the increase rate of the shape complexity of the built-up areas in the Central Plains Urban Agglomeration decelerated, and the overall expansion speed in the latter stage was slower than that in the previous stage. Before 2014, the data fluctuated greatly, showing that the expansion of the built-up areas in this stage was not stable from the perspective of spatial characteristics.

The AI generally rose from 2000 to 2020. This shows that the fragmentation degree of the built-up areas in the urban agglomeration tended to be low, and the links between cities and towns became increasingly close throughout the study period. However, the AI fluctuated violently and showed a trend of rising, falling and rising. From 2001 to 2014, the fluctuations were extremely sharp, indicating that the urban expansion in this stage was unstable. From 2015 to 2020, the fluctuations decelerated, indicating that the urban expansion in this stage tended to be stable, and the Central Plains Urban Agglomeration gradually transitioned from high-speed development to high-quality development.

According to the development plan of the Central Plains Urban Agglomeration, the agglomeration can be divided into six main regions: the Zhengzhou Metropolitan Region, Core Development Region, Eastern Industrial Transfer Demonstration Region, Southern High-efficiency Ecological Economy Demonstration Region, Western Transformation and Innovation Development Demonstration Region and Northern Cross-Regional Collaborative Development Demonstration Region (Figure 1(a)). Considering a five-years cycle, five landscape pattern indicators, including the TPA, NP, APA, LSI and AI, were adopted to analyse the status of the built-up areas in these six main regions (Figure 7). Figure 7 shows that the status of the built-up areas in various regions exhibited obvious differences. The APA and AI of the Zhengzhou Metropolitan Region had the highest values, reflecting the highest development degree of the region. The five indicators of the Western Transformation and Innovation Development Demonstration Region were not high, and this region was relatively underdeveloped overall. Figure 7(a) (b) and (d) shows that the TPA, NP and LSI of the Core Development Region and Northern Cross-Regional Collaborative Development Demonstration Region were relatively high. This reflects the large-scale and rapid development of the built-up areas in these two regions, providing is an important basis for supporting the development of the Central Plains Urban Agglomeration. Although the development speed and scale of the Eastern Industrial Transfer Demonstration Region and Southern High-efficiency Ecological Economy Demonstration Region were not as fast as those of the Core Development Region and Northern Cross-Regional Collaborative Development Demonstration Region, the APA and AI of these two regions have performed well in recent years, reflecting the obvious trend of the integration of the built-up areas in these two regions.

Figure 7. The changes in the typical landscape pattern indicators in the six main regions of the Central Plains Urban Agglomeration: (a) TPA, (b) NP, (c) APA, (d) LSI, and (e) AI.

Table 7 shows the calculation results of the landscape pattern indicators of the built-up areas in the Shandong Peninsula Urban Agglomeration. The TPA increased rapidly from 2000 to 2020, with an average annual growth rate of 6.8% and a total increase of 3.73 times. The change trends of the NP and PD were the same and increased at average annual growth rates of 8.82% and 8.86%, respectively. The Shandong Peninsula Urban Agglomeration is a dual-core urban agglomeration containing Jinan and Qingdao. The built-up areas of these two cities alternately became the largest patches, causing the LPI of the urban agglomeration to fluctuate, but the overall trend was upwards, with an average annual growth rate of 5.7%.

Table 7. The calculation results of the landscape pattern indicators of the built-up areas in the Shandong Peninsula Urban Agglomeration.

Year	TPA (km^2)	NP	PD (per 100 km^2)	LPI (%)	TEL (km)	AED (km per 100 km^2)	LSI	AI (%)
2000	1828.00	205	0.13	0.0906	2911.0	1.868	17.1860	80.7362
2001	2149.75	250	0.16	0.1126	3442.0	2.208	18.7634	80.5784
2002	2463.00	314	0.20	0.1030	4133.0	2.652	21.0804	79.5129
2003	2844.50	338	0.22	0.1163	4645.0	2.980	22.0280	80.0373
2004	3348.75	386	0.25	0.1156	5651.5	3.626	24.7328	79.2680
2005	3467.50	367	0.24	0.1193	5584.5	3.583	23.9873	80.2756
2006	3597.75	385	0.25	0.1678	5775.0	3.705	24.3958	80.3272
2007	3735.50	334	0.21	0.1647	5425.5	3.481	22.4286	82.2869
2008	3866.50	432	0.28	0.1517	6269.5	4.022	25.5783	80.0541
2009	3901.00	351	0.23	0.1862	5711.5	3.664	23.1560	82.1080
2010	4166.00	580	0.37	0.1468	7703.5	4.942	30.1506	77.1689
2011	4203.75	620	0.40	0.1752	7463.5	4.788	29.0538	78.1420
2012	4316.50	572	0.37	0.1296	7401.5	4.749	28.5437	78.8614
2013	5127.75	957	0.61	0.1660	9603.5	6.161	33.9617	76.7767
2014	5510.75	911	0.58	0.1803	9753.0	6.257	33.3165	78.0813
2015	5530.25	947	0.61	0.1788	9928.5	6.370	33.8154	77.7467
2016	5632.25	940	0.60	0.1894	9856.0	6.323	33.2326	78.3229
2017	6052.25	950	0.61	0.2254	10243.0	6.572	33.3141	79.0421
2018	6288.75	1033	0.66	0.2177	10818.0	6.941	34.4937	78.6946
2019	6441.00	1035	0.66	0.2523	11066.0	7.100	34.7733	78.7623
2020	6814.00	1112	0.71	0.2743	11724.5	7.522	35.8248	78.7250
Annual growth rate	6.8%	8.82%	8.86%	5.7%	7.21%	7.21%	3.74%	−0.13%

The change trends of the TEL, AED and LSI increased each year, with average annual growth rates of 7.21%, 7.21% and 3.74%, respectively. These results show that the shape complexity of the built-up areas in the Shandong Peninsula Urban Agglomeration was gradually increasing, and the scope of the built-up areas was gradually expanding. The average growth rate of the LSI was 5.37% before 2013 and 0.89% after 2013. This shows that the growth rate of the shape complexity of the built-up areas in the Shandong Peninsula Urban Agglomeration decelerated over time. Moreover, before 2013, the TEL, AED and LSI values fluctuated greatly, indicating that the urban expansion at this stage was not stable from the perspective of spatial characteristics and that the scope of the built-up areas was in a period of rapid expansion.

The AI fluctuated sharply and showed a trend of rising first, then falling and then rising again. However, from 2000 to 2020, the overall trend was downwards. This shows that the fragmentation degree of the built-up areas in the Shandong Peninsula Urban Agglomeration was rising, and the spatial pattern of the built-up areas was more scattered than the original. The fluctuations between 2001 and 2014 were extremely sharp, and this further explained the instability of urban expansion at this stage. From 2015 to 2020, the fluctuations decelerated and an overall upwards trend was observed, indicating that urban expansion was gradually stable and tended towards agglomeration during this stage.

In the Development Plan of Shandong Peninsula Urban Agglomeration (2016-2030), the Shandong Peninsula Urban Agglomeration is divided into six major areas, namely, the Jinan Metropolitan Region, Qingdao Metropolitan Region, Yanwei Metropolitan District, Dongbin Metropolitan District, Jizaohe Metropolitan District and Linri Metropolitan District, and the overall pattern of ‘two regions and four districts, network development’ could be constructed. Figure 8 shows the changes in the TPA, NP, APA, LSI and AI values in these six areas from 2000 to 2020. From Figure 8(a) (b) and (d), it can be seen that the six areas were in a rapid development trend since 2000, and among them, the Jinan Metropolitan Region and Qingdao Metropolitan Region developed fastest. As shown in Figure 8(c) and (e), the APA and AI of the six areas fluctuated greatly and showed downwards trends. These results indicate that the six areas were still in a decentralized development state and that the integration and interconnection of the urban built-up areas was not sufficient.

Figure 8. The changes in the landscape pattern indicators of the six areas in the Shandong Peninsula Urban Agglomeration from 2000 to 2020: (a) TPA, (b) NP, (c) APA, (d) LSI, and (e) AI.

6. Conclusions

This study proposed a deep learning method to extract the built-up areas using NTL data, and applied the method to analyze the spatial and temporal changes of the built-up areas in Chinese two urban agglomerations from 2000 to 2020. The built-up areas for extracting by the proposed method were in good agreement with the reference data. Compared with BAEM method, our method had a higher accuracy and higher operational efficiency.

During the past few decades, the gravity centre of the built-up areas in the Central Plains Urban Agglomeration migrated south-eastward, and the gravity centre of the built-up areas in the Shandong Peninsula Urban Agglomeration migrated south-westward, with these gravity centres gradually approaching the geometric centres of the corresponding urban agglomerations. The built-up areas in the Central Plains and Shandong Peninsula Urban Agglomerations grew rapidly, increasing by 4.14 times and 3.73 times from 2000 to 2020, respectively. The built-up areas in the Central Plains Urban Agglomeration expanded faster, while the urban development degree of the Shandong Peninsula Urban Agglomeration was higher. The urban distributions and development modes of these two urban agglomerations were quite different. The Central Plains Urban Agglomeration tended to further agglomerate, while the Shandong Peninsula Urban Agglomeration tended to disperse.

This study opens new research avenues for the long-term and reliable monitoring of built-up areas in large regions using NTL data. The proposed method shows great potentials in extracting the built-up areas. In future research, multisource data, such as Landsat OLI, land surface temperature (LST), normalized difference vegetation index (NDVI), and point-of-interest (POI) data, could be integrated to improve the accuracy of the proposed method. In terms of NTL data, other NTL data could also be used, such as LuoJia1-01 Jilin-1 data.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

Henan Provincial Science and Technology Research Project; National Natural Science Foundation of China;