Mapping Impervious Surface Areas Using Time-Series Nighttime Light and MODIS Imagery

Mapping impervious surface area (ISA) dynamics at the regional and global scales is an important task that supports the management of the urban environment and urban ecological systems. In this study, we aimed to develop a new method for ISA percentage (ISA%) mapping using Nighttime Light (NTL) and MODIS products. The proposed method consists of three major steps. First, we calculated the Enhanced Vegetation Index (EVI)-adjusted NTL index (EANTLI) and performed intra-annual and inter-annual corrections on the DMSP-OLS data. Second, based on the geographically weighted regression (GWR) model, we built a consistent NTL product from 2000 to 2019 by performing an intercalibration between DMSP-OLS and VIIRS images. Third, we adopted a GA-BP neural network model to monitor ISA% dynamics using NTL imagery, MODIS imagery, and population data. Taking the Guangdong–Hong Kong–Macao Greater Bay as the study area, our results indicate that the ISA% in our study area increased from 7.97% in 2000 to 17.11% in 2019, with a mean absolute error (MAE) of 0.0647, root mean square error (RMSE) of 0.1003, Pearson’s coefficient of 0.9613, and R2 (R-squared) of 0.9239. Specifically, these results demonstrate the effectiveness of the proposed method in mapping ISA and investigating ISA dynamics using temporal features extracted from consistent NTL and MODIS products. The proposed method is feasible when generating ISA% at a large scale at high frequency, given the ease of implementation and the availability of input


Introduction
The world has experienced rapid urbanization in the past decades. As of 2014, more than 54% of the global population resided in urban areas, and by 2050, 66% of the world's population are expected to be living in urban areas [1]. The rapid growth of urban population and economic activities indicate that China is experiencing the largest mass urbanization in human history [2,3]: the percent of the urban population having increased from 11.18% in 1950 to 36.22% in 2000 and to 60.6% in 2019 [4]. Such a rapid urbanization process has caused enormous impacts on urban climate [5][6][7], plant phenology [8,9], habitat dynamics [10], biogeochemical cycles [11], and hydrological patterns [12][13][14]. Therefore, long-term urban expansion monitoring is necessary for sustainable urban development. Impervious surfaces have been recognized as the most direct embodiment of urbanization, a key attribute that can represent urban environmental quality and is one of the predominant indicators of human settlements. Hence, numerous studies have been conducted to characterize impervious surface distributions with an array of remote sensing based methods that vary at the global, continental, national, regional, and city scales [15][16][17][18][19][20]. The authors of [21][22][23] provided overviews of the methods that estimate impervious surfaces using remote sensing imagery.
A better understanding of impervious surface dynamics contributes to the knowledge of the urbanization process and its potential influence on the urban ecological environment. data. Li et al. [52] used the VIIRS and DMSP-OLS archives to build a power function model that generates VIIRS-like DMSP-OLS data. Zheng et al. [53] developed a geographically weighted regression model (GWRc) for NTL data cross-sensor calibration.
In this study, we propose a novel approach to map the dynamics of ISA% from multisource data, including DMSP-OLS and VIIRS NTL composite data, MODIS products, and population data. Specifically, the objectives of this study are summarized below: (1) We adopted the Enhanced Vegetation Index (EVI)-adjusted NTL index (EANTLI) to mitigate saturation and blooming in DMSP-OLS imagery, improve the light intensity variations and spatial details in urban cores, and reduce the overestimation of impervious surfaces in suburban areas. (2) We adopted Seasonal and Trend decomposition using Loess (STL) to decompose monthly VIIRS NTL data and estimate annual VIIRS NTL composites to match the temporal resolution of DMSP-OLS data. (3) We built a cross-sensor calibration model to generate consistent NTL time series data by combining EANTLI and annual VIIRS NTL images. (4) We developed a multisource-driven ISA% estimation model for impervious surface mapping.
The rest of this paper is organized as follows. Section 2 describes the data sources used in this study and presents our approach for mapping ISA%. Section 3 presents the results. Section 4 presents the discussion, and Section 5 summarizes this study.

Materials and Methods
To evaluate the performance of our proposed model, we took Guangdong-Hong Kong-Macao Greater Bay ( Figure 1) as a study case. Located in southern China, the geographic location of Guangdong-Hong Kong-Macao Greater Bay is between 111 • 12 E-115 • 35 E and 21 • 25 N-24 • 30 N, with a total population of 70 million in 2017. The Greater Bay contains the Pearl River Delta, the largest urban agglomeration in the world, and two special administrative regions of Hong Kong and Macao. The study area is one of the most open and economically dynamic regions in China. Given its rapid growth, this area experienced a long period of unbalanced internal development [54], evidenced by its notably heterogeneous spatial urbanization patterns and distribution of impervious surfaces. Figure 2 shows the workflow of the proposed ISA% mapping model. In this study, we used three datasets: NTL, MODIS, and population data. These datasets capture three urban elements: economic activity, ecology, and population. The MODIS EVI data and DMSP-OLS data from 2000 to 2013 were used to calculate EANTLI. Intra-annual and interannual corrections were performed to improve the performance of EANTLI. In addition, we utilized the STL algorithm to decompose the monthly VIIRS NTL images and calculate the annual VIIRS NTL values. We adopted a GWR model for intercalibration of the EANTLI and VIIRS NTL data, and utilized the consistent NTL data, Land Surface Temperature (LST), Improved Normalized Difference Built-up Index (NDBVI), and population data for ISA% mapping.  Figure 2 shows the workflow of the proposed ISA% mapping model. In this study, we used three datasets: NTL, MODIS, and population data. These datasets capture three urban elements: economic activity, ecology, and population. The MODIS EVI data and DMSP-OLS data from 2000 to 2013 were used to calculate EANTLI. Intra-annual and interannual corrections were performed to improve the performance of EANTLI. In addition, we utilized the STL algorithm to decompose the monthly VIIRS NTL images and calculate the annual VIIRS NTL values. We adopted a GWR model for intercalibration of the EANTLI and VIIRS NTL data, and utilized the consistent NTL data, Land Surface Temperature (LST), Improved Normalized Difference Built-up Index (NDBVI), and population data for ISA% mapping.   The NTL data we used included DMSP-OLS and NPP-VIIRS images. The Version-4 DMSP-OLS NTL data set contains cloud-free annual composites collected from 1992 to 2013, covering an area from 180 • E to 180 • W and 65 • S to 75 • N. The 34 annual NTL images over a 22-year period are all in the 30 arc-second grids (approximately 1 km at the equator), with DN values ranging from 0 to 63. The DMSP/OLS NTL products we used are the version 4 "stable lights" series compiled over a 22-year span (1992-2013). The version 4 DMSP/OLS stable lights product already excludes sunlight, glare, moonlight, cloud coverage, and lighting [55,56], as well as ephemeral events such as wildfires. In this study, we retrieved the DMSP-OLS NTL dataset with a temporal coverage from 2000 to 2013.
The NPP-VIIRS NTL data, with mitigated saturation and blooming issues, was produced from the Visible Infrared Imaging Radiometer Suite (VIIRS) Day/Night Band (DNB) [57]. We used the Version 1 monthly composite VIIRS data with a 15 arc-second spatial resolution and a temporal coverage from April 2012 to December 2019. VIIRS DNB and DMSP-OLS Nighttime Light products were retrieved from the Earth Observation Group, Payne Institute for Public Policy, Colorado School of Mines (https://eogdata.mines.edu/ products/dmsp/ (accessed on 8 September 2020); https://eogdata.mines.edu/products/ vnl/ (accessed on 10 September 2020).

MODIS Data
The MODIS data that we used in this study contained Enhanced Vegetation Index (EVI), Normalized Difference Vegetation Index (NDVI), Land Surface Temperature (LST), and the surface spectral reflectance of Terra MODIS Bands 2 and 6. The MOD11A2 products include average LST values under clear-sky during an 8-day period and are stored on a Sinusoidal grid at a spatial resolution of 1 km. We averaged the 8-days period LST values per pixel for the whole year to derive an annual composite LST value. The MOD13A3 product provides monthly EVI and monthly NDVI data with a spatial resolution of 1 km in a Sinusoidal projection. The EVI is an improvement on traditional NDVI, as EVI avoids the saturation problem found in the NDVI for high vegetation coverage areas and considers the influence of the soil background on vegetation index. Moreover, the EVI uses the blue band to remove residual atmosphere contamination caused by smoke and thin clouds. We derived the annual average EVI and annual average NDVI to further mitigate the sensitivity to seasonal and inter-annual fluctuations. We used the annual average EVI and DMSP-OLS to calculate the EANTLI to mitigate saturation and blooming in the DMSP-OLS data. The MOD09A1 product provides an estimation of the surface spectral reflectance corrected for atmospheric conditions in the Terra MODIS Bands 1 through 7. The value of each pixel is selected from all the acquisitions within the 8-day composite period at a spatial resolution of 500 m.
The Normalized Difference Built-up Index (NDBI) is calculated as the reflectance ratio of short-wave infrared band and near-infrared band. Its improved version, the Improved Normalized Difference Built-up Index (NDBVI), outperforms the original NDBI in mapping urban land-use areas [58]. NDBI can be calculated using the following equation: where R SWIR and R NIR are the surface reflectance of the short wave infrared band (Band 6) and near-infrared band (Band 2) derived from the MOD09A1 product, respectively.
where NDBI avg describes the average annual NDBI derived from Equation (1), NDVI avg represents the annual average NDVI derived from MOD13A3 monthly NDVI product. MODIS data were downloaded from https://ladsweb.modaps.eosdis.nasa.gov/ (accessed on 28 September 2020), reprojected to the WGS 1984 geographic coordinate system, and resampled to the same cell size as the NTL data.

Population Data
The population data were downloaded in Geotiff format at a resolution of 30 arcsecond (approximately 1 km at the equator) in the WGS 1984 geographic coordinate system from the WorldPop (https://www.worldpop.org (accessed on 15 October 2020). WorldPop, collaborating with many national statistical offices and ministries of health, was initiated in 2011. These data mainly focus on improving the spatial demographic knowledge for low and middle-income countries. WorldPop develops data integration and disaggregation methods and utilizes census, survey, satellite, and cell phone data to produce consistent gridded outputs to supplement traditional data sources, such as the census. The WorldPop program produces a standardized and spatiotemporally consistent set of gridded products at the global scale, providing an open-access archive of highresolution population distribution data [59][60][61]. We collected the 2000-2019 population data from the "Global mosaics 2000-2020" datasets.

Auxiliary Data
We choose Landsat 5 TM images and Landsat 8 OLI/TIRS images from the United States Geological Survey website (https://earthexplorer.usgs.gov/ (accessed on 10 January 2021) to make a reference dataset of impervious surface coverage. The Landsat 5 TM achieves global coverage once every 16 days at a spatial resolution of 30 m (visible, NIR, MIR bands), 120 m (thermal band). The Landsat 8 satellite includes the Operational Land Imager (OLI) and the Thermal Infrared Sensor (TIRS) and provides seasonal coverage of the global landmass at a spatial resolution of 30 m (visible, NIR, SWIR bands), 100 m (thermal bands), and 15 m (panchromatic band). Several preprocessing steps were involved, including radiation calibration, atmospheric correction, and reprojection to the WGS 1984 geographic coordinate system. A random sampling technique was used to collect samples for evaluating the accuracy of the ISA% estimation in a full dynamic range (ISA% from 0% and 100%). The sampling unit for each reference sample is 0.008333 • × 0.008333 • in the WGS 1984 geographic coordinate system. In order to improve the classification accuracy, we utilized the visual interpretation method in ENVI 5.3 to distinguish all the Landsat pixels for a corresponding location (ISA or non-ISA). The ISA% was calculated for each reference sample.

NTL Data Correction
The strongly negative correlation between impervious surfaces and vegetation has been widely used to address the saturation and blooming effect seen in NTL data [34,48,62]. Light saturation reduces light intensity variations, as well as eliminating the spatial details in urban cores. The blooming effect in suburban areas results in overestimation of the lighted area. Hence, the EANTLI is used to mitigate the saturation problem in urban cores and the blooming effect in suburban areas [48]. Mathematically, the EANTLI can be written as: where NTL norm denotes the normalized DMSP-OLS NTL value, NTL is the original DMSP-OLS NTL value, and EVI is derived from the MOD13A3 annual average. More details can be found in Zhou et al. [48]. After the correction using Equation (3), the inconsistencies in lit pixels still exist, as DMSP-OLS data are composed of multi-year and multi-satellite image series. To take Remote Sens. 2021, 13, 1900 7 of 20 advantage of two images from the same year derived from two satellites, we derived an intra-annual composition of the NTL images: where DN intra (n,i) represents the DN value of pixel i after the intra-annual composition in the year n, DN 1 (n,i) and DN 2 (n,i) represent the DN values of pixel i from two EANTLI images in the year n, and n denotes the number of years from 2000 to 2013.
Further, we derived an inter-annual composition by assuming a continuously increasing trend in the DN values, given global urbanization processes [34,63]: where DN inter (n,i) denotes the DN value of the pixel i after the inter-annual correction in the year n. DN (n−1,i) , DN (n,i) , DN n+1,i are the DN values of pixel i from the intra-annual composited NTL images in the years n − 1, n, and n + 1, respectively.

VIIRS Time Series Decomposition
DMSP-OLS data are an annual composite, while VIIRS is a monthly composite. Therefore, we aggregated VIIRS into an annual composite product to match the temporal resolution of the DMSP-OLS data, as seasonal effects can appear in monthly VIIRS. Cities' nighttime brightness tends to fluctuate with snow coverage and phenology of vegetation. Snow coverage and reduced vegetation cover lead to increased NTL intensity, while increased vegetation cover can shield luminosity, leading to reduced NTL intensity [64].
The annual VIIRS composite was estimated, assisted by the STL time series decomposition algorithm. We removed the seasonal and noisy signals from the monthly VIIRS composite before aggregation. The STL algorithm decomposes the VIIRS time series (Y t ) of each pixel into a trend component (T t ) that represents an upward or downward movement of the time series, a seasonal component (S t ) that describes repetitive pattern over time, and an irregular component (e t ) that denotes the remaining VIIRS variation [65,66]: where t = 1, . . . , N. N (N = 93) denotes the total number of VIIRS observations. STL is an iterative non-parametric regression method involving the subjective selection of smoothing parameters. The seasonal window widths (s.window) should be an odd integer greater or equal to 7. The trend window widths (t.window) is determined as follows: where the frequency is set to 12, as VIIRS time series data are monthly. The stl() uses loess regression to calculate seasonal decomposition of Y t by determining T t , then computes S t (and e t ) from the difference Y t − T t , which sets t.window to the minimum value of Equation (7). The function stl(time series, s.window = "periodic") is referred to as the "standard approach". More information regarding the STL method can be found in Cleveland et al. [67]. The annual composite VIIRS was aggregated for specific years by calculating the mean value of VIIRS intensity of the corresponding trend component, which removes seasonal and other irregular components.

GWR-Based Intercalibration Model
The GWR model was used to perform intercalibration between EANTLI and VIIRS data to create consistent long-term NTL data. Given that the DMSP-OLS data fall short in spatial, temporal, and radiometric resolutions, VIIRS requires a degrading procedure before the intercalibration model can be established. In the previous step, the STL algorithm already aggregated the temporal resolution of VIIRS from monthly to annual, aiming to match the temporal resolution of DMSP-OLS. The VIIRS data were resampled to 1 km (VIIRS 1km ). Each VIIRS 1km pixel consisted of roughly four VIIRS pixels. A 5 × 5 window of Gaussian filter was applied to VIIRS in order to reduce the radiometric resolution. Parameter σ was set as 1.75 according to Li et al. [52]. A GWR model combines spectral and spatial information via the following equation: represent constant term and regression coefficient of the i-th pixel, respectively; ε i denotes the residual of the model. A Gaussian function was selected as the kernel function, and a corrected Akaike Information Criterion (AICc) was adopted to determine the best bandwidth setting. The GWR model was applied to build a regression model between EANTLI est and VIIRS 1km in 2012 and 2013, respectively. The results showed that the residuals of the two regression models were similar, and the difference between the residuals was not statistically significant (p > 0.05). Hence, we conclude that the regressive relationship of the study area was stable over time. We used VIIRS 1km and GWR model in 2013 to generate EANTLI-like NTL data. We obtained a consistent annual NTL time series from 2000 to 2019, formed by EANTLI in 2000-2013 and EANTLI-like VIIRS data in 2014-2019.

ISA% Estimation Model
To detect ISA% and calculate accuracy, we utilized an optimizing back-propagation (BP) neural network algorithm based on a genetic algorithm (GA), i.e., GA-BP neural network model. GA is an evolutionary algorithm based on population-based evolution. GA deploys an objective function to evaluate the fitness of each chromosome in the population. After a series of operations that include selection, crossover, and mutation, GA eliminates chromosomes with lower fitness to obtain a new population. These operations are repeated until the chromosome reaches a certain standard. The structure of the BP neural network consists of input layers, hidden layers, and output layers. It sets the initial weights randomly in the forward direction of the training process. The output data can be obtained according to certain rules. The weights are modified in the backward process based on the difference between the output data and ground-truth values. A repetitive procedure is implemented of the forward and backward processes until the difference between the output data and the ground-truth values is satisfactory. The training of the GA-BP neural network model starts with the GA. The GA performs a global search for the optimal weights, and the BP algorithm starts the training process with the best initial weights to approach the optimum solution. Therefore, the combination of GA and BP avoids local optimums and reduces the learning time [68,69]. We used MAE, RMSE, Pearson's coefficient, and R 2 as indicators to evaluate the accuracy.

Spatiotemporal Dynamics of ISA%
In this study, we adopted the global and local Moran's I indices to describe the spatiotemporal dynamics of ISA% from 2000 to 2019 across the study area. Ranging from −1 to 1, the global Moran's I index is an overall measure of spatial autocorrelation. A positive value indicates the level of similarity, while a negative value denotes the degree of difference. The global Moran's I index can be written as: where N represents the number of cities, ω ij denotes the matrix of the spatial weight of ith and jth cities based on the Queen's contiguity, x i and x j denote the ISA% for the ith and jth cities, and x is the average ISA% of the entire study area. The local Moran's I index denotes the degree of spatial autocorrelation between each sample and its neighbors: where The local Moran's I index describes four types of local spatial autocorrelation: (1) a High-High cluster (a high value is surrounded by high values); (2) a Low-Low cluster (a low value is surrounded by low values); (3) a High-Low cluster (a high value is surrounded by low values); (4) a Low-High cluster (a low value is surrounded by high values).
The SLOPE index can be applied to examine the variation trends of cities' ISA% from 2000 to 2019: where n (n = 20) denotes the number of years, y t is the tth year from 1 to 20. FIS i_t denotes the ISA% in the ith city of tth year. A positive SLOPE i represents an increasing trend, while a negative value denotes otherwise.

Evaluation of the NTL Data Correction
The saturation effect exists in NTL imagery for cities at different socioeconomic and development levels. Figures 3 and 4 demonstrate the saturation correction effects, intraannual correction, and inter-annual correction on DMSP-OLS data for the Guangdong-Hong Kong-Macao Greater Bay region. As displayed in Figure 3a, DMSP-OLS images show little variation in the DN value in certain areas, especially around the cities of Guangzhou, Dongguan, Shenzhen, Hong Kong, Foshan, and Zhongshan. Meanwhile, as shown in Figure 3b, the extended value range in the EANTLI image leads to improved DN value distinguishability. Figure 4 presents four major cities where the DN values from DMSP-OLS were capped at 63. The corresponding EANTLI images, however, display mitigated saturation results, evidenced by the heterogeneity.
Hong Kong-Macao Greater Bay region. As displayed in Figure 3a, DMSP-OLS images show little variation in the DN value in certain areas, especially around the cities of Guangzhou, Dongguan, Shenzhen, Hong Kong, Foshan, and Zhongshan. Meanwhile, as shown in Figure 3b, the extended value range in the EANTLI image leads to improved DN value distinguishability. Figure 4 presents four major cities where the DN values from DMSP-OLS were capped at 63. The corresponding EANTLI images, however, display mitigated saturation results, evidenced by the heterogeneity.    Figure 5 demonstrates the results of VIIRS time series decomposition with STL, given different changing patterns of land coverage. Seasonal, trend, and irregular components are successfully separated. Bars to the right side of each plot in Figure 5 have the same length but have relatively different sizes as they are plotted on different scales. The long bars in the seasonal plots indicate that the time series decompose the prominence seasonality signal, which illustrates the necessity to strip away the seasonal components to obtain the trend components. The seasonal and irregular signals can be clearly observed, as the sum of seasonal and irregular components accounted for 23.70% and 27.50% of the original VIIRS time series in Figures 5a,b, respectively. These irregular components (i.e., the abnormal increase in values in June 2018 and February 2019) affected monthly VIIRS time series, presumably leading to calculation errors. STL is able to detect trend components and seasonal components, thus removing the impacts from these outliers in the original time series. The above results suggest that removing seasonal and irregular components from monthly VIIRS data via STL is a necessary step to obtain consistent yearly VIIRS products.  Figure 5 demonstrates the results of VIIRS time series decomposition with STL, given different changing patterns of land coverage. Seasonal, trend, and irregular components are successfully separated. Bars to the right side of each plot in Figure 5 have the same length but have relatively different sizes as they are plotted on different scales. The long bars in the seasonal plots indicate that the time series decompose the prominence seasonality signal, which illustrates the necessity to strip away the seasonal components to obtain the trend components. The seasonal and irregular signals can be clearly observed, as the sum of seasonal and irregular components accounted for 23.70% and 27.50% of the original VIIRS time series in Figure 5a,b, respectively. These irregular components (i.e., the abnormal increase in values in June 2018 and February 2019) affected monthly VIIRS time series, presumably leading to calculation errors. STL is able to detect trend components and seasonal components, thus removing the impacts from these outliers in the original time series. The above results suggest that removing seasonal and irregular components from monthly VIIRS data via STL is a necessary step to obtain consistent yearly VIIRS products.

Evaluation of the Intercalibration Model
Based on the GWR model, we constructed an intercalibration model linking EANTLI and VIIRS to obtain EANTLI-like VIIRS data. The Guangdong-Hong Kong-Macao Greater Bay region has experienced rapid development and urbanization in the past decades, which led to significant differences in the socioeconomic development of the cities in this area. Thus, we constructed GWR models according to the city size and its economic development. The study area was divided into seven zones, including Foshan, Huizhou, Jiangmen, Zhaoqing, Zhongshan_Zhuhai_Macao, Shenzhen_Hong Kong, and Guangzhou_Dongguan. Table 1 shows the performance of GWR models for seven zones, as well as for the entire study area. Indicators to evaluate the performance consist of MAE, RMSE, Pearson's coefficient, and R 2 . Using the intercalibration models, we derived the EANTLI-like VIIRS data in 2014-2019. Figure 6 shows the DN value distribution between the EANTLI-like VIIRS data and the EANTLI data along the latitudinal transect. The results demonstrate that EANTLI-like VIIRS data present similar DN values and are capable of expressing considerable heterogeneity, suggesting that the intercalibration models we constructed are effective.

Evaluation of the Intercalibration Model
Based on the GWR model, we constructed an intercalibration model linking EANTLI and VIIRS to obtain EANTLI-like VIIRS data. The Guangdong-Hong Kong-Macao Greater Bay region has experienced rapid development and urbanization in the past decades, which led to significant differences in the socioeconomic development of the cities in this area. Thus, we constructed GWR models according to the city size and its economic development. The study area was divided into seven zones, including Foshan, Huizhou, Jiangmen, Zhaoqing, Zhongshan_Zhuhai_Macao, Shenzhen_Hong Kong, and Guangzhou_Dongguan. Table 1 shows the performance of GWR models for seven zones, as well as for the entire study area. Indicators to evaluate the performance consist of MAE, RMSE, Pearson's coefficient, and R 2 . Using the intercalibration models, we derived the EANTLI-like VIIRS data in 2014-2019. Figure 6 shows the DN value distribution between the EANTLI-like VIIRS data and the EANTLI data along the latitudinal transect. The results demonstrate that EANTLI-like VIIRS data present similar DN values and are capable of expressing considerable heterogeneity, suggesting that the intercalibration models we constructed are effective.

ISA% Mapping
In this study, the GA-BP model was applied to map annual ISA% from 2000 to 2019. For BP neural network, NTL, LST, NDBVI, and population data served as input feature datasets. The reference dataset of impervious surface coverage derived from Landsat images was used as the output. A three-layer neural network was applied with four input layers, six hidden layers, and one output layer. In the genetic algorithm, the population size, crossover rate, and mutation rate were set as 50, 0.2, and 0.1, respectively. The yearly ISA% distribution over the study period was divided into eleven grades, as shown in Figure 7.

ISA% Mapping
In this study, the GA-BP model was applied to map annual ISA% from 2000 to 2019. For BP neural network, NTL, LST, NDBVI, and population data served as input feature datasets. The reference dataset of impervious surface coverage derived from Landsat images was used as the output. A three-layer neural network was applied with four input layers, six hidden layers, and one output layer. In the genetic algorithm, the population size, crossover rate, and mutation rate were set as 50, 0.2, and 0.1, respectively. The yearly ISA% distribution over the study period was divided into eleven grades, as shown in Figure 7. Cities with a lower level of economic development, such as Huizhou, Zhaoqing, Jiangmen, Zhuhai, and Zhongshan, experienced increased ISA% especially in urban cores, and the average annual growth rate of ISA% in these cities have increased dramatically. In Guangzhou, Shenzhen, Foshan, and Dongguan, zones with a particularly high speed of urbanization, impervious surfaces expanded from urban cores to suburban areas, so the average annual growth rate of ISA% maintained relatively stable. In comparison, Hong Kong and Macau are in a mature period of impervious surface expansion, as they already reached a high level of urbanization before 2000. Impervious surfaces in Hong Kong and Macau demonstrated a slight expansion during the investigated period, with a relatively slow average annual growth rate of ISA% and unchanged ISA spatial patterns. Remote Sens. 2021, 13, x FOR PEER REVIEW 14 of 21 Cities with a lower level of economic development, such as Huizhou, Zhaoqing, Jiangmen, Zhuhai, and Zhongshan, experienced increased ISA% especially in urban cores, and the average annual growth rate of ISA% in these cities have increased dramatically. In Guangzhou, Shenzhen, Foshan, and Dongguan, zones with a particularly high speed of urbanization, impervious surfaces expanded from urban cores to suburban areas, so the average annual growth rate of ISA% maintained relatively stable. In comparison, Hong Kong and Macau are in a mature period of impervious surface expansion, as they already reached a high level of urbanization before 2000. Impervious surfaces in Hong Kong and Macau demonstrated a slight expansion during the investigated period, with a relatively slow average annual growth rate of ISA% and unchanged ISA spatial patterns.

Accuracy Assessment
We compared our estimated ISA% maps with the GAIA 1985-2018 and FROM-GLC10 products developed by Gong et al. [28,70]. In their researches, Gong et al. [28] utilized Landsat images and ancillary datasets that consist of NTL data and the Sentinel-1 Synthetic Aperture Radar data to develop an automatic mapping procedure on the Google Earth Engine (GEE) platform to provide annual global artificial impervious areas (GAIA) from 1985 to 2018 at a 30 m resolution. Additionally, Gong et al. [70]  To quantitatively assess the performance of our model in estimating ISA% maps, 500 samples from historical Landsat images were interpreted to calculate MAE, RMSE, Pearson's coefficient, and R 2 between the ground-truth values and predicted values derived from different products. Table 2 displays the accuracy assessment results. As shown in Table 2

Accuracy Assessment
We compared our estimated ISA% maps with the GAIA 1985-2018 and FROM-GLC10 products developed by Gong et al. [28,70]. In their researches, Gong et al. [28] utilized Landsat images and ancillary datasets that consist of NTL data and the Sentinel-1 Synthetic Aperture Radar data to develop an automatic mapping procedure on the Google Earth Engine (GEE) platform to provide annual global artificial impervious areas (GAIA) from 1985 to 2018 at a 30 m resolution. Additionally, Gong et al. [70] developed FROM-GLC10 by transferring a training sample set of 2015 at 30 m resolution to classify 10 m resolution images of 2017.
To quantitatively assess the performance of our model in estimating ISA% maps, 500 samples from historical Landsat images were interpreted to calculate MAE, RMSE, Pearson's coefficient, and R 2 between the ground-truth values and predicted values derived from different products. Table 2 displays the accuracy assessment results.

Spatiotemporal Analysis of ISA% during 2000-2019
As Figure    As shown in Figure 8, the degree of spatial autocorrelation of ISA% in Guangzhou and Foshan was significantly higher than in the other cities. As the capital of Guangdong Province, Guangzhou demonstrated a relatively stable trend, while the degree of spatial autocorrelation of ISA% in Foshan increased from 0.7717 in 2000 to 0.8693 in 2019. Huizhou and Zhaoqing witnessed a higher growth rate compared to cities. The Global Moran's I of Huizhou increased from 0.4441 to 0.6561, whereas the Global Moran's I of Zhaoqing increased from 0.3769 to 0.5574 over these 20 years. Jiangmen also showed an upward trend in the Global Moran's I value, while the value in Zhongshan changed from 0.4928 to 0.5959. The degree of spatial autocorrelation in Hong Kong demonstrated a relatively steady trend. The degree of spatial autocorrelation in Dongguan jumped in the first eight years of the study period and then showed a slight downward tendency. Similar patterns can be seen in Shenzhen. Compared to all other cities, Zhuhai stands out, given its continuously declining degree of spatial autocorrelation of ISA% during the study period. Figure 9 shows Local Moran's I indices of the Greater Bay in 2000 and 2019, with detected spatial cluster modes. NS indicates insignificant cluster mode. HH represents the spatial clustering of high ISA% values. HL suggests that a grid with a high ISA% is surrounded by ones with a low ISA%. LH suggests that a grid with low ISA% is surrounded by ones with high ISA%. LL indicates the spatial clustering of low ISA% values. We can observe the HH clusters are distributed in Guangzhou, Foshan, Zhongshan, Dongguan, Shenzhen, and Hong Kong, while the LL clusters exist in Zhaoqing, Jiangmen, Huizhou, and the north of Guangzhou.
As shown in Figure 8, the degree of spatial autocorrelation of ISA% in Guangzhou and Foshan was significantly higher than in the other cities. As the capital of Guangdong Province, Guangzhou demonstrated a relatively stable trend, while the degree of spatial autocorrelation of ISA% in Foshan increased from 0.7717 in 2000 to 0.8693 in 2019. Huizhou and Zhaoqing witnessed a higher growth rate compared to cities. The Global Moran's I of Huizhou increased from 0.4441 to 0.6561, whereas the Global Moran's I of Zhaoqing increased from 0.3769 to 0.5574 over these 20 years. Jiangmen also showed an upward trend in the Global Moran's I value, while the value in Zhongshan changed from 0.4928 to 0.5959. The degree of spatial autocorrelation in Hong Kong demonstrated a relatively steady trend. The degree of spatial autocorrelation in Dongguan jumped in the first eight years of the study period and then showed a slight downward tendency. Similar patterns can be seen in Shenzhen. Compared to all other cities, Zhuhai stands out, given its continuously declining degree of spatial autocorrelation of ISA% during the study period. Figure 9 shows Local Moran's I indices of the Greater Bay in 2000 and 2019, with detected spatial cluster modes. NS indicates insignificant cluster mode. HH represents the spatial clustering of high ISA% values. HL suggests that a grid with a high ISA% is surrounded by ones with a low ISA%. LH suggests that a grid with low ISA% is surrounded by ones with high ISA%. LL indicates the spatial clustering of low ISA% values. We can observe the HH clusters are distributed in Guangzhou, Foshan, Zhongshan, Dongguan, Shenzhen, and Hong Kong, while the LL clusters exist in Zhaoqing, Jiangmen, Huizhou, and the north of Guangzhou. Distinctive changes in the HH cluster appear in ISA% for Guangzhou, Foshan, Jiangmen, Zhongshan, and Huizhou, suggest a considerably increased ISA%. In addition, we calculated the Local Moran's I for Huizhou, Guangzhou, and Zhuhai. We selected these three cities because of their representativeness, given their different urbanization and levels of spatial autocorrelation in the ISA%. The changes in their spatial clustering of ISA% over time are displayed in Figure 10  Distinctive changes in the HH cluster appear in ISA% for Guangzhou, Foshan, Jiangmen, Zhongshan, and Huizhou, suggest a considerably increased ISA%. In addition, we calculated the Local Moran's I for Huizhou, Guangzhou, and Zhuhai. We selected these three cities because of their representativeness, given their different urbanization and levels of spatial autocorrelation in the ISA%. The changes in their spatial clustering of ISA% over time are displayed in Figure 10. We observed a notable change in HH clusters from 2000 to 2019 in Huizhou. The expansion of HH clusters in Guangzhou demonstrated an uptrend from 2000 to 2010 but maintained a stable level from 2010 to 2019. In comparison, changes in Zhuhai were not as clear as in the other two cities. Figure 11 shows the temporal variation of the ISA% in different cities from 2000 to 2019. Rapid growth occurred in Guangzhou and Foshan, corresponding to their high rate of urbanization during these 20 years. As the three cities with the lowest per capita gross domestic product in the Greater Bay in 2018, Zhaoqing and Jiangmen experienced moderate growth, and Huizhou demonstrated relatively stable growth in ISA%. Remote Sens. 2021, 13, x FOR PEER REVIEW 17 of 21   of urbanization during these 20 years. As the three cities with the lowest per capita gross domestic product in the Greater Bay in 2018, Zhaoqing and Jiangmen experienced moderate growth, and Huizhou demonstrated relatively stable growth in ISA%.

Conclusions
The complexity of impervious surface features impedes mapping at the regional and global scale. NTL images, as an emerging data source, have been shown to be a reliable data source for mapping ISA. Nevertheless, we need to properly deal with saturation and blooming problems, as well as the inconsistency in DMSP-OLS data. In addition, calibrations are needed before a long-term, consistent NTL time series can be formed to ensure the compatibility and comparability of images from DMSP-OLS and VIIRS satellite sensors. Furthermore, ISA extraction based on single-source remote sensing data has widely recognized limitations that include low adaptability and low universality. Therefore, we made the following improvements. (1) We used EANTLI to address the saturation and blooming issues in DMSP-OLS data. Intra-annual and inter-annual corrections were performed to address the inconsistencies in lighted pixels. (2) We generated long-term and consistent NTL data from 2000 to 2019 by combining EANTLI and VIIRS. The trend components of monthly VIIRS data were extracted and used to build a cross-sensor calibration model based on the GWR model. (3) We utilized multisource temporal features extracted from consistent NTL data, MODIS, and population data to establish an ISA% estimation model.
Taking the Guangdong-Hong Kong-Macao Greater Bay region as a study area, the proposed method was utilized to map ISA% from 2000 to 2019. The results indicate that ISA% in the study area increased from 7.97% in 2000 to 17.11% in 2019. An accuracy assessment shows that our proposed method achieves great performance in ISA% mapping, with MAE of 0.0647, RMSE of 0.1003, Pearson's coefficient of 0.9613, and of 0.9239. A spatiotemporal analysis of ISA% reveals significant spatial autocorrelation in the entire Greater Bay, as well as in each city during the investigated 20-year period. Thus, we infer that cities experienced divergent patterns in ISA% change over time.

Conclusions
The complexity of impervious surface features impedes mapping at the regional and global scale. NTL images, as an emerging data source, have been shown to be a reliable data source for mapping ISA. Nevertheless, we need to properly deal with saturation and blooming problems, as well as the inconsistency in DMSP-OLS data. In addition, calibrations are needed before a long-term, consistent NTL time series can be formed to ensure the compatibility and comparability of images from DMSP-OLS and VIIRS satellite sensors. Furthermore, ISA extraction based on single-source remote sensing data has widely recognized limitations that include low adaptability and low universality. Therefore, we made the following improvements. (1) We used EANTLI to address the saturation and blooming issues in DMSP-OLS data. Intra-annual and inter-annual corrections were performed to address the inconsistencies in lighted pixels. (2) We generated long-term and consistent NTL data from 2000 to 2019 by combining EANTLI and VIIRS. The trend components of monthly VIIRS data were extracted and used to build a cross-sensor calibration model based on the GWR model. (3) We utilized multisource temporal features extracted from consistent NTL data, MODIS, and population data to establish an ISA% estimation model.
Taking the Guangdong-Hong Kong-Macao Greater Bay region as a study area, the proposed method was utilized to map ISA% from 2000 to 2019. The results indicate that ISA% in the study area increased from 7.97% in 2000 to 17.11% in 2019. An accuracy assessment shows that our proposed method achieves great performance in ISA% mapping, with MAE of 0.0647, RMSE of 0.1003, Pearson's coefficient of 0.9613, and R 2 of 0.9239. A spatiotemporal analysis of ISA% reveals significant spatial autocorrelation in the entire Greater Bay, as well as in each city during the investigated 20-year period. Thus, we infer that cities experienced divergent patterns in ISA% change over time.
Our method can be used to generate high-frequency ISA% maps at a large scale due to its ease of implementation and the availability of data sources. These high-frequency ISA% maps can be used to analyze the socioeconomic, environmental dynamics, and urban processes at a regional or global scale. As more NTL satellite sensors are starting to emerge