Monitoring the Spatiotemporal Trajectory of Urban Area Hotspots Using the SVM Regression Method Based on NPP-VIIRS Imagery

: Urban area hotspots are considered to be an ideal proxy for spatial heterogeneity of human activity, which is vulnerable to urban expansion. Nighttime light (NTL) images have been extensively employed in monitoring current urbanization dynamics. However, the existing studies related to NTL images mainly concern detection of urban areas, leaving inner spatial differences in urban NTL luminosity poorly explored. In this study, we propose an innovative approach to explore the spatiotemporal trajectory of urban area hotspots using monthly Visible Infrared Imaging Radiometer Suite (VIIRS) NTL images. Firstly, multi-temporal VIIRS NTL intensity was decomposed by time-series analysis to obtain annual stable components after data preprocessing. Secondly, the support vector machine (SVM) regression model was utilized to identify urban area hotspots. In order to ensure the model accuracy, the grid search and cross-validation method was integrated to achieve the optimized model parameters. Finally, we analyzed the spatiotemporal migration trajectory of urban area hotspots by the center of gravity method (i.e., shift distance and angle of urban area hotspot centroid). The results indicate that our method successfully captured urban area hotspots with a regression coefﬁcient over 0.8. Meanwhile, the ﬁndings give an intuitive understanding of coupling interaction between urban area hotspots and socioeconomic indicators. This study provides important insights for further decision-making regarding sustainable urban planning.


Introduction
Urbanization is a global concern associated with massive population shift, urban expansion, and industrial structure adjustment [1][2][3]. Especially as a developing country, China has experienced rapid urbanization since the reform and opening-up policy was implemented [4,5]. The rapid urban expansion led to intensive population increases and regional development differences. Therefore, timely and comprehensive understanding of urban growth is vital for sustainable urbanization. With advances in remote sensing and Geographic Information Systems (GIS) technology, satellite imagery provides convenient and spatially explicit ways of mapping urban expansion. Compared with other remote sensing products, the artificial nighttime light (NTL) image, as a recorder of anthropogenic luminosity [6][7][8], provides a more specific perspective of urban dynamics at different scales (i.e., local, regional, and global).
To date, two types of widely used types of NTL imagery at the global scale can be acquired from different satellite sensors: the Defense Meteorological Satellite Program/Operational Linescan System (DMSP/OLS) and the Visible Infrared Imaging Radiometer Suite (VIIRS) Day/Night Band (DNB) onboard the Suomi National Polar-Orbiting Partnership (SNPP) satellite. These types of imagery have been increasingly exploited in studies of urbanization dynamics [9][10][11][12] and estimation of socioeconomic statistics [13][14][15]. In particular, the VIIRS images (2012-present), as an advanced NTL product, are superior to the former coarse-resolution DMSP/OLS (1992-2013) in many respects, for instance, higher spatial resolution (~500 m) and time resolution, without blooming and oversaturated pixels [16][17][18]. Furthermore, equipped with an onboard calibration system, the VIIRS images capture NTL intensity at a higher radiometric resolution than DMSP data [19][20][21][22][23]. Accordingly, there is a growing interest in taking advantage of VIIRS images in urban area extent extraction and socioeconomic estimation [24][25][26][27][28][29]. Due to the capacity of detecting lower light levels, raw VIIRS images contain some temporary light pixels and background noises; many studies have also focused on filtering algorithms or other processes to reduce the impact of noises. In addition to these studies, some researchers have paid attention to the relationship between urban NTL illumination and development disparity [30][31][32]. With the development of light emitting diode (LED) technology, the global transition to LED lights has not only allowed for an increase in light emission but also dramatically changed lighting spectra. Although, global spectral shift likely affects the ability of existing sensors such as DMSP and VIIRS to quantify artificial lights from space [16,17,33]. Given that they are not measuring incoming light in the blue band, the NTL illumination provides an intuitive approach to monitoring human activities.
From the view of nighttime lights, an urban area hotspot can be described as the region of nighttime light hotspot involving multiple bright pixels in an urban area as an ideal proxy of the area with intensified human activities [30,34], which can be regarded spatially as the urban center. NTL illumination mostly derives from artificial lights of urban areas, such as buildings, parking lots, and other elements of the built environment, which represent the level of urbanization development and are closely related to human activity. Therefore, the spatial distribution of NTL illumination intuitively reflects the regional spatial heterogeneity of urbanization [30]. Tracing urban-area hotspots not only reveals the reciprocity between urban expansion and human activities but also provides new insight into sustainable urbanization for government policymakers.
In recent years, some scholars have explored urban area hotspots using coarseresolution DMSP NTL images. Xiao et al. [31] calculated centroids of urban agglomerations from NTL imagery and analyzed the dynamic shift of the centroids. Zhao et al. [32] studied the spatial differentiation and morphologic characteristics of urban core zones in China based on geomorphologic partition theory. Cai et al. [34] identified the structure of polycentric cities using multi-source geospatial big data. In addition, Zheng et al. [30] utilized the Gaussian volume model to depict urban NTL hotspot features quantitatively. Nevertheless, studies regarding identification of urban area hotspots by VIIRS images have not been reported. Given the superiority of satellite sensors, the VIIRS image is more sensitive in detecting low and high NTL intensity [35,36], which offers the possibility to identify urban area hotspots from the perspective of a good data source. However, due to the differences in city size and urban expansion pattern, the spatial distribution of VIIRS NTL intensity is not always suitable for the Gaussian volume model. Therefore, it is necessary to find a universal approach to identify the polycentric urban area hotspots.
The Support Vector Machine (SVM) method proposed by Vapnik [37] is a promising machine learning technique describing the nonlinear relationship between model inputs and outputs with good generalization ability [38]. It has been increasingly used in remote sensing applications as classifier or regression models [39][40][41][42]. The SVM regression model serves as a sample statistical theory which can be used for function fitting for scattered points in the transformation space by a generalized linear function. Thus, it is expected to obtain a continuous prediction output. For instance, Zhan et al. [43] and Yang et al. [44] utilized the SVM regression model to depict urban heat islands (UHI). The UHI is a phenomenon that occurs in urban areas and is associated with significantly higher ambient temperature within a city compared to the surrounding non-urban area. In contrast, an urban area hotspot represents a region that includes a substantial number of lit pixels with higher NTL intensity. We will employ the SVM regression method to investigate the urban area hotspots due to the similarities between the spatial characteristics of UHI and urban area hotspots.
In this study, we designed a practicable framework to explore the space-time dynamics of urban area hotspots by monthly VIIRS NTL data. The primary objectives of our study were as follows: (1) to identify the urban area hotspots by using the SVM regression model, and (2) to analyze the spatiotemporal trajectory of the urban area hotspots. utilized the SVM regression model to depict urban heat islands (UHI). The UHI is a phenomenon that occurs in urban areas and is associated with significantly higher ambient temperature within a city compared to the surrounding non-urban area. In contrast, an urban area hotspot represents a region that includes a substantial number of lit pixels with higher NTL intensity. We will employ the SVM regression method to investigate the urban area hotspots due to the similarities between the spatial characteristics of UHI and urban area hotspots.

Study Area and Datasets
In this study, we designed a practicable framework to explore the space-time dynamics of urban area hotspots by monthly VIIRS NTL data. The primary objectives of our study were as follows: (1) to identify the urban area hotspots by using the SVM regression model, and (2) to analyze the spatiotemporal trajectory of the urban area hotspots.

Study Area
Wuhan (29°58 ′ − 31°22 ′ N,113°41 ′ − 115°05 ′ E), the capital city of Hubei province and the largest city in central China, plays a crucial part in the economy, industry, and transportation of the whole country. Wuhan is located in the eastern part of Jianghan plain and has a total land area of 8569.15 km 2 divided among 13 municipal districts (Jianghan, Jiangxia, Qiaokou, Jiang'an, Huangpi, Xinzhou, Hongshan, Wuchang, Hanyang, Caidian, Hannan, Dongxihu, Qingshan) ( Figure 1). Over recent years, Wuhan has grown to a total population and urban built-up area of 10.22 million and 534.28 km 2 in 2013, and 11.08 million and 723.74 km 2 in 2018, respectively. The rapid urbanization in Wuhan resulted in a striking increase in artificial light illumination from 2013 to 2018. Thus, Wuhan is regarded as a suitable area for our study.

Nighttime Light Data and Preprocessing
Monthly VIIRS NTL images for 70 months from January 2013 to December 2018 were acquired from NOAA/NGDC (http://ngdc.noaa.gov/eog/viirs/download_viirs_ntl.html accessed on 1 June 2020) (see Table 1). Two versions of monthly composite VIIRS products, the "vcm" (any data impacted by stray light were excluded) and "vcmsl" (stray lightimpacted data were corrected) data, were used in this study. The preprocessing of monthly VIIRS images involved the following three steps [26]: (1) since the monthly VIIRS NTL imagery is an initial product containing lit pixels from cities, towns, and other luminous sources, the annual "vcm-orm-ntl" composite of 2015 was utilized to filter temporary light pixels and background noises. In addition, given the unavailability of the "vcm" and "vcmsl" data in June 2018, we used the corresponding data in June 2017 to represent them.
(2) All zero values in the "vcm" data were replaced by "vcmsl" data from the corresponding month to gain more reliable monthly VIIRS data. Subsequently, a Savitzky-Golay (SG) filter was adopted to remove abnormal observations in sequential monthly VIIRS NTL data while retaining their overall trends. It is a smoothing algorithm based on local polynomial least-squares fitting [45], which has been extensively applied in time-series analysis of remote sensing imagery [46,47]. (3) All preprocessed monthly VIIRS NTL images were re-projected to the Asia_North_Albers_Equal_Area_Conic projection and resampled at 480 m. The preprocessed datasets were used to generate annual stable VIIRS NTL data by quantifying the temporal changes of VIIRS NTL brightness based on the time-series analysis method (Section 3.1).

Auxiliary Data
In contrast to traditional statistical survey data, point of interest (POI) data provides a more precise location-based representation of urban spatial characteristics related to human activities. Therefore, POI data were used to explore the correlation between the spatial distribution of hotspots and POIs. In our research, a total of 328,803 POIs from Wuhan in 2018 were obtained using Baidu Map, which consisted of more than 15 types of artificial surface features, such as restaurants, schools, residences, cinema, shopping, and parks ( Figure 2). A kernel density analysis of POIs in the study area can identify the location of the city center, the higher degree of aggregation of POIs indicating more-intensive human activities. It can be used as a valuable reference for urban area hotspots. Considering the effects of urbanization on urban area hotspots, we selected a builtup urban area and socioeconomic statistics (i.e., population and GDP) as indicators to quantitatively explore the interaction between urban hotspots variation and economic statistics. The socioeconomic indicators from 2013 to 2018 were obtained from the Wuhan statistical yearbook from the Wuhan Bureau of Statistics (http://tjj.wuhan.gov.cn/tjfw/tjnj, accessed June 2020). Due to a lag in statistical data, the Wuhan statistical yearbook 2014 reflected the socioeconomic status of Wuhan city in 2013. Because the Wuhan statistical yearbook 2019 was not available, we used the statistics data obtained from Wuhan statistical yearbook 2018 to represent not only 2017 but also 2018.

Methodology
The methods used in this study included the following three steps: data processing of monthly VIIRS images, SVM regression modeling for urban area hotspot identification, and spatiotemporal trajectory analysis of urban area hotspots ( Figure 3).
Above all, given the unavailability of annual VIIRS NTL images each year, time-series analyses were conducted on monthly VIIRS images to obtain stable annual VIIRS images. Then, we divided the stable annual NTL pixel intensities for each year into a training To ensure the consistency of geographic coordinate systems, we further converted the coordinates of POI to Asia_North_Albers_Equal_Area_Conic projected coordinates. Kernel density analysis was conducted to generate a POI density map to further validate our method for identifying urban area hotspots.
Considering the effects of urbanization on urban area hotspots, we selected a builtup urban area and socioeconomic statistics (i.e., population and GDP) as indicators to quantitatively explore the interaction between urban hotspots variation and economic statistics. The socioeconomic indicators from 2013 to 2018 were obtained from the Wuhan statistical yearbook from the Wuhan Bureau of Statistics (http://tjj.wuhan.gov.cn/tjfw/tjnj, accessed on 1 June 2020). Due to a lag in statistical data, the Wuhan statistical yearbook 2014 reflected the socioeconomic status of Wuhan city in 2013. Because the Wuhan statistical yearbook 2019 was not available, we used the statistics data obtained from Wuhan statistical yearbook 2018 to represent not only 2017 but also 2018.

Methodology
The methods used in this study included the following three steps: data processing of monthly VIIRS images, SVM regression modeling for urban area hotspot identification, and spatiotemporal trajectory analysis of urban area hotspots ( Figure 3). set and testing set for training and validating the SVM regression model for identification of urban area hotspots. To ensure the accuracy and robustness of our SVM regression model, grid search and cross-validation were combined to retrieve the optimal parameters. Finally, the space-time trajectory of the centroid of the urban area hotspot was depicted by the differential distance and angle to reveal its spatiotemporal dynamics. The specific details of our proposed procedures are described in the following sections.

Time-Series Analysis of VIIRS NTL Intensity
Previous studies demonstrated the existence of monthly variation of VIIRS NTL intensity [48,49]. The NTL brightness is extremely sensitive to the influence of land surface albedo change, showing a regular fluctuation during a year. To identify the urban NTL hotspots accurately, we needed to obtain the annual stable component of VIIRS data (annual stable VIIRS data) by the time-series analysis method for 12 monthly VIIRS images. The seasonal variation of NTL intensity was quantified using Equation (1) [26]: where represents the time-series observations of NTL intensity of a pixel for month it, refers to annual stable NTL data, and T denotes 12 months of one year.
Considering the possible existing phenomenon of pixel-level NTL intensity seasonal variation, we performed the curve fitting of seasonal variations test for 2013. As illustrated in Figure 4, the test achieved a good performance at the pixel scale, with a high goodness of fit (R 2 ). The R 2 values for the intra-city pixels (Figure 4a-d) in Wuhan in 2013 were above 0.85, indicating a more apparent NTL intensity seasonality in urban areas. Obviously, the time-series analysis model can capture the general trend of monthly VIIRS NTL intensity and remove the influence of seasonal variation effectively from the observed NTL radiance intensity. Using the time-series analysis method, we obtained reliable annual stable NTL data during 2013-2018. The NTL intensity of annual stable NTL data were split into training sets and test sets for further identifying urban area hotspots using the machine learning method of SVM [43,44]. Above all, given the unavailability of annual VIIRS NTL images each year, time-series analyses were conducted on monthly VIIRS images to obtain stable annual VIIRS images. Then, we divided the stable annual NTL pixel intensities for each year into a training set and testing set for training and validating the SVM regression model for identification of urban area hotspots. To ensure the accuracy and robustness of our SVM regression model, grid search and cross-validation were combined to retrieve the optimal parameters. Finally, the space-time trajectory of the centroid of the urban area hotspot was depicted by the differential distance and angle to reveal its spatiotemporal dynamics. The specific details of our proposed procedures are described in the following sections.

Time-Series Analysis of VIIRS NTL Intensity
Previous studies demonstrated the existence of monthly variation of VIIRS NTL intensity [48,49]. The NTL brightness is extremely sensitive to the influence of land surface albedo change, showing a regular fluctuation during a year. To identify the urban NTL hotspots accurately, we needed to obtain the annual stable component of VIIRS data (annual stable VIIRS data) by the time-series analysis method for 12 monthly VIIRS images. The seasonal variation of NTL intensity was quantified using Equation (1) [26]: where NTL it represents the time-series observations of NTL intensity of a pixel for month it, NTL stable refers to annual stable NTL data, and T denotes 12 months of one year. Considering the possible existing phenomenon of pixel-level NTL intensity seasonal variation, we performed the curve fitting of seasonal variations test for 2013. As illustrated in Figure 4, the test achieved a good performance at the pixel scale, with a high goodness of fit (R 2 ). The R 2 values for the intra-city pixels (Figure 4a-d) in Wuhan in 2013 were above 0.85, indicating a more apparent NTL intensity seasonality in urban areas. Obviously, the time-series analysis model can capture the general trend of monthly VIIRS NTL intensity and remove the influence of seasonal variation effectively from the observed NTL radiance intensity. Using the time-series analysis method, we obtained reliable annual stable NTL data during 2013-2018. The NTL intensity of annual stable NTL data were split into training sets and test sets for further identifying urban area hotspots using the machine learning method of SVM [43,44].

SVM Regression Model
As a machine learning method, the support vector machine (SVM) method has achieved good results in nonlinear regression model fitting [50,51]. The Gaussian radial basis function (RBF) kernel has been generally utilized for the optimal SVM solution for its advantage of nonlinear mapping [52]. Therefore, we employed the SVM regression model with RBF kernel function to identify Wuhan's urban area hotspots from 2013 to 2018 in this study.
The basic theory of SVM regression is demonstrated as follows. Given the form of training data (x ,y ), x represents the training vectors and y denotes their target values. Supposing that the linear regression function is f(x) = W • x + b, the SVM method can be transformed into searching for the extreme value with the following equations: where W represents the weight vector, C denotes regularization constant, ξ , ξ * are positive slack variables, and ε is the parameter of the ε-intensive objective function. The Lagrange multiplier method was used to transform the original problem into a dual problem. After solving the dual problem, the fitting function can be simplified as where m and m * are the corresponding support vectors of samples and k(x, x ) represents the kernel function.
Specifically, the SVM regression model has two significant parameters, the penalty parameter C and kernel coefficient γ. The choice of an appropriate model parameter influences the regression performance directly. As an essential tool for parameter optimization, the grid search method can thoroughly search until the acquisition of model-optimized parameters. In addition, the partition of the training data also impacts the evalua-

SVM Regression Model
As a machine learning method, the support vector machine (SVM) method has achieved good results in nonlinear regression model fitting [50,51]. The Gaussian radial basis function (RBF) kernel has been generally utilized for the optimal SVM solution for its advantage of nonlinear mapping [52]. Therefore, we employed the SVM regression model with RBF kernel function to identify Wuhan's urban area hotspots from 2013 to 2018 in this study.
The basic theory of SVM regression is demonstrated as follows. Given the form of training data (x i ,y i ), x i represents the training vectors and y i denotes their target values. Supposing that the linear regression function is f(x) = W T ·x + b, the SVM method can be transformed into searching for the extreme value with the following equations: where W represents the weight vector, C denotes regularization constant, ξ i , ξ i * are positive slack variables, and ε is the parameter of the ε-intensive objective function. The Lagrange multiplier method was used to transform the original problem into a dual problem. After solving the dual problem, the fitting function can be simplified as where m i and m i * are the corresponding support vectors of samples and k(x, x i ) represents the kernel function.
Specifically, the SVM regression model has two significant parameters, the penalty parameter C and kernel coefficient γ. The choice of an appropriate model parameter influences the regression performance directly. As an essential tool for parameter optimization, the grid search method can thoroughly search until the acquisition of model-optimized parameters. In addition, the partition of the training data also impacts the evaluation results of the parameter combination. Therefore, grid search and k-fold cross-validation (K-CV) methods have been combined to reduce the training error caused by sampling randomness and ensure the accuracy of the model [53].
Referring to Zhan et al. [43] and Yang et al., [44], the SVR model with RBF kernel function can effectively build the non-liner relationship between geographic coordinate vector (longitude and latitude) and surface urban heat island intensity. Therefore, the longitude and latitude served as the inputs of SVR to build the NTL intensity prediction model in this study. The initial C and γ ranged from 2 −15 to 2 15 , and the step interval was 1.
To perform a five-fold CV method, the input dataset was divided into five subsets for SVM regression. Consequently, based on annual stable VIIRS data using the SVM regression algorithm, we obtained the urban area hotspots of Wuhan in 2013, 2015, and 2018.

Spatiotemporal Migration Analysis of Urban Area Hotspot
The morphology of urban area hotspots exhibits different shapes and magnitudes, making it difficult to separately analyze cluster dynamics for each urban area hotspot. Therefore, to explore the relationship between urban area hotspots and socioeconomic indicators, this study adopted the center of gravity method to analyze the spatiotemporal migration patterns of urban area hotspots from the perspective of the whole study area. Over the selected years, the shift of the gravity center reveals the variation of NTL illumination intensity with urbanization level, thereby reflecting the spatial heterogeneity of human activity. The position of the gravity center of the urban area hotspot was calculated as follows: where X t and Y t represent the longitude and latitude of the gravity center of an urban area hotspot at time t, respectively; X i and Y i denote the longitude and latitude of the urban area hotspot at the ith pixel scale, respectively; C ti is the NTL radiance intensity of the ith pixel at time t. Additionally, the shift distance (L) and angle (α) of the gravity center were calculated by Equations (7) and (8), respectively. The magnitude of shift distance indicates the urban expansion intensity during different periods, and the shift angle reveals the main directions of urban expansion.
where L t+1 represents the shift distance of the gravity center of the urban area hotspot in the period t-t + 1.
where α t+1 denotes the shift angle of the gravity center of the urban area hotspot from t to t + 1.  (Figure 6b), which is mainly influenced by the seasonally changed albedo. In the winter, snowfall significantly increased and vegetation decreased, leading to an increase of NTL brightness. In the summer, increased vegetation resulted in a decrease in NTL illumination.

Nighttime Light Data Processing
As Figure 7 illustrates, the pixels in urban areas have clear seasonal variations during 2013-2018. In contrast, background (non-urban area) pixels are not correlated with the seasons. This is because the albedo of land surface shows more apparent and regular changes in the urban areas than those in the non-urban areas. We also calculated the Pearson correlation coefficient between annual stable VIIRS data in 2015 and annual "vcmorm-ntl" composite data of 2015. A high correlation (0.98) was obtained, which supports the validity of the annual stable VIIRS data after taking seasonal variations into consideration. seasonality at the city scale ( Figure 6a) and pixel scale (Figure 6b), which is mainly influenced by the seasonally changed albedo. In the winter, snowfall significantly increased and vegetation decreased, leading to an increase of NTL brightness. In the summer, increased vegetation resulted in a decrease in NTL illumination.     changes in the urban areas than those in the non-urban areas. We also calculated the Pearson correlation coefficient between annual stable VIIRS data in 2015 and annual "vcmorm-ntl" composite data of 2015. A high correlation (0.98) was obtained, which supports the validity of the annual stable VIIRS data after taking seasonal variations into consideration.  Figure 8c) were conducted separately. Specifically, we used the average VIIRS NTL intensity in March-May, June-August, September-November, and December-February to represent spring, summer, autumn, and winter, respectively. As Figure 8 demonstrates, a fairly good linear relationship of NTL intensity between winter and the other seasons (spring, summer, autumn) was observed, with a slope close to 1. They all show a high correlation (R 2 > 0.9), with a p-value of <0.01. The results show that seasonal curve fitting can effectively reduce the impact of the seasonal variation of NTL illumination. To verify the performance of seasonal curve fitting, linear regression analyses of pixellevel VIIRS NTL intensity between winter and spring (Figure 8a), summer (Figure 8b), and autumn ( Figure 8c) were conducted separately. Specifically, we used the average VIIRS NTL intensity in March-May, June-August, September-November, and December-February to represent spring, summer, autumn, and winter, respectively. As Figure 8 demonstrates, a fairly good linear relationship of NTL intensity between winter and the other seasons (spring, summer, autumn) was observed, with a slope close to 1. They all show a high correlation (R 2 > 0.9), with a p-value of <0.01. The results show that seasonal curve fitting can effectively reduce the impact of the seasonal variation of NTL illumination. VIIRS NTL intensity in March-May, June-August, September-November, and December-February to represent spring, summer, autumn, and winter, respectively. As Figure 8 demonstrates, a fairly good linear relationship of NTL intensity between winter and the other seasons (spring, summer, autumn) was observed, with a slope close to 1. They all show a high correlation (R 2 > 0.9), with a p-value of <0.01. The results show that seasonal curve fitting can effectively reduce the impact of the seasonal variation of NTL illumination.     Figure 9b,c shows the district-scale relationship between the total population growth and total NTL Inc , and the GDP growth per capita and mean NTL Inc , respectively. We have detected the possible noise of data (such as big outliers) according to the Pauta criterion (3σ criterion) [54]. The result reveals that all data points are within the allowable error range; thus, the obtained linear regression model is acceptable. There is a certain positive correlation between NTL growth, on the one hand, and GDP and population, on the other.  Figure 9b,c shows the district-scale relationship between the total population growth and total NTL , and the GDP growth per capita and mean NTL , respectively. We have detected the possible noise of data (such as big outliers) according to the Pauta criterion (3σ criterion) [54]. The result reveals that all data points are within the allowable error range; thus, the obtained linear regression model is acceptable. There is a certain positive correlation between NTL growth, on the one hand, and GDP and population, on the other.  Figure 10 shows the parameter optimization results of the SVM regression model in 2013. The grid search combined with five-fold CV serves as an effective strategy to optimize the parameters. As Figure 10 shows, the regression accuracy of C and γ in a specific range is very high, and the mean square error (MSE) of the fitted result in a specific range is very low. If the parameters were not selected properly, the accuracy would be too low. Therefore, we further narrowed the search range, and the values of C (1.6818) and γ (16384) with the highest accuracy were regarded as the best parameters for identifying the urban area hotspot in 2013. Subsequently, the optimal parameters for the other years  Figure 10 shows the parameter optimization results of the SVM regression model in 2013. The grid search combined with five-fold CV serves as an effective strategy to optimize the parameters. As Figure 10 shows, the regression accuracy of C and γ in a specific range is very high, and the mean square error (MSE) of the fitted result in a specific range is very low. If the parameters were not selected properly, the accuracy would be too low. Therefore, we further narrowed the search range, and the values of C (1.6818) and γ (16384) with the highest accuracy were regarded as the best parameters for identifying the urban area hotspot in 2013. Subsequently, the optimal parameters for the other years (2015, 2018) were also acquired by following the above steps. After parameter optimization, the SVM regression model had good generalization ability for simulating the NTL intensity. In the study, all SVM regression models achieved good performance with a R 2 ≈ 0.8 in 2013, 2015, and 2018. According to the simulated NTL intensity distribution, we can depict the urban area hotspots involving multiple bright pixels with high NTL intensity. Figure 11  According to the simulated NTL intensity distribution, we can depict the urban area hotspots involving multiple bright pixels with high NTL intensity. Figure 11 shows the urban area hotspots from 2013 to 2018. The urban hotspots in 2013 (Figure 11a) were mainly concentrated in the Qingshan, Wuchang, and Jiang'an municipal districts. In 2015, the urban area hotspots expanded toward the southwest of Wuhan (Figure 11b). In 2018, the distribution of urban area hotspots was more balanced in the central area of Wuhan (Figure 11c), which is consistent with the distribution of POI density in 2018 (Figure 11d).  Table 2 shows the shift distance (L) and angle (α) of the urban area hotspots' gravity center during the study period. The gravity center of urban area hotspots shifted to the southwest (SW) direction by 3.85 km from 2013 to 2015. Furthermore, the gravity center of urban area hotspots kept moving toward the SW direction by 2.76 km during 2015-2018.   Table 2 shows the shift distance (L) and angle (α) of the urban area hotspots' gravity center during the study period. The gravity center of urban area hotspots shifted to the southwest (SW) direction by 3.85 km from 2013 to 2015. Furthermore, the gravity center of urban area hotspots kept moving toward the SW direction by 2.76 km during 2015-2018.  (Figure 12c). These comparisons reveal the consistency between the growth of urban area hotspots and urban expansion, which supports the validity of identification of urban area hotspots in this study. To some extent, the shift distance (L) and angle (α) of the urban area hotspots' gravity center reveal the imbalance between different districts in Wuhan. During 2013-2018, the SW direction regions of Wuhan achieved more rapid development than other regions. On the one hand, it may be related to the prosperity of industrial parks located in southwest Wuhan, such as the Wuhan Economic and Technological Development Zone (EDZ) and the Dongxihu Development Zone (DDZ) [56]. On the other hand, the portion of the population from southwest districts of Wuhan kept increasing gradually during the study period [57]. Densely populated areas with intensified human activities tend to have higher nighttime light radiance, while rural areas tend to be dimly illuminated.

Spatiotemporal Dynamic Analysis of Urban Area Hotspots
To further reveal the interaction between the spatiotemporal dynamics of urban hotspots and socioeconomic development related to human activities, we quantitatively explored the coupling between the movement of urban area hotspots and socioeconomic statistics indicators (i.e., population and GDP) during 2013-2018. Referring to the Wuhan statistical yearbooks, Figure 13 shows the movement of the population and GDP gravity center of Wuhan during 2013-2018. We utilized two metrics, the space-overlaps as S and changes in consistency as θ [58], to quantify the coupling relationship. S measures the spatial distance between the two gravity centers in the same year. The farther the spatial distance, the lower the spatial coupling. θ refers to the cosine of the angle between the vectors of the displacement of the two gravity centers relative to the previous time point. The range of θ is [−1, 1]. When θ is equal to 1, it means that the two gravity centers move in the same direction. While θ is equal to -1, it means that the two gravity centers move in opposite directions. Therefore, we defined S and S to represent the space-overlaps between the moving trajectory of urban area hotspots and population gravity center, GDP, respectively. Meanwhile, θ and θ expressed changes in consistency between the moving trajectory of urban area hotspots and population gravity center, respectively. To some extent, the shift distance (L) and angle (α) of the urban area hotspots' gravity center reveal the imbalance between different districts in Wuhan. During 2013-2018, the SW direction regions of Wuhan achieved more rapid development than other regions. On the one hand, it may be related to the prosperity of industrial parks located in southwest Wuhan, such as the Wuhan Economic and Technological Development Zone (EDZ) and the Dongxihu Development Zone (DDZ) [56]. On the other hand, the portion of the population from southwest districts of Wuhan kept increasing gradually during the study period [57]. Densely populated areas with intensified human activities tend to have higher nighttime light radiance, while rural areas tend to be dimly illuminated.
To further reveal the interaction between the spatiotemporal dynamics of urban hotspots and socioeconomic development related to human activities, we quantitatively explored the coupling between the movement of urban area hotspots and socioeconomic statistics indicators (i.e., population and GDP) during 2013-2018. Referring to the Wuhan statistical yearbooks, Figure 13 shows the movement of the population and GDP gravity center of Wuhan during 2013-2018. We utilized two metrics, the space-overlaps as S and changes in consistency as θ [58], to quantify the coupling relationship. S measures the spatial distance between the two gravity centers in the same year. The farther the spatial distance, the lower the spatial coupling. θ refers to the cosine of the angle between the vectors of the displacement of the two gravity centers relative to the previous time point. The range of θ is [−1, 1]. When θ is equal to 1, it means that the two gravity centers move in the same direction. While θ is equal to -1, it means that the two gravity centers move in opposite directions. Therefore, we defined S u−p and S u−G to represent the space-overlaps between the moving trajectory of urban area hotspots and population gravity center, GDP, respectively. Meanwhile, θ u−p and θ u−G expressed changes in consistency between the moving trajectory of urban area hotspots and population gravity center, respectively.  Table 3 illustrates the S and θ values in 2013, 2015, and 2018. During the study period, the higher θ values revealed that the population had more influence on the increment of NTL intensity. With the development of urbanization, the θ increased faster, indicating the growing effect of GDP on the formation of urban area hotspots. The decreasing S indicated that the effect of densely populated areas on the moving trajectory of urban area hotspots increased significantly. Moreover, the S first decreased and then rose, showing that the effect of GDP on the moving trajectory of urban area hotspot tended to rise first and then fall. In conclusion, the results manifest a good spatial coupling degree between urban area hotspots and socioeconomic development in Wuhan from 2013 to 2018. Table 3. Space-overlaps (S) and changes in consistency (θ) of city-level urban area hotspot and population and GDP gravity center.

Year
Urban

Discussion
NTL luminosity provides a unique opportunity to capture NTL clusters representing urban area hotspots closely related to human activities. Time-series nighttime images can  Table 3 illustrates the S and θ values in 2013, 2015, and 2018. During the study period, the higher θ u−p values revealed that the population had more influence on the increment of NTL intensity. With the development of urbanization, the θ u−G increased faster, indicating the growing effect of GDP on the formation of urban area hotspots. The decreasing S u−p indicated that the effect of densely populated areas on the moving trajectory of urban area hotspots increased significantly. Moreover, the S u−G first decreased and then rose, showing that the effect of GDP on the moving trajectory of urban area hotspot tended to rise first and then fall. In conclusion, the results manifest a good spatial coupling degree between urban area hotspots and socioeconomic development in Wuhan from 2013 to 2018. Table 3. Space-overlaps (S) and changes in consistency (θ) of city-level urban area hotspot and population and GDP gravity center.

Discussion
NTL luminosity provides a unique opportunity to capture NTL clusters representing urban area hotspots closely related to human activities. Time-series nighttime images can be employed to monitor the spatiotemporal trajectory of urban area hotspots in a more direct and consistent fashion, compared to traditional census data. In this study, we proposed an applicable framework to monitor the spatiotemporal variation of urban area hotspots using VIIRS NTL images. Considering the existing phenomenon of NTL illumination seasonal variations, we employed the time-series analysis model to obtain the annual stable VIIRS data, and then the SVM regression method was used to identify urban NTL clusters. After parameter optimization by the grid search and K-CV techniques, nonlinear fitting with the SVM regression model presented a good ability to simulate urban area hotspots.
To validate the identification of urban area hotspots of Wuhan, we delineated the POIs kernel density map to recover the spatial distribution characteristics of urban area hotspots. By comparing the POIs map, the identified urban area hotspot well represented the inner urban spatial structure. Meanwhile, we compared the space-time trajectory of urban area hotspots with movement of the gravity center of urban areas extracted from the annual stable VIIRS data from 2013 to 2018. The results demonstrated the high consistency of direction of movement between the urban sprawl and urban area hotspot. The shift distance and angle of urban area hotspots' centroid further revealed the discrepancy of urban expansion during different periods, which will be beneficial to more sustainable urban planning and management. The coupling between the gravity center of the urban hotspot and socioeconomic statistics was analyzed for population and GDP. The results demonstrate that high GDP and dense population have a positive correlation with the increase of NTL intensity.
However, despite the merits of this study, a few data limitations should be mentioned. Firstly, existing NTL products (DMSP/OLS, VIIRS/DNB, Luojia1-01, etc.) do not include the blue channel, and the increasing emission of blue light from LED sources may result in less total NTL emissions estimation. That is, VIIRS has a panchromatic band, making it incapable of measuring incoming light in the blue band, which may underestimate the total NTL emissions. Meanwhile, due to the unavailability of annual standard VIIRS products each year released by NOAA, we recognize the uncertainty of the annual stable VIIRS images obtained by time-series analysis method. Secondly, given that the population statistics represented the permanent resident population obtained from the statistical yearbook, the uncounted floating population and inaccurate statistics may affect the coupling relationship between the hotspot and the population. With the development of big data technology, it is possible to employ mobile internet data, such as mobile signaling [59], taxi-accessible tracks [60,61], check-in data in Sina web [62], and GPS data [63] as an effective reference for urban area hotspot categories. In a future study, we will explore urban area hotspot categories based on high-resolution NTL imagery combining more multi-source data.
Serving as a useful proxy for human activity, urban area hotspots reflect the spatial disparities among urban clusters due to urban expansion. A timely and explicit understanding of urban area hotspots is important for urban planning and environmental management. Due to the finer spatial and temporal resolution of VIIRS NTL imagery, the spatial distribution of NTL intensity surface may not always follow Gaussian surface exploited to simulate the urban hotspot from the DMSP NTL images [30]. Furthermore, the Gaussian volume model needs to treat each potential urban cluster separately, making the process of identifying urban area hotspots of polycentric cities more tedious. Referring to the UHI identification application, the SVM regression method is innovatively adopted for identification of urban area hotspots in this study. For irregular and complex spatial patterns of NTL intensity clusters distribution, the SVM regression method performs well as a fast, effective tool in simulating urban area hotspots, possessing good learning ability combined with grid search and k-fold CV techniques. Furthermore, we exhibited the diversified morphology of the urban area hotspots in the whole city over time. From the spatiotem-poral analysis of urban area hotspots, we detected the coupling between the urban area hotspots and socioeconomic indicators (i.e., population and GDP). The finding shows the diversity of the urban NTL hotspots related to imbalanced human activities, although these statistical data may have some errors. This study could provide significant insights for decision-making analysis in urban development. Considering the urban expansion discrepancy and complexity of the hotspot migration in different polycentric cities, we will investigate the genericity and applicability of the whole method to other contexts in monitoring urban area hotspot migration in a future study.

Conclusions
This paper presents an innovative framework to identify urban area hotspots using VIIRS NTL images. Urban area hotspots were identified by the optimal SVM regression model combined with the grid search and k-fold cross-validation (K-CV) methods. A case study was used to depict the spatiotemporal dynamics of urban area hotspots within Wuhan during 2013-2018. The results confirm that the optimal SVM regression models successfully simulated urban area hotspots with a regression coefficient over 0.8. The results also prove that the proposed approach effectively depicts the location, spatial morphology, and migration trajectory of urban area hotspots, reflecting the explicit spatial variation of human activity with urban expansion.
Due to the rapid socioeconomic development and intensified human activity associated with urbanization, the spatiotemporal trajectory of urban area hotspots was consistent with the trend of urbanization development. This study explored the coupling between urban area hotspots and socioeconomic indicators (i.e., population and GDP). The results demonstrate that the relationship between the growth of urban hotspots and socioeconomic development manifests a high degree of spatial coupling. The NTL intensity increased in direct response to the increase of GDP and population. These findings prove that our proposed framework provides a novel perspective on the interaction between urban sprawl and human activities.
In future studies, the following issues should be considered: (1) the incorporation of more multi-source data to investigate urban area hotspots (such as cellular signaling data); (2) the possibility of urban area hotspot category identification based on high-resolution NTL imagery (such as Luojia 1-01 data).
Author Contributions: Yanhong Zou, Minghui Chen and Yuling Ruan conceived and designed the study. Yanhong Zou, Yuling Ruan, and Jingya Shen collected and processed the data, performed the analysis, and wrote the paper. Minghui Chen contributed to data analysis and writing. Yanhong Zou and Yuling Ruan edited the draft, and approved the submitted manuscript. All authors have read and agreed to the published version of the manuscript.