Monitoring the Distribution and Variations of City Size Based on Night-Time Light Remote Sensing: A Case Study in the Yangtze River Delta of China

Abstract

Effectively monitoring the size of a city in real time enables the scientific planning of urban development. Models that utilize the distribution and variations in city size generally use population data as inputs, which cannot be obtained in a timely and rapid manner. However, night-time light (NTL) remote sensing may be an alternative method. A case study was carried out on the Yangtze River Delta (YRD) in China, and the rank–size rule, the law of primate cities, and the Gini coefficient were employed to monitor the variation in city size in the study area. The urban areas extracted based on NTL remote sensing were utilized instead of the traditionally used population data to evaluate the variations in city size from 2012 to 2017. Considering the empiricism and subjectivity of the thresholding method, urban areas were extracted from NTL data combined with the normalized differential vegetation index and land-surface temperature data based on the artificial neural network algorithm. Based on the results, the YRD did not fit the distribution of the primate cities from 2012 to 2017. However, this region satisfied the rank–size rule well, which indicated that the development of medium–small cities was more prominent than that of larger cities, and the dispersed force was larger than the concentrated force. Notably, the city size reached a relatively balanced level in the study area. Further, sensitivity analysis revealed that the relatively low extraction accuracy of urban areas of few small cities had little effect on the results of city size variations. Moreover, the validation of city size computed from statistical population data and its comparison with results calculated based on the statistical data of urban areas aligned with the results of this study, which indicates the rationality and applicability of monitoring the variations in city size using the urban areas extracted from NTL remote sensing instead of population data.

1. Introduction

Dramatic urbanization has occurred globally, and is mainly due to rapid economic development and population growth [1]. According to statistical data, more than 54% of land is populated by humans worldwide. In fact, the total population in urban areas is predicted to exceed 2 billion by 2050 [2]. Rapid urbanization can lead to several environmental issues, including air pollution, urban heat islands, climate change, shortage of resources, and pressure for sustainability worldwide [3,4,5,6]. Rapid urbanization may also result in unhealthy development in urban systems, especially developing countries. For instance, the development of big cities was better than that of small and medium cities. Accordingly, the development of city size in an urban system was recognized to be unbalanced. A comprehensive understanding of variations in city size is thus useful for urban planning and decision making, and monitoring the distribution and variations in city size, which is also helpful for understanding the development model and rational planning of urban systems, is of great significance.

Several popular measurement methods for analyzing the distribution and variations in city size include the rank–size rule [7], the law of the primate city [8], and the Gini coefficient [9]. Briefly, the urban population or urban area is input into the above models [10], and the distribution and variation trends of city size are analyzed based on the calculated indices. Although the most direct and effective monitoring method is population data, there are some problems with this method [11,12,13]. First, the statistical standards for population data vary among countries. For instance, population statistics in China are mainly based on administrative units, as urban and rural areas are divided according to the administrative boundaries of urban districts and towns. Socioeconomic units in other countries are generally based on functional areas, such as metropolitan areas in the United States. Moreover, considerable time and labor are required to obtain population data. Night-time light (NTL) remote sensing, which has emerged in recent years, has enabled new directions for research on city monitoring, as it can provide real-time information about the earth with the same observation standard. Several studies employed remote-sensing data to monitor the urbanization of cities and assess socioeconomic activities [14,15,16]. Nitsch [17] showed that urban area and population have similar distribution trends when measuring urban size distribution. Therefore, in this study, urban areas extracted from NTL remote sensing were used to replace the traditionally used population data in these city size analysis models to monitor the distribution and variations in city size.

Many satellite images of different resolutions can be used to extract urban areas for cities that are suitable for different scenarios. High- and medium-high-resolution images, such as Landsat series [18], are suitable for exploring urbanization at the city scale, but are not suitable for large scales, such as urban and global levels. Mapping urban areas on a large scale requires several high-resolution images with cloud-free images, which also requires significant time and labor for processing. Coarse-resolution satellite images, such as MODIS and NTL data, are suitable for large-scale urban area extraction [1,19,20]. However, NTL data can be utilized to monitor human activities, as it captures city light at night, which separates urban areas from the surrounding suburbs [21,22]. According to prior studies, the NTL value has a highly close connection with population, gross domestic product, and built-up areas [23,24]. Thus, this value is more suitable for selecting NTL data for the analysis of variations in city size compared to other optimal satellite images.

Currently, commonly used NTL data are obtained from two NTL satellites, the Defense Meteorological Satellite Program-Operational Linescan System (DMSP-OLS) and Suomi National Polar-Orbiting Partnership Visible Infrared Imaging Radiometer Suite (NPP-VIIRS). NPP-VIIRS has better quality than DMSP-OLS, as it decreases the problem of the saturation effect that exists in DMSP-OLS [25,26,27]. Thus, NPP-VIIRS is better choice for extracting urban built-up areas. Recently, methods for extracting urban areas based on NPP-VIIRS data mainly include two types: thresholding and classification. The thresholding method has been used by many researchers, as it is simple and easy to conduct [28,29,30,31,32]. Nevertheless, this method is too subjective and empirical to determine the thresholds for extracting urban areas, owing to considerable uncertainty among cities at different development levels [33]. To resolve such limitations of the thresholding method, some researchers have utilized machine learning methods for the extraction of urban areas. For instance, Jing et al. [34] applied four machine learning algorithms to extract urban areas of eastern China. Based on their results, the machine learning methods achieved good accuracy and lower sensitivity than the thresholding method. Xu et al. [35] extracted urban areas from NTL data by utilizing an artificial neural network algorithm, which led to a more accurate quality result than the thresholding method. Hence, extracting urban areas based on machine learning is better than extraction achieved with the traditional thresholding approach. The selection of training samples is another technical problem in urban area extraction using a classification method based on machine learning. Extracting training samples from NTL data alone may lead to the overestimation or underestimation of urban areas. To improve the extraction of urban areas, some researchers have added other remote-sensing data for sample selection, including normalized differential vegetation index (NDVI) and land-surface temperature (LST) data [34,35]. Notably, a better extraction result was achieved when the NTL data were combined with multi-source remote-sensing data to select training samples to extract urban areas based on the machine learning method.

The objectives of this study were to extract urban built-up areas in the study region using a combination of NTL data and NDVI and LST data based on the machine learning algorithm, and to monitor the distribution and variations in city size in the study area by utilizing the models of rank–size rule, the law of primate city, and the Gini coefficient, using extracted urban areas instead of the traditional population data.

2. Methodology

2.1. Study Area

The Yangtze River Delta (YRD) is located in eastern China and is adjacent to the East China Sea (Figure 1). This region has a subtropical monsoon climate, with an annual mean temperature of 16.9 °C. The region also has three provinces, namely Jiangsu, Zhejiang, and Anhui, and one international city, Shanghai city, with a total of 26 cities. The YRD covers an area of 217,700 km² and has a population of 150 million. Owing to its unique geographical location and natural resources, the YRD urban agglomeration has become the most developed and largest economic zone in China, and it is one of six metropolitan agglomerations in the world. Further, this area has experienced tremendous urban expansion over the past few decades.

2.2. Data Collection and Preprocessing

Multiple datasets were used in this study (

Table 1. Data used in this study.

Data	Description	Resolution	Time
NPP VIIRS	Night-time light	500 m/Month
MOD 13A1	Normalized Differential Vegetation Index	500 m/16 Days	2012–2017
MOD 11A2	Land-surface temperature	1 km/8 Days
MOD 44w	Global land water mask	250 m/Year	2012–2017
MCD 12Q1	Land-cover type	500 m/Year	2012–2017
Population		Year	2012–2017
Statistical urban areas		0.01 km2	2012–2017

), including NTL data and four MODIS datasets. Version 1 NPP-VIIRS NTL data, on a monthly basis from 2012 to 2017, were obtained from the National Oceanic and Atmospheric Administration (NOAA) (https://ngdc.noaa.gov/eog/viirs/download_dnb_composites.html (accessed on 11 November 2020)). Additionally, the four MODIS data from the same time series, which were 16-day Normalized Differential Vegetation Index (NDVI) product (MOD13A1), 8-day MODIS land-surface temperature (LST) product (MOD11A2), annual global land water mask product (MOD 44w), and annual MODIS Land Cover Type product (MCD12Q1), were obtained from National Aeronautics and Space Administration (NASA) (https://ladsweb.modaps.eosdis.nasa.gov/search/ (accessed on 17 November 2020)). The NPP-VIIRS, MOD 13A1 NDVI, and MOD 11A2 LST data were used to select samples for urban area extraction for the study area. MCD12Q1 land-cover data were used to determine the optimal combination of coefficients for urban area extraction. The MOD 44w data were employed to remove the light intensity of water bodies from the NPP-VIIRS images. In this study, all remote-sensing data were reprojected to the Albers Conical Equal Area projection and resampled to a resolution of 500 m for consistency with the NTL and land-cover data. The annual NPP-VIIRS and LST data were reproduced by the mean value composition in this study, while the annual NDVI data were reproduced based on the maximum value composition to decrease the influence of cloud contamination [36]. The population and actual urban area data from 2012 to 2017 were obtained from the statistical yearbooks of provinces and cities. The population data were utilized to calculate the indices for city size analysis, and the statistical urban areas were used for comparison with urban areas extracted using our NTL remote-sensing-based method. In this study, the permanent urban population at the end of the year was regarded as the urban population of each city in the study area.

2.3. Monitoring of the Distribution and Variations in City Size

In this study, the rank–size rule, the law of the primate city, and the Gini coefficient were employed to analyze the distribution and variations in city size for the study area. The rank–size rule is usually represented by Zipf’s index, and the law of the primate city is calculated using the urban primacy index. Zipf’s index, the urban primacy index, and the Gini coefficient were calculated from the acreage of the extracted urban areas. The three indices were also computed from the population data to validate the results obtained from the extracted urban areas. If the variation trend of the three indices calculated from both data sources is the same, urban areas extracted from NTL remote sensing can replace the population to analyze the variation in the differences in city size.

2.3.1. Rank–Size Rule

The rank–size rule was utilized to reveal the size distribution in an urban system based on the relationship between the size of cities and the rank of city sizes, which can also reflect the development level of a country or an urban system. In 1949, Zipf provided a theoretical foundation [7], which can be expressed as follows:

$$\lg P_i=\lg P_1-q\lg R_i$$

(1)

where lg is the basic logarithm with base 10, $R_i$ is the rank of city $i$ according to city size among all cities, $P_1$ is the city size of the largest city, $P_i$ is the city size of city $i$, and $q$ is Zipf’s index.

The ideal Zipf index should be 1, indicating balanced urban development within the urban system. However, in general, Zipf’s index is not equal to one. An index greater than 1 indicates that the development of big cities is prominent, while that of small and medium-sized cities is underdeveloped. On the contrary, an index lower than 1 indicates that the development of big cities is not prominent, whereas that of small- and medium-sized cities is better. An increasing trend in Zipf’s index indicates that the concentrated force in the distribution of city size is larger than the dispersed force, which means that the development of city size in large cities is faster than that in small cities. Conversely, a decrease in Zipf’s index indicates that the dispersed force in the urban system is larger than the concentrated force, and the development of city size in small cities is faster than in big cities.

2.3.2. Law of Primate City

The urban primacy index was utilized to describe the extent of concentration of the urban population in the primate city in an urban system, which can largely reflect the development features of the urban system. This index was proposed by Mark [8], which indicates that the size of the largest city is markedly greater than that of other cities in an urban system. The index is equal to the ratio of the population of the largest city to that of the second largest city in the study area, which can be described by the following formula:

$$S_2=P_1/P_2$$

(2)

where $S_2$ is a 2-city index and $P_1$ and $P_2$ are the populations of the largest and second-largest cities, respectively.

According to some researchers, a 2-city index is too one-sided. Thus, the 4-city index and 11-city index were proposed, and can be calculated as follows:

$$S_4=P_1/(P_2+P_3+P_4)$$

(3)

$$S_{11}=2P1/{\displaystyle \sum _{i=2}^{11}Pi}$$

(4)

where $S_4$ and $S_{11}$ are the 4-city index and 11-city indices, respectively, $P_1$ is the population of the largest city in the urban system, and $P_2$ to $P_{11}$ is the population of the second to the eleventh city based on city size.

The threshold of the 2-city index should be 2, while the thresholds of the 4-city and 11-city indices should be 1. The greater the number of cities involved in the calculation of the urban primacy index, the more reliable the results. Accordingly, the 11-city index was selected to compute the urban primacy index for the study area. In this study, the urban primacy index was calculated from the age of the extracted urban areas and the population in the study area.

2.3.3. Gini Coefficient

The Gini coefficient was originally used to measure the income gap between residents in a country or region, as proposed by Gini in 1912 [9]. Marshall applied the Gini model to cities to analyze city development based on differences in size by analyzing the degree of population aggregation on the whole urban agglomeration, and proposed the concept of the city Gini coefficient for the first time, which can be calculated as follows:

$$G=T/2S(n-1)$$

(5)

where n is the total number of cities in the urban system, S is the total population of these cities in the urban system, and T is the sum of the absolute differences in population size among these cities in the urban system.

As an indicator of the relative difference, the Gini coefficient can reflect the balance degree of the population distribution in various cities, ranging from 0 to 1. The closer the Gini coefficient to the value of 0, the more dispersed the city size, while the closer the coefficient to the value of 1, the more concentrated the city size. The value range of the Gini coefficient can be divided into five groups: values below 0.2, 0.2, 0.3, 0.3, 0.4, 0.4, 0.5, and greater than 0.5. These groups represent the absolute balance of city size, relative balance of city size, relative reasonable city size, significant difference in city size, and abnormal disparity, respectively.

2.4. Urban Area Extraction

In this study, NPP-VIIRS NTL, NDVI, and LST data were utilized to extract the urban areas of the study area based on an artificial neural network (ANN) algorithm. The main steps for extracting urban areas are described below. First, the pixels occupied by water bodies were removed using global land water mask data, and potential urban areas were regarded as pixels with NTL values greater than 1, according to Yang et al. [37].

The thresholds of the three datasets utilized to select samples were subsequently determined based on the standard deviation from the mean value. Notably, urban areas were extracted city by city in the study area, as the development and economic levels markedly vary in the study area, and the accuracy of the extracted results may be low if urban areas of different cities are extracted as a whole. The equations are as follows:

$$\left\{\begin{array}{l}{T}_j^{VIIRS}=Mea{n}_j^{VIIRS}+aSt{d}_j^{VIIRS}\\ T_j^{NDVI}=Mea{n}_j^{NDVI}+bSt{d}_j^{NDVI}\\ T_j^{LST}=Mea{n}_j^{LST}+cSt{d}_j^{NDVI}\end{array}\right.$$

(6)

where $T_j^{VIIRS}$, $T_j^{NDVI}$, and $T_j^{LST}$ are the thresholds of the NPP-VIIRS, NDVI, and LST data for city j, respectively; $Mea{n}_j^{VIIRS}$, $Mea{n}_j^{NDVI}$, and $Mea{n}_j^{LST}$ are the mean values of the three datasets for city j; $St{d}_j^{VIIRS}$, $St{d}_j^{NDVI}$, and $St{d}_j^{NDVI}$ are the standard deviations of these data for city j; and a, b, and c are coefficients, which are normally set to −1, −0.5, 0, 0.5, and 1, respectively. Different values of $T_j^{VIIRS}$, $T_j^{NDVI}$, and $T_j^{LST}$ can be obtained according to the values of a, b, and c.

Based on the features of high NTL values, high temperatures, and low vegetation coverage in urban areas, the training samples of urban and non-urban areas were selected based on the following formulas:

$$Sample(i,j)=\left\{\begin{array}{c}VIIRS(i,j)>T_j^{VIIRS}\&NDVI(i,j)<T_j^{NDVI}\&LST(i,j)>T_j^{LST},Urban\\ VIIRS(i,j)<T_j^{VIIRS}\&NDVI(i,j)>T_j^{NDVI}\&LST(i,j)<T_j^{LST},Non-urban\end{array}\right.$$

(7)

where $Sample(i,j)$ represents the selected samples of pixel i of city j, and $VIIRS(i,j)$,$NDVI(i,j)$, and $LST(i,j)$ are the pixel values of the NPP-VIIRS, NDVI, and LST data for pixel i of city j, respectively.

Therefore, there are 125 combinations of coefficients a, b, and c. The optimal set of combinations was selected based on the following standard:

$$\min \left|S_k-S_0\right|\&\max \left\{g\_mea{n}_k\right\}$$

(8)

where $S_k$ and $S_0$ are the number of total pixels of the extracted urban areas based on the NTL data and the number of total pixels of the urban areas in the land-cover data, respectively, and $g\_mea{n}_k$ is the geometric mean precision, which is introduced below.

Based on the optimal set of coefficients confirmed above, the optimal thresholds for the three datasets used to select samples for urban area extraction were determined in the third step. Thereafter, the samples were selected and randomly divided into approximately 10% training samples and 90% test samples. An ANN algorithm was employed to extract the urban areas of the study area based on the training samples.

ANN is a type of machine learning model that has been widely utilized in many studies owing to its self-adapting, self-organizing, and self-learning advantages. A feedforward network structure composed of an input layer, a hidden layer, and an output layer was utilized in this study [35]. The input layer consists of three nodes: NTL, LST, and NDVI. Based on a sensitivity analysis, the number of hidden nodes was set to five. The output nodes have two possible results: urban and non-urban areas. The pixels that met the expectations of urban and non-urban areas were outputted as (1, 0) and (0, 0), respectively.

After training, the test samples were utilized to validate the accuracy of urban area extraction for the study area. In this study, in addition to the Kappa coefficient used, the geometric mean accuracy was also used [38]. The geometric mean accuracy was calculated as follows:

$$g\_mean=\sqrt{\begin{array}{ccc}\frac{d}{m+d}& \times & \frac{d}{n+d}\end{array}}$$

(9)

where $g\_mean$ is the geometric mean accuracy, $d$ is the number of matched urban area pixels between NTL data and land-cover data, $m$ is the number of urban area pixels for which NTL data do not match the land-cover data, and $n$ is the number of urban area pixels for which the land-cover data do not match the NTL data.

3. Results and Discussion

3.1. Urban Expansion of the YRD

3.1.1. Accuracy Assessment of Urban Area Extraction

In this study, the kappa coefficient and geometric mean accuracy were used to validate the accuracy of the urban area extraction of the YRD.

Table 2. Mean accuracy of urban area extraction in the YRD.

Accuracy	2012	2013	2014	2015	2016	2017
g-mean (%)	82.10	84.88	83.06	79.21	83.94	84.10
Kappa	0.793	0.828	0.805	0.761	0.816	0.817

shows the average geometric mean accuracy and kappa coefficient for the YRD. As depicted in the table, the urban area extraction of YRD based on the method proposed in this study achieved high accuracy, with an average geometric mean accuracy and kappa coefficient ranging from 79.21% to 84.88% and 0.761 to 0.828, respectively. The geometric mean accuracy of all years was between 82% and 85%, except for 2015, which had a value of 79.21%. In 2012 and 2015, the Kappa coefficients were 0.793 and 0.761; however, in other years, this value was approximately 0.82. Figure 2 and Figure 3 show the distribution of geometric mean accuracy and Kappa coefficient of urban area extraction for all cities in the YRD from 2012 to 2017. This result indicates that most cities in the YRD have a high geometric mean accuracy and Kappa coefficient for six years, with values above 90% and 0.80, respectively. However, the accuracy of several medium- and small-sized cities, such as Anqing, Chizhou, Xuancheng, and Chuzhou of the Anhui province, which had relatively lower economic development, was relatively low, with a geometric mean accuracy and Kappa coefficient of less than 60% and 0.60, respectively.

3.1.2. Spatial–Temporal Variations of Urban Areas

Figure 4 displays the distribution of the extracted urban areas of the study area from 2012 to 2017 and highlights an obvious spatial discrepancy in urban areas. This finding indicates that most urban areas within the study area were distributed in international metropolises (Shanghai, China), metropolitan circles (Suzhou, Changzhou, and Wuxi, China), and provincial capitals (Nanjing, Ningbo, Hangzhou, and Hefei, China) during the six years. On the contrary, only a small proportion of urban areas were scattered in medium- and small-sized cities. As shown in Figure 5, the acreage of all cities in the study area increased gradually from 2012 to 2017. The urban areas of all cities in the study area extended along the central urban areas to the surrounding suburbs (Figure 4). However, the increase in the urban areas of big cities was small, such as Shanghai, Suzhou, Wuxi, and Nanjing, the urban areas of which increased from 2012 to 2014. In contrast, the extension force of the urban areas of medium- and small-sized cities, such as Xuancheng, Shaoxing, and Nantong, increased in recent years, especially from 2015 to 2017.

3.2. Variations in City Size in the YRD

3.2.1. Variations in Rank–Size

Table 3. Fitting results of rank–size in the YRD from 2012 to 2017 (y represents Zipf’s index; x in the left and right columns represent urban area and urban population, respectively).

Year	Urban Areas		Population
	Fitting Equations	R2	Fitting Equations	R2
2012	y = −1.048x + 3.782	0.796	y = −0.884x + 3.347	0.840
2013	y = −1.034x + 3.784	0.801	y = −0.880x + 3.352	0.838
2014	y = −1.019x + 3.784	0.795	y = −0.875x + 3.357	0.835
2015	y = −1.008x + 3.788	0.798	y = −0.856x + 3.351	0.860
2016	y = −1.010x + 3.793	0.797	y = −0.852x + 3.357	0.860
2017	y = −0.982x + 3.791	0.789	y = −0.845x + 3.361	0.858

shows the fitting results of the double logarithm of the rank–size model in the YRD from 2012 to 2017. These results indicate that both the urban areas extracted from remote-sensing data and the population data can fit the model well, with R² values of approximately 0.79 and greater than 0.83 from 2012 to 2017, respectively. Although the value of R² fitted from urban areas was lower than that fitted from population size, the value met the requirements of the rank–size model. Therefore, a subsequent analysis of the variations in rank size of the YRD from 2012 to 2017 was performed.

Zipf’s index was calculated from urban areas and the population size of the YRD from 2012 to 2017, respectively; the values are shown in Figure 6. The Zipf’s index calculated from the extracted urban areas from 2012 to 2017 was below the value of one, which indicates that the development of big cities in the YRD was not prominent, while that of small- and medium-sized cities was better. In addition, a decreasing Zipf’s index trend calculated from the extracted urban areas was found, with a value of 0.704 in 2012, which decreased to 0.651 in 2017. Such findings indicate that the dispersed force in the YRD was larger than the concentrated force. Thus, the development of the city size of small cities is faster than that of big cities in the YRD. However, from 2010 to 2013, the urban area in the Yangtze River Delta expanded rapidly, while the resident population increased slowly, especially in small and medium-sized cities. Thus, the Zipf’s index calculated from urban area is bigger than that calculated by population.

The Zipf’s index calculated from the population data of the YRD displayed the same distribution as that calculated from the extracted urban areas, which indicates the feasibility of analyzing the variation in rank–size by utilizing NTL remote sensing.

3.2.2. Variations in Primate City

In this study, the 11-city index of urban primacy was utilized to analyze the variations in urban primacy in the study area. Figure 7 shows the distribution of the YRD urban primacy index from 2012 to 2017. As depicted in the figure, the urban primacy index calculated from the extracted urban areas of the YRD was lower than one from 2012 to 2017, which indicates that the YRD agglomeration did not belong to the distribution of the primate city during this period. Additionally, the urban primacy index of the study area decreased gradually, declining from 0.506 in 2012 to 0.471 in 2017. Such a decrease indicates that the size distribution of the primate city in the YRD economic region had a relatively decreasing trend compared with other cities, and the status and function of the primate city in the study area gradually reduced during this period. Such findings indirectly illustrate that the size of medium- and small-sized cities in the study area increased gradually during this period.

The urban primacy index computed from population data also displayed a similar distribution and variation trend to that calculated from the extracted urban areas, which was also below the value of one and decreased gradually from 2012 to 2017. This finding indicates that the application of NTL remote sensing for the analysis of variations in urban primacy is a feasible method.

3.2.3. Variations in the Gini Coefficient

In this study, the Gini coefficient for cities was used to analyze the degree of balance in city development in the study area. Figure 8 shows the variations in the Gini coefficient computed from the extracted urban areas using NTL remote-sensing and population data from 2012 to 2017. The value of the Gini coefficient calculated from the extracted urban areas ranged from 0.2 to 0.3, which indicates that the development of the city size in the study area occurred at a relatively balanced level during this period. Furthermore, the Gini coefficient displayed a gradually decreasing tendency from 2012 to 2017, with the value declining from 0.226 in 2012 to 0.215 in 2017. Such a decrease indicates that the city size of the YRD displayed a scattering trend from 2012 to 2017. Thus, the dispersed force in the study area was greater than the concentrated force, implying that the development of medium- and small-sized cities was faster than that of big cities, and the development of big cities reached saturation in the study area during this period.

The Gini coefficient calculated from the population data was used to validate the value computed from the extracted urban areas of the study area. The value of the Gini coefficient calculated from population data also ranged from 0.2 to 0.3, and displayed a decreasing tendency from 2012 to 2017 (Figure 8), aligning with the variations in the Gini coefficient computed from the extracted urban areas. This result also highlights the feasibility of using extracted urban areas instead of population data to study the distribution and variation in city size.

3.3. Sensitive Analysis

In the accuracy assessment of urban area extraction presented in Section 3.1.1, five small to medium-sized cities with lower levels of economic development, namely Chizhou, Xuancheng, Anqing, Chuzhou, and Yancheng, were found to have relatively low accuracy, with Kappa coefficients and geometric mean accuracy below 60% and 0.6, respectively. This finding was mainly due to cities with lower levels of economic development having relatively lower signals of NTL and LST images, which may affect the urban extraction results. However, the purpose of this study was to prove that urban areas extracted using remote-sensing technology can be utilized to analyze variations in city size instead of statistical population data. Whether the extraction accuracy can affect the results of the city size in the study area must be determined. Therefore, we adjusted the urban areas of the five cities by increasing or decreasing 5% of the extracted urban areas to explore the effect of extraction accuracy on city size monitoring. When the urban areas of one city were adjusted, those of the other four cities remained unchanged. Zipf’s index, the urban primacy index, and the Gini coefficient were computed from the originally extracted urban areas; the urban areas increased by 5%, and urban areas decreased by 5%. The top 11 cities ranked by the acreage of urban areas remained unchanged regardless of whether they decreased or increased by 5%, which caused very slight changes in the Gini coefficient (between 0.0001–0.0003). Thus, only Zipf’s index and the Gini coefficient were calculated (Figure 9a–e) and (Figure 10a–e).

As depicted in Figure 9a–e, the Zipf’s index calculated from the urban areas increasing by 5% or decreasing by 5% was consistent with that computed from the original urban areas, with values below 1.0 from 2012 to 2017, presenting a gradually decreasing tendency. This finding indicates that the development of small and medium-sized cities was more prominent than that of big cities in the YRD, and the dispersed force was larger than the concentrated force in the YRD. Further, as depicted in Figure 10a–e, the Gini coefficient calculated from both urban areas increasing by 5% or decreasing by 5% was consistent with that computed from the original urban areas, with all values ranging from 0.2 to 0.3, highlighting a gradually decreasing tendency from 2012 to 2017. This finding indicates that the development of YRD has reached a relatively balanced level. Further, the city size of YRD tended to be scattered from 2012 to 2017. Regardless of a decrease or increase by 5% in the extracted urban areas, the difference between the two indices calculated from the original urban areas and the adjusted urban areas was very small, with values below 0.001 (the three lines of the two indices calculated from the original and adjusted urban areas were very close). Consequently, the relatively low accuracy of urban area extraction of some cities has little effect on the monitoring of variations in city size in this study.

3.4. Comparison with Other Results

In this study, Zipf’s law, the urban primacy index, and the Gini coefficient were computed from the statistical urban areas, and the distribution and variations in city size were compared with those calculated from the urban areas extracted from NTL remote sensing (Figure 11, Figure 12 and Figure 13). As shown in Figure 11, Zipf’s index calculated from statistical urban areas was consistent with that computed from the urban areas extracted from remote-sensing technology, with values ranging from 0.6 to 0.72 and a gradually decreasing tendency. Moreover, in Figure 12, the urban primacy index calculated from the statistic urban areas showed a decreasing trend, with values ranging from 0.41 to 0.51, thereby aligning with the values calculated from the extracted urban areas. The Gini coefficient calculated from statistical urban areas and extracted urban areas also showed some variation tendency in Figure 13, with values ranging from 0.2 to 0.3. Therefore, all three indices calculated from statistics and extracted urban areas revealed that the YRD did not fit the distribution of the primate city, the dispersed force was stronger than the concentrated force in the study area, and the development of the YRD reached a relatively balanced level. Consequently, the city size analyzed from the extracted urban areas had the same variation trend as that analyzed from the statistical urban areas, which proved the feasibility of utilizing NTL remote sensing to analyze the distribution and variation of city size instead of statistics.

4. Conclusions

The YRD of China was selected as the study area for this assessment and multi-source remote-sensing data based on an ANN algorithm were employed to extract urban areas for the analysis of variations in city size from 2012 to 2017 in the study area, instead of the traditional approach of using population data. The main conclusions are as follows.

By combining NTL data with land-surface temperature and normalized differential vegetation index data, the method employed in this study could be used to efficiently extract the urban areas of the study area based on the ANN algorithm. The geometric mean accuracy and Kappa coefficient ranged from 79.21% to 84.88% and 0.761 to 0.828, respectively. In addition, the urban areas of the study area expanded gradually from 2012 to 2017, with the entire region and all cities in the region undergoing expansion. Moreover, the rank–size rule, urban primacy index, and Gini model were utilized to assess the development of city size in the study area from 2012 to 2017. Based on the results, the study area did not fit the distribution of the primate city, and the status and function of the primate city in the YRD gradually decreased during this period. The development of big cities was not prominent compared with that of small and medium-sized cities, and the dispersed force was greater than the concentrated force in the study area. The development of city size in the YRD urban agglomeration reached a relatively balanced level during this period. Additionally, the sensitivity analysis revealed that the relatively low extraction accuracy of urban areas for a few small cities had little effect on the results of city size in this study. The three indices calculated from the statistical population data and the statistical urban areas also displayed the same distribution and variation tendency as those computed from the extracted urban areas via NTL remote sensing. These findings indicate the applicability and feasibility of analyzing the variation in city size using NTL remote sensing.