An Improved Correction Method of Nighttime Light Data Based on EVI and WorldPop Data

Defense Meteorological Satellite Program’s Operational Linescan System (DMSP/OLS) data has the shortcomings of discontinuous and pixel saturation effect. It was also incompatible with the Soumi National Polar-Orbiting Partnership Visible Infrared Imaging Radiometer Suite (NPP/VIIRS) data. In view those shortcomings, this research put forward the WorldPop and the enhanced vegetation index (EVI) adjusted nighttime light (WEANTL) using EVI and WorldPop data to achieve intercalibration and saturation correction of DMSP/OLS data. A long time series of nighttime light images of china from 2001 to 2018 was constructed by fitting the DMSP/OLS data and NPP/VIIRS data. Corrected nighttime light images were examined to discuss the estimation ability of gross domestic product (GDP) and electric power consumption (EPC) on national and provincial scales, respectively. The results indicated that, (1) after correction, the nighttime light (NTL) data can guarantee the growth trend on national and regional scales, and the interannual volatility of the corrected NTL data is lower than that of the uncorrected NTL data; (2) on the national scale, compared with the established model of NTL data and GDP data (NTL-GDP), the determination coefficient (R2) and the mean absolute relative error (MARE) are 0.981 and 8.518%. The R2 and MARE of the established model of NTL data and EPC data (NTL-EPC) were 0.990 and 4.655%; (3) on the provincial scale, the R2 and MARE of NTL-GDP model under the provincial units are 0.7386 and 38.599%. The R2 value and MARE of NTL-EPC model are 0.8927 and 29.319%; (4) on the provincial scale, the R2 and MARE of NTL-GDP model on time series are 0.9667 and 10.877%. The R2 and MARE of NTL-GDP model on time series are 0.9720 and 6.435%; the established TNL-GDP and TNL-EPC models with 30 provinces data all passed the F-test at the 0.001 level; (5) the prediction accuracy of GDP and EPC on time series was nearly 100%. Therefore, the correction method provided in this research can be applied in estimating the GDP and EPC on multiple scales reliably and accurately.


Introduction
Night light remote sensing, as an active branch of remote sensing [1], can obtain the information of visible light and near infrared electromagnetic wave emitted from the ground surface on cloudless nights [2] It can reflect the urban night lighting and capture the night lighting of fishing boats [3], natural gas combustion [4], and forest fire light [5]. It intuitively reflects the human activities and social development rules [6] and is widely applied in the estimation of socioeconomic parameters such as GDP [7], population density [8], electric power consumption (EPC) [9], greenhouse gas emissions [10] and poverty index [11]; urban expansion [12] and urbanization research [13]; major events assessment such as energy crisis [14], earthquake [15] and war [16]; environmental effect analysis of urban expansion [17]; light pollution [18] and effect analysis [19]. Consequently, it has become the main data source for monitoring human socioeconomic activities and natural phenomena.
The Defense Meteorological Satellite Program's Operational Linescan System (DMSP/OLS) and Soumi National Polar-Orbiting Partnership Visible Infrared Imaging Radiometer Suite (NPP/VIIRS) are the most widely used to obtain the nighttime light data. In this research, it provides nighttime light (NTL) data from 1992 to 2013 and 2012 to 2018, which are named as DMSP/OLS and NPP/VIIRS, respectively.
The Operational Linescan System is one of the principal sensors on the DMSP satellite platform. The radiation value detected by OLS sensors at night is four orders of magnitude lower than that of other sensors [20]. Due to the superior performance of OLS sensors, the NTL data is often used for exploring human activities. However, DMSP/OLS data also has a series of shortcomings [21]. Light saturation in the city center, blooming effect, and lack of airborne calibration are the major obstacles to the application of DMSP/OLS data. The light saturation phenomenon is caused by the low 6-bit radiation resolution of OLS sensor. The pixel value of DMSP/OLS data ranges from 0 to 63. When the radiation brightness of a light source in the city center exceeds to the detection upper limit (10 −8 W/ (cm 2 *sr*um)), the DN values of DMSP/OLS data are all forced to be 63 [22]. Therefore, it is impossible to distinguish the spatial difference of actual lighting intensity in the urban central area where the lighting intensity exceeds the detection limit [23]. Furthermore, the correlation between the detected NTL and economic activity is reduced, which increases errors in the established model [24,25], and limits the application of NTL data [26]. At the same time, DMSP/OLS dataset is made up of NTL data collected by different satellite sensors in different years. In the absence of airborne calibration [4], the night stable light data was affected by the difference of sensors, the difference of satellite cross time, and the influence of the sensor performance degradation [27]. Therefore, the data cannot be directly compared [28]. In addition, the illumination area is overestimated due to the reflection of urban light to urban edge and internal water body, and then increases the error of distribution estimation [29].
In order to improve the accuracy of nighttime light data in representing social and economic activities, many studies have attempted to solve the problems in DMSP/OLS data. Aiming at the saturation problem of DMSP/OLS data, scholars have proposed a series of methods to correct the saturation. The correction method is mainly divided into only using the NTL data and using other satellite data to correct the NTL data. There are three calibration methods that only use the NTL data: (1) radiation calibration desaturation method. This method is used to adjust the sensor gain to make it lower than the typical work settings and obtain the nighttime light change in the city center [30]. The limited nighttime light data obtained can be combined with the nighttime light data obtained under the high gain setting to generate the unsaturated nighttime light data [23]. Although this method can effectively reduce the saturation effect of NTL data, it is labor intensive and cost intensive. Therefore, calibration data only can be obtained with a limited number of years (1996, 1999, 2000, 2002, 2004, 2006, and 2010). (2) Desaturation method based on frequency distribution of the DN value of light pixels [31]. This method is based on the assumption that the variation trend of the DN value of night light in saturated regions which is consistent with that in unsaturated regions, the number of pixels in saturated regions can be predicted by establishing a linear correction model of the number of pixels in unsaturated regions and the DN value [32]. Letu et al. used the cubic regression model to construct the DN value frequency of light pixels to correct the saturation problem of light images on the basis of the linear correction model [33]. This method is a desaturation method on the regional scale which is not suitable for the pixel scale. Okinawa as invariant target areas [34]. This correction method can reduce the saturation of the pixel to a certain extent, but the saturated pixel is corrected to the same DN value through the regression model, which means that the corrected pixels lack spatial difference and cannot reflect the real nighttime light data. At the same time, saturation correction will lead to distortion of unsaturated pixels in suburban and rural areas of underdeveloped areas [35]. There are also three calibration methods that use other satellite data to correct the NTL data: (1) correction of nighttime light data based on vegetation index. Vegetation adjusted normalized urban index (VANUI) [36,37] was proposed by combining normalized difference vegetation index (NDVI) with DMSP/OLS data based on the principle that urban characteristics (including nighttime light) should be inversely proportional to vegetation coverage [38]. Enhanced vegetation index adjusted nighttime light index (EANTLI) [39,40] was proposed combined with the advantages of enhanced vegetation index (EVI) [41] to make it impossible for VANUI to effectively highlight the difference of nighttime light intensity in saturated areas under the rapid urbanization process in China and other countries. EVI data is the vegetation index obtained by the calculation of infrared band, near-infrared band, and blue light band [42] in remote sensing data. It is an important numerical expression method for monitoring the growth information of surface vegetation, which can effectively describe the growth and biomass information of green plants. However, at present, this method is mostly used to correct DMSP/OLS images of specific years, and rarely to correct long time series DMSP/OLS datasets. Moreover, rare studies have been conducted on the consistency correction of DMSP/OLS data based on EVI correction and NPP/VIIRS data. (2) NTL data was corrected based on surface temperature and vegetation index spatial distribution difference of vegetation was a natural factor and the influence of human factors on the evolution of urban internal structure. LST and EVI Regulated NTL City Index (LERNCI) based on EVI, land surface temperature (LST), and nighttime light data was proposed [43]. LERNCI based on the double correction of LST and EVI can better improve the spatial difference caused by single index correction. (3) Other data to correct the NTL data: GDP data can effectively reflect the urban development vitality. Therefore, a linear model between the GDP grid data of corresponding time and the NTL data of unsaturated areas is constructed, and the NTL data of saturated areas can be calculated by the model to eliminate the saturation phenomenon [44].
In view of the lack of airborne calibration, the invariant region method was used to realize intercalibration of DMSP/OLS data [45]. Assuming that there are relatively unchangeable pixels in the multi time nighttime light image, and building calibration model of invariant pixels in different images, the calibration model is used to correct the corresponding night light image to get comparable night light images. Elvidge et al. [4] assumed that the NTL data of Sicily in Italy remained basically unchanged from 1994 to 2008. Based on the baseline image in F121999, the regression model of other images from 1994 to 2008 was constructed to achieve intercalibration. Liu et al. [46] used the second order regression model to realize intercalibration of Chinese DMSP/OLS data from 1992 to 2008 by taking Jixi city of Heilongjiang province as the invariant target area. The above two intercalibration methods solve the problem of incomparability of DMSP/OLS data, but do not effectively correct the saturation problem in the image. How to effectively solve the problem of data saturation in DMSP/OLS data while solving the problem of data incomparability in DMSP/OLS data is still a hot research topic.
With the attenuation and failure of DMSP/OLS, NPP/VIIRS data became the new generation of NTL data [47]. As NPP/VIIRS data with higher spatial and radiating quality than DMSP/OLS data [48], it addresses the shortcomings existing in DMSP/OLS data (saturation of city center light data and lack of airborne calibration), but the published NPP/VIIRS monthly composites is a preliminary product which has not been filtered to screen out lights from aurora, fires, boats, and other temporal lights. To solve this problem, Li et al. [26] proposed to extract the corresponding pixel in NPP/VIIRS data as the stable pixel in NPP/VIIRS data by using the bright pixel in DMSP/OLS data of 2010. Wu et al. [49] used the bright pixel of an annual composite of 2015 to extract the corresponding pixel in the NPP/VIIRS data from 2015 to 2017 as the stable pixel in the NPP/VIIRS data of respective years. Although the above methods can remove some unstable pixels in the NPP/VIIRS data, the first method to extract the are cloud-free coverage, average visible light, and stable light. Among the three types of data, the stable light data contains light from cities, towns, and other locations with continuous lighting, while fires, volcanoes, background noise, and other transient events are discarded [27]. The 20 images were captured by four different DMSP satellites, F14 (2001-2003), F15 (2001-2007), F16 (2004-2009), and F18 (2010-2013). In this paper, the radiometric calibration product of F162010 (https://eogdata.mines.edu/dmsp/download_radcal.html) was used to compare the influence of calibration method on the saturation effect of nighttime lighting data.
NPP/VIIRS data from 2012 to 2018 was also collected from the website of Paynes Institute for Public Policy, Colorado School of Mines (https://eogdata.mines.edu/download_dnb_composites.html). The images used for this study were monthly VIIRS DNB composite from 2012 to 2018 and two annual composites of 2015 and 2016. However, the NPP/VIIRS monthly composite was not filtered to eliminate light detection associated with gas flares, fires, volcanoes or auroras, and the data set was not processed to eliminate background noise. Temporary lighting and background noise were eliminated in the annual composite image.
EVI data is the EVI monthly composite of the moderate resolution imaging spectroradiometer (MODIS). The data set was provided by Geospatial Data Cloud site, Computer Network Information Center, Chinese Academy of Sciences. (http://www.gscloud.cn).
The WorldPop dataset was developed by the WorldPop Project (https://www.worldpop.org). The dataset provides annual gridded population data from the 2000-2020 period, with a spatial resolution of 100 m. The input variables of WorldPop include the most recent official census population data and a wide range of spatial ancillary datasets. The spatial datasets include settlement locations and extents, nighttime satellite images, land cover, roads, building maps, health facility locations, vegetation, topography, and refugee camps. A random forest regression tree-based mapping approach was used to generate a predictive weighting layer to reallocate population counts into gridded pixels [55]. The WorldPop dataset has two products, one is the number of people per hectare, the other is the number of people per grid, and the latter dataset was used in this paper.
Gas burning zones are generated due to insufficient infrastructure in the oil producing areas to make full use of natural gas. The data of the gas burning area used in this paper can be downloaded from the website of Paynes Institute for Public Policy, Colorado School of Mines (https://eogdata.mines.edu/download_global_flare.html).
Vector maps of China's county-level administrative boundaries and rivers and lakes data were obtained from 1:400 million database in the National Fundamental Geographic Information System. All nighttime lighting data were extracted according to China's administrative boundary.
The socioeconomic statistics from 2002 to 2019 were provided by the China statistical yearbook which was published by China National Bureau of statistics (http://www.stats.gov.cn/tjsj/ndsj/). These include two attributes for reference: gross domestic product (GDP) and electric power consumption (EPC). The unit of the GDP is billion RMB, and the unit of EPC is billion·kilowatt·h −1 .
Agriculture is an important part of social and economic activities, especially in a traditional agricultural country like China. However, most of the agricultural production is located in the darker area of the nighttime light image, so the nighttime light data cannot represent the agricultural part of economic activities [56]. Therefore, it is necessary to remove the contribution of agricultural sector to economic activities when using nighttime data to predict GDP. The GDP data in the following paragraphs are all subtracting the data of the agricultural sector.

Methodology
This section introduces the method of NTL data processing and result verification. The method is divided into three steps ( Figure 1): (1) data preprocessing: preprocessing of DMSP/OLS data, NPP/VIIRS data, and EVI data. (2) Calibration method: the calibration method is divided into DMSP/OLS data calibration method and NPP/VIIRS data calibration method. The DMSP/OLS data and NPP/VIIRS data of the same year after calibration are fitted and analyzed by linear regression model to obtain a long time series NTL data in China. (3) Result verification: the accuracy of the results was quantitatively verified by comparing the corrected pixels and analyzing the linear relationship between GDP, EPC, and NTL data.  [4]. To avoid the influence of image grid deformation on night image data and facilitate the calculation of light pixel area, all night light intensity images were first extracted in the administrative region of China. Then, the extracted image coordinate system was transformed to Lambert equal area projection coordinate system and the nearest-neighbor resampling algorithm (NEAREST) was used to resample the mesh in the image to 1000 m.  [4]. To avoid the influence of image grid deformation on night image data and facilitate the calculation of light pixel area, all night light intensity images were first extracted in the administrative region of China. Then, the extracted image coordinate system was transformed to Lambert equal area projection coordinate system and the nearestneighbor resampling algorithm (NEAREST) was used to resample the mesh in the image to 1000 m.

Removal of Water and Gas Burning Areas
The reliability of NTL data can be further improved by removing NTL data from water and natural gas combustion areas in DMSP/OLS and NPP/VIIRS images [57].

Preprocessing of NPP/VIIRS Data
Radiation value of night light image data in January 2013 is from −0.72 to 1171.82 ( Figure 2), but the image is a preliminary data set which is not filtered to eliminate light detection associated with auroras, fires, ships, and other temporary lights [22]. Therefore, it is necessary to reduce the impact of unstable night lighting and background noise on nighttime lighting images. Usually, an appropriate threshold interval is set to remove abnormal values (high and low values) to extract effective nighttime lighting images. In general, threshold interval can be set by percentage method, total light area increase method, and economic intensity method. In this paper, the economic intensity method was used to extract the effective night light image [58]. Beijing, Shanghai, and Guangzhou are the three most developed cities in China. Therefore, the pixel values of all regions in China should theoretically be less than the maximum values of these three cities. Hence, the maximum radiation value of these three cities was taken as the upper limit value of the threshold interval. For example, in January 2013, the maximum radiation value of Beijing, Shanghai, and Guangzhou was 255.68, so 255.68 was used as threshold value to correct abnormal value of China nighttime light image data in January 2013. The pixel points higher than the threshold value in the nighttime light image were

Removal of Water and Gas Burning Areas
The reliability of NTL data can be further improved by removing NTL data from water and natural gas combustion areas in DMSP/OLS and NPP/VIIRS images [57].

Preprocessing of NPP/VIIRS Data
Radiation value of night light image data in January 2013 is from −0.72 to 1171.82 ( Figure 2), but the image is a preliminary data set which is not filtered to eliminate light detection associated with auroras, fires, ships, and other temporary lights [22]. Therefore, it is necessary to reduce the impact of unstable night lighting and background noise on nighttime lighting images. Usually, an appropriate threshold interval is set to remove abnormal values (high and low values) to extract effective nighttime lighting images. In general, threshold interval can be set by percentage method, total light area increase method, and economic intensity method. In this paper, the economic intensity method was used to extract the effective night light image [58]. Beijing, Shanghai, and Guangzhou are the three most developed cities in China. Therefore, the pixel values of all regions in China should theoretically be less than the maximum values of these three cities. Hence, the maximum radiation value of these three cities was taken as the upper limit value of the threshold interval. For example, in January 2013, the maximum radiation value of Beijing, Shanghai, and Guangzhou was 255.68, so 255.68 was used as threshold value to correct abnormal value of China nighttime light image data in January 2013. The pixel points higher than the threshold value in the nighttime light image were identified as the abnormal pixel points. As the abnormal pixel points may gather, in order to ensure the complete elimination of the abnormal pixel points, the average value of 24 pixels around the outliers was assigned to the outliers with average filtering method [59]. After processing, the corrected NTL data in January 2013 were all lower than 255. 68 identified as the abnormal pixel points. As the abnormal pixel points may gather, in order to ensure the complete elimination of the abnormal pixel points, the average value of 24 pixels around the outliers was assigned to the outliers with average filtering method [59]. After processing, the corrected NTL data in January 2013 were all lower than 255.68. Table A1 in the appendix shows the upper limit of nighttime thresholds for each month from 2012 to 2018. In general, NPP/VIIRS monthly composite image has negative radiation value pixels. This is because in the dark current correction process, the average night ocean data in the northern hemisphere is used and there is system noise which makes the pixels in the dark area appear as negative values [60]. Therefore, the pixels with negative radiation value are assigned as zero [61]. After the above two processing steps, the corrected nighttime light image of China in January 2013 was obtained as shown in Figure 3. Compared with Figures 2 and 3, it can be found that after the above processing, the abnormally highlighted pixels and the pixels with radiation value less than zero were processed. Therefore, 60 NPP/VIIRS monthly corrected nighttime light images from 2012 to 2018 were obtained by the above correction methods, and then formula (1) was used to process the NPP/VIIRS monthly corrected nighttime light images to obtain the preliminary corrected NPP/VIIRS annual composite images from 2012 to 2018.  In general, NPP/VIIRS monthly composite image has negative radiation value pixels. This is because in the dark current correction process, the average night ocean data in the northern hemisphere is used and there is system noise which makes the pixels in the dark area appear as negative values [60]. Therefore, the pixels with negative radiation value are assigned as zero [61]. After the above two processing steps, the corrected nighttime light image of China in January 2013 was obtained as shown in Figure 3. Compared with Figures 2 and 3, it can be found that after the above processing, the abnormally highlighted pixels and the pixels with radiation value less than zero were processed. Therefore, 60 NPP/VIIRS monthly corrected nighttime light images from 2012 to 2018 were obtained by the above correction methods, and then Formula (1) was used to process the NPP/VIIRS monthly corrected nighttime light images to obtain the preliminary corrected NPP/VIIRS annual composite images from 2012 to 2018.
where, DN y is the annual average nighttime light value, and DN i is the monthly average nighttime light value in ith month. Due to the serious lack of NPP/VIIRS data from May to July, they are not involved in the calculation, hence: i 5, 6, 7; n = 9.

Preprocessing of EVI Data
In order to reduce the sensitivity of seasonal and grade fluctuations [37], formula (2) was used in this paper to process monthly products to obtain annual average EVI data instead of annual maximum EVI data. Due to the lack of monthly EVI data in 2004 and 2008, the average EVI data in 2003 and 2005 was assigned to the average EVI data in 2004, and the average EVI data in 2007 and 2009 was assigned to the average EVI data in 2008. The 13 EVI images calculated from 2001 to 2013 were extracted from the administrative region of China by ArcGIS, and then the extracted image coordinate system was transformed to Lambert equal area projection coordinate system and the nearest-neighbor resampling algorithm (NEAREST) was used to resample the mesh in the image to 1000 m.
where, y EVI is the annual average EVI value, and i EVI is the monthly average EVI value in ith month.

DMSP/OLS Data Calibration Method
Since DMSP/OLS data was not radiometric calibrated, the existing DMSP/OLS nighttime light image data set was obtained by six different sensors. The image data of different years cannot be directly compared due to the differences in the performance of different sensors, the decline of their detection performance with the use time, and the different imaging environment. Moreover, due to the upper limit of sensor detection, pixel saturation exists in the city center. Based on the principle that urban characteristics (including nighttime light) should be inversely proportional to vegetation coverage [38] and directly proportional to population distribution data [53,54], WorldPop and the enhanced vegetation index adjusted nighttime light (WEANTL) constructed by EVI, WorldPop, and

Preprocessing of EVI Data
In order to reduce the sensitivity of seasonal and grade fluctuations [37], Formula (2) was used in this paper to process monthly products to obtain annual average EVI data instead of annual maximum EVI data. Due to the lack of monthly EVI data in 2004 and 2008, the average EVI data in 2003 and 2005 was assigned to the average EVI data in 2004, and the average EVI data in 2007 and 2009 was assigned to the average EVI data in 2008. The 13 EVI images calculated from 2001 to 2013 were extracted from the administrative region of China by ArcGIS, and then the extracted image coordinate system was transformed to Lambert equal area projection coordinate system and the nearest-neighbor resampling algorithm (NEAREST) was used to resample the mesh in the image to 1000 m.
where, EVI y is the annual average EVI value, and EVI i is the monthly average EVI value in ith month.

DMSP/OLS Data Calibration Method
Since DMSP/OLS data was not radiometric calibrated, the existing DMSP/OLS nighttime light image data set was obtained by six different sensors. The image data of different years cannot be directly compared due to the differences in the performance of different sensors, the decline of their detection performance with the use time, and the different imaging environment. Moreover, due to the upper limit of sensor detection, pixel saturation exists in the city center. Based on the principle that urban characteristics (including nighttime light) should be inversely proportional to vegetation coverage [38] and directly proportional to population distribution data [53,54], WorldPop and the enhanced vegetation index adjusted nighttime light (WEANTL) constructed by EVI, WorldPop, and DMSP/OLS data can restore the light intensity in the saturated area to a certain extent, alleviate the saturation phenomenon of light intensity in the urban center, increase the spatial difference of light intensity, and make the nighttime light data better reflect social and economic activities. Therefore, in this paper, the invariant region method was used to conduct intercalibration of DMSP/OLS data, and then EVI and WorldPop data were used to conduct saturation calibration of DMSP/OLS data, so as to obtain higher quality and analyzable long time series nighttime light images.

Intercalibration
Obtaining historical stable images in China, it is necessary to consider the image data collected by the same sensor in different years and different sensors in the same year. The stable bright pixels obtained by different sensors in the same year should be consistent. The bright pixel of the previous year must be the stable bright pixel of the next year [62].
Processing of nighttime light images in China from 2001 to 2013: in ArcGIS software, the DN value of pixels which is larger than zero is assigned as one and the DN of pixels which is equal to zero remains unchanged by using the grid calculator tool. By assuming that all the pixels is which DN values are equal to one in the processed nighttime light image of F182013 are all stable light pixels, after image intersection processing between the nighttime light image of F182013 and nighttime light image of F182012, the pixels in which DN values are equal to one in the nighttime light image is the stable bright pixel of 2012. The processed nighttime light image is called the stable nighttime light image mask. If there are multiple nighttime light images in the same year, such as the nighttime light image of F162007 and F152007, the stable nighttime light image mask of 2008 and the nighttime light image of F162007 and F152007 processed by binary method are intersecting, processed to obtain the stable nighttime light image mask of 2007. By that analogy, the stable nighttime light image masks from 2001 to 2013 are obtained.
Obtaining the stable night light image of China: each stable nighttime light image mask is used to remove the unstable pixels from the nighttime light image of the corresponding year, and get the stable nighttime light image of China after removing the unstable pixels.
Sample data area and reference calibration nighttime light image selected: according to the principle of small interannual fluctuation of DN values and wide dynamic range distribution in the year [46], Jixi city, Heilongjiang province [63] was selected as the sample area. F162007 nighttime light image with high total gray value and good continuity was used as the reference calibration nighttime light image [62].
Correction model determined: the single quadratic regression model shown in Formula (3) is used to fit and regress the image data of Jixi City in each year with the reference nighttime light image (F162007). The regression parameters of each image are shown in Table A2 in the Appendix A. The regression model of the corresponding year is used to regress and correct the obtained DN value of regional stable nighttime light images in China.
where, DN correct is the DN value after intercalibration; DN is the original DN value; a, b, and c are the regression coefficients of quadratic regression equations.

Saturation Calibration
DMSP/OLS nighttime light images of different years after intercalibration are comparable, but intercalibration does not diminish the saturation effect of pixels in the urban center. Therefore, based on the principle of negative correlation trend between vegetation and human activities [38,[64][65][66] and reducing the impact of soil background and atmosphere on vegetation index and overcoming the NDVI saturation problem [66], and on the basis of the nighttime light data proportional to the population distribution data [53,54], WEANTL was put forward based on correction of EVI, WorldPop, and NTL data.
WEANTL is defined as: where, WorldPop and EVI stand for the corresponding values of pixels to be corrected, respectively.WorldPopAVG S and EVIAVG S are average values of all the saturated NTL pixels that represent the mean averages of WorldPop S and EVI S of the saturated pixels. DN correct is the DN value after intercalibration. Firstly, the DMSP/OLS image was divided into saturated and unsaturated regions [67]. In this paper, the region of pixel point with DN value equal to DN correctmax was defined as the saturated region, and the region of pixel point with DN value less than DN correctmax was defined as the unsaturated region. In this paper, only the pixel points in the saturated area were corrected for saturation calibration, and the pixel points in the unsaturated area remained as the original value.
Therefore, the DN value after saturation calibration is noted as DN WEANTL , then DN WEANTL can be calculated by Formula (5): The WEANTL is based on the assumption that the NTL data are positively related to population data [53,54] and negatively related to EVI data [39,40]. Here, WEANTL is proposed based on the common observation that, as a result of high population densities and decreased vegetation cover, the nighttime light is brighter. We expect that WEANTL can better improve the problem of spatial difference caused by single index correction, distinguish the spatial difference of nighttime light data in saturated areas more effectively, and retain the original night light data in unsaturated areas.

Intra-Annual Composition
In order to make full use of the data information of nighttime light images obtained by two different sensors in the same year and solve the problem of inconsistency of luminous pixels [68], intra-annual composition of nighttime light images after intercalibration and saturation calibration was conducted by Formula (6). Years requiring integration include 2001-2007. Since the unstable pixels in the nighttime light image were eliminated in the intercalibration process, it is only necessary to take the average value of the two nighttime light images as the fusion result of the year.
where, n is equal to 2001, 2002 . . . 2007; DN a (n,i) and DN b (n,i) represent the DN value of ith pixel obtained by two different sensors in nth years nighttime light image after intercalibration and saturation calibration. DN (n,i) is the DN value of ith pixel in nth years nighttime light images after intra-annual composition.

Interannual Series Correction
There is no reverse urbanization process in China [69,70]. Therefore, it is assumed that the bright pixels detected in the nighttime light image of the previous year will not disappear in the next year nighttime light image [71]. Therefore, Formula (7) is used to carry out interannual series correction for nighttime light images in China to ensure that the DN value of pixels in the image of the previous year is not greater than that in the next year image.
where, n is equal to 2001, 2002 . . . 2012; DN (n,i) and DN (n+1,i) represent the DN value of ith pixel in nth years and in (n+1)th year nighttime light image, respectively. After correction, the DN value of each pixel of the nighttime light image in 2013 did not change, but the DN value of the pixel after correction of the nighttime light image before 2013 may be significantly reduced, resulting in the DN value of the early image being less than the DN value after intra-annual composition.
Therefore, Formula (8) [72] was used to carry out the interannual series correction in the opposite direction for the nighttime light image in China, so as to ensure that the DN value of the pixel in the image of the next year is not less than the DN value of the pixel in the image of the previous year.
where, n is equal to 2002, 2003 . . . 2013; DN (n−1,i) and DN (n,i) represent the DN value of ith pixel in (n−1)th years and in nth year nighttime light image, respectively. This adjustment sequence can avoid excessively weakening the DN value of the pixels in the early annual image, but it will also cause the DN value of the pixels in the later annual image to be excessively enlarged. Therefore, the average DN value of the pixels in the same year after the two times interannual series correction is processed by Formula (9) [73], so as to obtain the nighttime light image closer to the actual situation.
where, n is equal to 2001, 2002 . . . 2013; DN c (n,i) represents the DN value of ith pixel in nth year nighttime light image corrected by Formula (7); DN d (n,i) represents the DN value of ith pixel in nth year nighttime light image corrected by Formula (8).

NPP/VIIRS Data Calibration Method
In the data preprocessing process, the abnormal values and negative values in NPP/VIIRS nighttime light image have been processed by threshold method, and monthly composite have been synthesized into annual composite. However, the short-term data such as forest fire, aurora, or volcanic eruption still exist in the annual composite after processing, so the mask method [63] was used to process the annual data synthesized in 2012-2018.
The time series from 2012 to 2018 was divided into two stages: from 2012 to 2015 and from 2016 to 2018. It is hypothesized that there will not be a large number of pixels changing from dim pixel to bright pixel in the two stages. Therefore, in this paper, the annual composite of 2015 NPP/VIIRS annual nighttime light images processed by binary-value method was multiplied by the annual composites from 2012 to 2015, so as to obtain stable annual composites from 2012 to 2015. Then the annual composite of 2016 NPP/VIIRS annual nighttime light images processed by binary-value method was multiplied by the annual composites from 2016 to 2018, so as to obtain stable annual composites from 2016 to 2018.

Calibration Method of Two NTL Data
First of all, the total nighttime lights on the county scale of 31 provinces by using Formula (10) for the corrected DMSP/OLS data and NPP/VIIRS data in 2012 and 2013 was calculated. The single quadratic regression model [52] shown in Formula (11) was chosen as the regression model to carry out the regression analysis on the total nighttime lights of DMSP/OLS data and NPP/VIIRS data in 2012 and 2013, respectively.
where, DN i is the DN value of ith pixel in the nighttime light image, n is the total number of pixels, and TNL is the total nighttime lights.
where, TNL DMSP/OLS represents the total nighttime lights of DMSP/OLS data in 2012 and 2013, TNL NPP/VIIRS represents the total nighttime lights of NPP/VIIRS data in 2012 and 2013, a, b, and c are the regression coefficients of single quadratic regression equations. The regression parameters of each province are shown in Table A3 in the Appendix A. The established regression model was used to fit the NPP/VIIRS data from 2014 to 2018 to the DMSP/OLS data scale. Finally, the nighttime light data set from 2001 to 2018 was calculated.

Calibration Results Verification
In this paper, the feasibility and accuracy of the above correction method were verified by means of quantitative test, which includes the quantitative verification of the pixel correction results and socioeconomic indicators and correction results. In this paper, three indicators were selected for quantitative verification of pixel correction results: total nighttime light (TNL), normalized difference index (NDI), and the sum of normalized difference index (SNDI) [54]. As shown in Formula (13), TNL refers to the sum of night light radiation values in the target area, which is a key indicator to measure the overall trend of regional economic development. NDI (Formula (12)) is an indicator used to evaluate the difference of TNL in two periods of nighttime light images in adjacent years, the greater the absolute value, the greater the difference of light brightness between years, and the greater the likelihood of abnormal fluctuations in the data. SNDI (Formula (13)) is the cumulative sum of the normalized difference indexes in all years, which mainly reflect the total variation trend of light brightness in a certain time range.
where, TNL 1 and TNL 2 are the total nighttime lights of the adjacent two years, respectively, p is the number of years.
In this paper, GDP [74] and EPC [75] were chosen to conduct a quantitative test on the correction results, since economic development patterns vary with natural environment and geographical location [12], the linear model (Equation (14)), quadratic polynomial model (Equation (15)), and power function model (Equation (16)) were used to uncover the optimal socioeconomic indicator dynamics models using nighttime data among China's different regions over time and then the best fitting equation was chosen to carry out fitting analysis, compare the fitting accuracy after the correction, and verify the accuracy of the correction algorithm.
where, G is the socioeconomic indicator data of the target region (GDP or EPC), TNL is the total nighttime light of the corresponding year, a, b, and c are coefficients. The regression model was used to predict the socioeconomic data (GDP or EPC) of the target region, and to evaluate the prediction ability of TNL to the socioeconomic data. For each region, the best fitting equation was used to predict the socioeconomic data, and the relative error (RE) [26] and the mean absolute relative error (MARE) [76] were used to evaluate ability of TNL to predict socioeconomic indicators.
where, G is the socioeconomic indicator data of the target region (GDP or EPC), G is the socioeconomic data calculated by the regression equation, n is the number of years or administrative regions.

Analysis of Pixel Comparison Results
In this paper, 30 provinces in China were divided into four economic regions. Pixel results were compared from economic regions and national scales. In order to illustrate the accuracy of the correction method, night light data before correction and that after correction by the correction method (which was abbreviated as reference correction method below) [63,77] were also compared. The comparison results are shown in Figures 4-6. In Figure 4, the abscissa represents different years, the ordinate represents the TNL, and the curve represents the TNL in different years under different treatments. The distribution results of NDI in different regions before and under the two correction methods are compared in Figure 5. SNDI in four economic regions and national scales is compared in Figure 6.
Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 35 where, G is the socioeconomic indicator data of the target region (GDP or EPC), TNL is the total nighttime light of the corresponding year, a, b, and c are coefficients. The regression model was used to predict the socioeconomic data (GDP or EPC) of the target region, and to evaluate the prediction ability of TNL to the socioeconomic data. For each region, the best fitting equation was used to predict the socioeconomic data, and the relative error (RE) [26] and the mean absolute relative error (MARE) [76] were used to evaluate ability of TNL to predict socioeconomic indicators.
where, G is the socioeconomic indicator data of the target region (GDP or EPC), ' G is the socioeconomic data calculated by the regression equation, n is the number of years or administrative regions.

Analysis of Pixel Comparison Results
In this paper, 30 provinces in China were divided into four economic regions. Pixel results were compared from economic regions and national scales. In order to illustrate the accuracy of the correction method, night light data before correction and that after correction by the correction method (which was abbreviated as reference correction method below) [63,77] were also compared. The comparison results are shown in Figures 4-6. In Figure 4, the abscissa represents different years, the ordinate represents the TNL, and the curve represents the TNL in different years under different treatments. The distribution results of NDI in different regions before and under the two correction methods are compared in Figure 5. SNDI in four economic regions and national scales is compared in Figure 6.   , the total nighttime lighting shows a downward trend on the national scale, which is inconsistent with the urbanization process in China. On the regional scale, it can be seen that from 2008 to 2009, the northeast region shows the opposite change of the total nighttime lighting with that of the national scale, while the remaining regions show a downward trend of the total nighttime lighting in these years. Moreover, from 2009 to 2010, the total nighttime lighting shows an increase trend, which was also inconsistent with the actual change of nighttime lights. The 2014 NPP/VIIRS data is much smaller in magnitude than the 2013 DMSP/OLS data, which results in the data drop from 2013 to 2014. Although NPP/VIIRS data is obtained from the same satellite, the total nighttime lighting also drops at national and regional scales from 2015 to 2016, and the total nighttime lighting drops at western regions from 2017 to 2018. Therefore, the uncorrected nighttime light data do not reflect the socioeconomic situation in China very well.    Figure 4c, it is the change of total nighttime lighting at the national and regional scale after the reference correction method. It can be seen on the national scale from 2001 to 2018, the total nighttime lighting shows the trend of growth, in accordance with urbanization process in China. However, on the regional scale, from 2013 to 2014, the total nighttime lighting in the eastern and northeast regions decreased.
NDI of the nighttime light data (including the DMSP/OLS and NPP/VIIRS data) from 2002 to 2018 were corrected under the correction method in this paper and the reference correction method has improved significantly and the interannual volatility of nighttime lighting data has been greatly  , the total nighttime lighting shows a downward trend on the national scale, which is inconsistent with the urbanization process in China. On the regional scale, it can be seen that from 2008 to 2009, the northeast region shows the opposite change of the total nighttime lighting with that of the national scale, while the remaining regions show a downward trend of the total nighttime lighting in these years. Moreover, from 2009 to 2010, the total nighttime lighting shows an increase trend, which was also inconsistent with the actual change of nighttime lights. The 2014 NPP/VIIRS data is much smaller in magnitude than the 2013 DMSP/OLS data, which results in the data drop from 2013 to 2014. Although NPP/VIIRS data is obtained from the same satellite, the total nighttime lighting also drops at national and regional scales from 2015 to 2016, and the total nighttime lighting drops at western regions from 2017 to 2018. Therefore, the uncorrected nighttime light data do not reflect the socioeconomic situation in China very well.  , the total nighttime lighting shows a downward trend on the national scale, which is inconsistent with the urbanization process in China. On the regional scale, it can be seen that from 2008 to 2009, the northeast region shows the opposite change of the total nighttime lighting with that of the national scale, while the remaining regions show a downward trend of the total nighttime lighting in these years. Moreover, from 2009 to 2010, the total nighttime lighting shows an increase trend, which was also inconsistent with the actual change of nighttime lights. The 2014 NPP/VIIRS data is much smaller in magnitude than the 2013 DMSP/OLS data, which results in the data drop from 2013 to 2014. Although NPP/VIIRS data is obtained from the same satellite, the total nighttime lighting also drops at national and regional scales from 2015 to 2016, and the total nighttime lighting drops at western regions from 2017 to 2018. Therefore, the uncorrected nighttime light data do not reflect the socioeconomic situation in China very well.  , the total nighttime lighting shows a downward trend on the national scale, which is inconsistent with the urbanization process in China. On the regional scale, it can be seen that from 2008 to 2009, the northeast region shows the opposite change of the total nighttime lighting with that of the national scale, while the remaining regions show a downward trend of the total nighttime lighting in these years. Moreover, from 2009 to 2010, the total nighttime lighting shows an increase trend, which was also inconsistent with the actual change of nighttime lights. The 2014 NPP/VIIRS data is much smaller in magnitude than the 2013 DMSP/OLS data, which results in the data drop from 2013 to 2014. Although NPP/VIIRS data is obtained from the same satellite, the total nighttime lighting also drops at national and regional scales from 2015 to 2016, and the total nighttime lighting drops at western regions from 2017 to 2018. Therefore, the uncorrected nighttime light data do not reflect the socioeconomic situation in China very well.  As shown in Figure 4b, it can be seen on the national scale from 2001 to 2018, the total nighttime lighting after calibration shows the trend of growth, in accordance with urbanization process in China. However, on the regional scale, from 2013 to 2014 and from 2017 to 2018, the total nighttime lighting in the northeast region decreased. The total nighttime lighting in the central region decreased from 2013 to 2015. Additionally, the total nighttime lighting in the western region decreased from 2014 to 2015. As shown in Figure 4c, it is the change of total nighttime lighting at the national and regional scale after the reference correction method. It can be seen on the national scale from 2001 to 2018, the total nighttime lighting shows the trend of growth, in accordance with urbanization process in China. However, on the regional scale, from 2013 to 2014, the total nighttime lighting in the eastern and northeast regions decreased.
NDI of the nighttime light data (including the DMSP/OLS and NPP/VIIRS data) from 2002 to 2018 were corrected under the correction method in this paper and the reference correction method has improved significantly and the interannual volatility of nighttime lighting data has been greatly reduced ( Figure 5). Compared with Figure 5b,c, it can be found that the interannual volatility under the correction method in this paper from 2013 to 2014 is lower than that of reference correction method. Therefore, the difference of nighttime light data within each province should be considered in the regression model of DMSP/OLS data and NPP/VIIRS data. The results of regression fitting for each province are better than those on the national scale.
It can also be seen that no matter the correction method in this paper and reference correction method can all improve the interannual fluctuation of nighttime light data ( Figure 6). Before correction, SNDI in northeast, central, western, and national areas are all greater than one except that    , the total nighttime lighting shows a downward trend on the national scale, which is inconsistent with the urbanization process in China. On the regional scale, it can be seen that from 2008 to 2009, the northeast region shows the opposite change of the total nighttime lighting with that of the national scale, while the remaining regions show a downward trend of the total nighttime lighting in these years. Moreover, from 2009 to 2010, the total nighttime lighting shows an increase trend, which was also inconsistent with the actual change of nighttime lights. The 2014 NPP/VIIRS data is much smaller in magnitude than the 2013 DMSP/OLS data, which results in the data drop from 2013 to 2014. Although NPP/VIIRS data is obtained from the same satellite, the total nighttime lighting also drops at national and regional scales from 2015 to 2016, and the total nighttime lighting drops at western regions from 2017 to 2018. Therefore, the uncorrected nighttime light data do not reflect the socioeconomic situation in China very well.  , the total nighttime lighting shows a downward trend on the national scale, which is inconsistent with the urbanization process in China. On the regional scale, it can be seen that from 2008 to 2009, the northeast region shows the opposite change of the total nighttime lighting with that of the national scale, while the remaining regions show a downward trend of the total nighttime lighting in these years. Moreover, from 2009 to 2010, the total nighttime lighting shows an increase trend, which was also inconsistent with the actual change of nighttime lights. The 2014 NPP/VIIRS data is much smaller in magnitude than the 2013 DMSP/OLS data, which results in the data drop from 2013 to 2014. Although NPP/VIIRS data is obtained from the same satellite, the total nighttime lighting also drops at national and regional scales from 2015 to 2016, and the total nighttime lighting drops at western regions from 2017 to 2018. Therefore, the uncorrected nighttime light data do not reflect the socioeconomic situation in China very well.  As shown in Figure 4b, it can be seen on the national scale from 2001 to 2018, the total nighttime lighting after calibration shows the trend of growth, in accordance with urbanization process in China. However, on the regional scale, from 2013 to 2014 and from 2017 to 2018, the total nighttime lighting in the northeast region decreased. The total nighttime lighting in the central region decreased from 2013 to 2015. Additionally, the total nighttime lighting in the western region decreased from 2014 to 2015. As shown in Figure 4c, it is the change of total nighttime lighting at the national and regional scale after the reference correction method. It can be seen on the national scale from 2001 to 2018, the total nighttime lighting shows the trend of growth, in accordance with urbanization process in China. However, on the regional scale, from 2013 to 2014, the total nighttime lighting in the eastern and northeast regions decreased.
NDI of the nighttime light data (including the DMSP/OLS and NPP/VIIRS data) from 2002 to 2018 were corrected under the correction method in this paper and the reference correction method has improved significantly and the interannual volatility of nighttime lighting data has been greatly reduced ( Figure 5). Compared with Figure 5b,c, it can be found that the interannual volatility under the correction method in this paper from 2013 to 2014 is lower than that of reference correction method. Therefore, the difference of nighttime light data within each province should be considered in the regression model of DMSP/OLS data and NPP/VIIRS data. The results of regression fitting for each province are better than those on the national scale.
It can also be seen that no matter the correction method in this paper and reference correction method can all improve the interannual fluctuation of nighttime light data ( Figure 6). Before correction, SNDI in northeast, central, western, and national areas are all greater than one except that As shown in Figure 4b, it can be seen on the national scale from 2001 to 2018, the total nighttime lighting after calibration shows the trend of growth, in accordance with urbanization process in China. However, on the regional scale, from 2013 to 2014 and from 2017 to 2018, the total nighttime lighting in the northeast region decreased. The total nighttime lighting in the central region decreased from 2013 to 2015. Additionally, the total nighttime lighting in the western region decreased from 2014 to 2015. As shown in Figure 4c, it is the change of total nighttime lighting at the national and regional scale after the reference correction method. It can be seen on the national scale from 2001 to 2018, the total nighttime lighting shows the trend of growth, in accordance with urbanization process in China. However, on the regional scale, from 2013 to 2014, the total nighttime lighting in the eastern and northeast regions decreased.
NDI of the nighttime light data (including the DMSP/OLS and NPP/VIIRS data) from 2002 to 2018 were corrected under the correction method in this paper and the reference correction method has improved significantly and the interannual volatility of nighttime lighting data has been greatly reduced ( Figure 5). Compared with Figure 5b,c, it can be found that the interannual volatility under the correction method in this paper from 2013 to 2014 is lower than that of reference correction method. Therefore, the difference of nighttime light data within each province should be considered in the regression model of DMSP/OLS data and NPP/VIIRS data. The results of regression fitting for each province are better than those on the national scale.
It can also be seen that no matter the correction method in this paper and reference correction method can all improve the interannual fluctuation of nighttime light data ( Figure 6). Before correction, SNDI in northeast, central, western, and national areas are all greater than one except that in the eastern region. After correction, SNDI are about 50% of that before correction. After correcting the nighttime light data with reference correction method, when establishing the regression model of DMSP/OLS and NPP/VIIRS data to fit the NPP/VIIRS data, then fitting the NPP/VIIRS data to the DMSP/OLS data scale in 2014, it will lead to large differences between the correction nighttime light data and DMSP/OLS data in 2013 among some provinces. As a result, the corrected SNDI of the reference correction method on the regional scale is slightly larger than that of the calibration method in this paper. When total nighttime lighting on the national scale is calculated, the influence of the difference in fitting of nighttime light data within the provincial region will be reduced.

Modeling Results for the Spatiotemporal Dynamics of the Socioeconomic Data
Using the extended temporal coverage nighttime light data, the spatiotemporal dynamics of the GDP from 2001 to 2018 at both the national and provincial scales were modeled. Three regression models were applied: a linear model, a quadratic polynomial model, and a power function model. Then, the best fitting equation was chosen to carry out fitting analysis, compare the fitting accuracy after the correction, and verify the accuracy of the correction algorithm.

Regression Result Analysis at National Scale
Referring to the approach proposed by Zhu et al. [22], through the three regression models mentioned above, the long-term relationships between the corrected TNL and GDP (Table 1), EPC ( Table 2) in China from 2001 to 2018 were analyzed on time series at the national scale. The R 2 values obtained from three models between the corrected TNL and GDP are 0.972 (linear), 0.984 (quadratic), and 0.981 (power), respectively, which all represent high correlations and have no marked difference. However, the MARE values vary among the different regression models. For example, the MARE value is 16.679% in the linear model and the MARE value is 10.056% in the quadratic polynomial model, whereas the MARE value is 8.518% in the power function model. This indicates that some details in the relationship between the TNL and GDP may be ignored when the optimal model is selected solely based on the R 2 values. Thus, these results indicate that all of the three regression models can appropriately reflect the long-term GDP dynamics in China using the nighttime light data from 2001 to 2018 and that the power function model is optimal at the national scale. Moreover, all of the absolute RE values of the power function model employing the extended temporal data are less than or equal to 19.94%. Table 1. Accuracy of the estimated statistical gross domestic product (GDP) on time series at the national scale through the three regression models. Note: GDP represents the statistical gross domestic product, PG represents the predicted GDP, RE represents the relative error in percentage, and MARE represents the mean absolute relative error in percentage. The unit of the GDP is billion RMB. The R 2 values obtained from three models between the corrected TNL and EPC are 0.990 (linear), 0.990 (quadratic), and 0.985 (power), respectively, which all represent high correlations and have no marked difference. Moreover, the MARE value is 4.702% in the linear model and the MARE value is 4.655% in the quadratic polynomial model, whereas the MARE value is 4.950% in the power function model. Thus, these results indicate that all of the three regression models can appropriately reflect the long-term EPC dynamics in China and that the quadratic polynomial model is optimal at the national scale. Moreover, all of the absolute RE values of the quadratic polynomial model employing the extended temporal data are less than or equal to 11.64%. Therefore, modeling the long-term dynamics of socioeconomic data in China from 2001 to 2018 using the extended temporal data is reliable. Table 2. Accuracy of the estimated statistical electric power consumption (EPC) on time series at the national scale through the three regression models. Note: EPC represents the statistical electric power consumption, PE represents the predicted EPC, RE represents the relative error in percentage, and MARE represents the mean absolute relative error in percentage. The unit of EPC is billion·kilowatt·h −1 . The predicted GDP and EPC were estimated by the established power function and quadratic polynomial model to further verify the accuracy of GDP modeling (Table 3) and EPC modeling (Table 4) based on the time series.

Year
The R 2 value obtained between the corrected TNL corrected in this paper and GDP is 0.981, whereas the R 2 value of that corrected under the reference correction method is 0.951. Moreover, the MARE value of the model between the nighttime light data corrected in this paper and GDP is 8.518%, while the MARE of the model between the nighttime light data corrected under the reference correction method and GDP is 24.761% (Table 3). Therefore, the correction method in this paper and the reference correction method can be implemented to GDP on a national scale for dynamic modeling for a long time. Moreover, the calibration method in this paper is slightly better than the reference correction method in modeling accuracy. By comparing the modeling results of nighttime light data before 2013 and that after 2013, it can be found that the MARE of the model between the nighttime light data corrected in this paper and GDP before 2013 is 8.365%, and the MARE of the model between the nighttime light data corrected in this paper and GDP after 2013 is 8.913%. While, the MARE of the model between the nighttime light data corrected under the reference correction method and GDP before 2013 is 21.720%, and the MARE of the model between the nighttime light data corrected under the reference correction method and GDP after 2013 is 32.666%. It can be found that the correction method in this paper has a certain contribution to the saturation correction of DMSP/OLS data compared with the reference correction method, which makes the NTL data better represent GDP. The R 2 value obtained between the corrected TNL in this paper and EPC is 0.990, whereas the R 2 value of that corrected under the reference correction method is 0.976. Moreover, the MARE value of the model between the nighttime light data corrected in this paper and EPC is 4.655%, while the MARE of the model between the nighttime light data corrected under the reference correction method and EPC is 7.834% (Table 4). Therefore, the correction method in this paper and the reference correction method can be implemented to EPC on a national scale for dynamic modeling for a long time. Moreover, the calibration method in this paper is slightly better than the reference correction method in modeling accuracy. By comparing the modeling results of nighttime light data before 2013 and that after 2013, it can be found that the MARE of the model between the nighttime light data corrected in this paper and EPC before 2013 is 5.238%, and the MARE of the model between the nighttime light data corrected in this paper and EPC after 2013 is 3.141%. While the MARE of the model between the nighttime light data corrected under the reference correction method and EPC before 2013 is 9.607%, and the MARE of the model between the nighttime light data corrected under the reference correction method and EPC after 2013 is 3.223%. By comparing Tables 3 and 4, it can be found that corrected nighttime light data has better modeling ability for EPC (4.655%) on a national scale for dynamic modeling for a long time than for GDP (8.518%), and can be used to predict EPC more accurately.

Regression Result Analysis Based on Provincial Units
The long-term relationships between the corrected TNL and GDP (Table 5), EPC (Table 6) in China from 2001 to 2018 were analyzed using the three regression models, and the R 2 and MARE values for each of the three models for each year in provincial units were calculated accordingly. As shown in Table 5, the mean R 2 values obtained from three models between the corrected TNL and GDP are 0.8068 (linear), 0.8258 (quadratic), and 0.7386 (power), respectively. Moreover, the MARE value is 45.239% in the linear model and the MARE value is 55.708% in the quadratic polynomial model, whereas the MARE value is 38.599% in the power function model. In addition, the MARE value of the power function is lower than that of the other two functions for each year. Therefore, considering the comparison results comprehensively, the power function model was chosen as the fitting function of GDP and TNL.  As shown in Table 6, the mean R 2 values obtained from three models between the corrected EPC and TNL are 0.8883 (linear), 0.8927 (quadratic), and 0.7315 (power), respectively. Moreover, the MARE value is 30.350% in the linear model and the MARE value is 29.319% in the quadratic polynomial model, whereas the MARE value is 30.830% in the power function model. In addition, the MARE value of the power function is lower than that of the other two functions almost every year. Therefore, considering the comparison results comprehensively, the quadratic polynomial model was chosen as the fitting function of EPC and TNL. It can be seen that the fitting function model in provincial units is the same to that in the national scale. The relationship between TNL before correction and under two correction methods and GDP and EPC in provincial units in different years was evaluated by power function model (Table 7) and quadratic polynomial model (Table 8), respectively.
As shown in Table 7, the minimum R 2 value of the model between TNL of DMSP/OLS images before correction and GDP appeared in 2012, which was 0.600. The maximum R 2 value of the model between TNL of DMSP/OLS images before correction and GDP appeared in 2003, which was 0.786. However, after 2014, R 2 value of the model between TNL of NPP/VIIRS images before correction and GDP was lower than 0.5. The minimum R 2 value appeared in 2017, which was only 0.165. The maximum R 2 value appeared in 2016, which was only 0.417. The R 2 values of the model between TNL of NPP/VIIRS annual composite in 2015 and GDP and EPC were 0.792 and 0.834, respectively. Therefore, there were still short-term lights in NPP/VIIRS annual composite images with background noise removed, so the fitting results between TNL and GDP and EPC were worse than those between annual composite images with unstable pixels removed and GDP and EPC.
Under the correction method in this paper, the R 2 values of the model between TNL and GDP in provincial units were higher than 0.65. The minimum R 2 value that appeared in 2017 was 0.661. While the maximum R 2 value that appeared in 2006 was 0.798. The average R 2 value was 0.7386. The regression result of the correction results under the reference correction method and GDP was also better than the regression result before correction. The minimum R 2 value that appeared in 2016 and 2017 was 0.490. While the maximum R 2 value appeared in 2006, which was 0.771. The average R 2 value was 0.6650. At the same time, the R 2 and MARE values of the model between TNL of the F162010 radiometric calibration nighttime light image and GDP were 0.795 and 34.97%, the R 2 and MARE values of that under the correction method were 0.746 and 38.01%, whereas the R 2 and MARE values of that under the reference correction method were 0.697 and 42.17%. Therefore, it can be found that the correction method in this paper plays a greater role in reducing the "saturation effect" of the image than the reference correction method. By comparing the results of 2015, it can be seen that the R 2 and MARE values (0.682 and 42.38%) of the model after correction in this paper is lower than that (0.792 and 32.23%) of annual composite image, but better than that (0.509 and 52.45%) corrected under the reference correction method. As shown in Table 8, under the correction method in this paper, the minimum R 2 value between TNL after correction in this paper and EPC that appeared in 2013 was 0.877. While the maximum R 2 value that appeared in 2016 was 0.909. The average R 2 value was 0.8927. Under the reference correction method, the minimum R 2 value that appeared in 2015 was 0.855. While the maximum R 2 value appeared in 2014, which was 0.892. The average R 2 value was 0.8792. At the same time, the R 2 and MARE values of the model between TNL of the F162010 radiometric calibration nighttime light image and EPC were 0.861 and 33.36%, the R 2 and MARE values of that under the correction method were 0.893 and 30.31%, whereas the R 2 and MARE values of that under the reference correction method were 0.874 and 30.68%. Therefore, it can be also found that the correction method in this paper plays a greater role in reducing the "saturation effect" of the image than the reference correction method. By comparing the results of 2015, it can be seen that the R 2 and MARE values (0.900 and 26.82%) after correction in this paper is better than that (0.855 and 26.77%) corrected under the reference correction method, and better than that (0.834 and 29.96%) of the annual composite image.
By comparing the relationship between TNL and EPC under different treatments, it can also be found that the regression results under the correction method introduced in this paper are better than those under the reference correction method, and the regression results under two correction methods are all better than the regression results before the correction. Therefore, the nighttime light data corrected by the correction method in this paper can reflect GDP and EPC better than in provincial units than that corrected under the reference correction method and that before correction.

Regression Result Analysis Based on Time Series
The long-term relationships between the corrected TNL and GDP (Table 9), EPC (Table 10) in China from 2001 to 2018 were analyzed using the three regression models, and the R 2 and MARE values for each of the three models for each province were calculated accordingly. As shown in Table 9, the mean R 2 values obtained from three models between the corrected TNL and GDP are 0.9529 (linear), 0.9664 (quadratic), and 0.9671 (power), respectively. Furthermore, the R 2 values of the three models are similar among most of the provinces, whereas the MARE values vary substantially with the different models for each province. It can be concluded that at the provincial level, all of the three models are reliable for modeling the long-term GDP dynamics the nighttime light data from 2001 to 2018. In general, the power function model is the best-fitting model, with the minimum mean MARE value of 11.181%, whereas the MARE value is 19.883% in the linear model and the MARE value is 13.361% in the quadratic polynomial model, which is supported by previous studies [13,22,78]. As shown in Table 10, the mean R 2 values obtained from three models between the corrected EPC and TNL are 0.9628 (linear), 0.9713 (quadratic), and 0.9675 (power), respectively. Moreover, the MARE value is 8.468% in the linear model and the MARE value is 6.915% in the quadratic polynomial model, whereas the MARE value is 7.066% in the power function model. It can be seen that the fitting function model for each province is the same as that in provincial units and the national scale.  To further ascertain the most effective model for each province, the optimal models were chosen through a comparison of the R 2 and MARE values of each model. The quadratic polynomial model is found to be optimal in eight provinces (Fujian, Gansu, Hainan, Jiangsu, Shaanxi, Shanghai, Xinjiang, and Yunnan) while the power function model is optimal for modeling the dynamics of the GDP from 2001 to 2018 for the remaining provinces. The quadratic polynomial model is found to be optimal in 12 provinces (Fujian, Gansu, Guizhou, Hainan, Heilongjiang, Hubei, Hunan, Jilin, Shaanxi, Shanghai, Xinjiang, Zhejiang), the linear model is optimal in five provinces (Beijing, Guangdong, Hebei, Liaoning, and Yunnan) while the power function model is optimal for modeling the dynamics of the EPC from 2001 to 2018 for the remaining provinces.
Consequently, although the power function model and quadratic polynomial model is generally the best-fitting model for long-term modeling of GDP and EPC dynamics using the nighttime light data from 2001 to 2018 for each province in time series, the optimal model varies for different provinces of China. The relationship between TNL before correction and under two correction methods and GDP and EPC in provincial units in different years was evaluated by the corresponding best-fitting model, respectively.
Referring to the method proposed by Shi et al. [61], the minimum R 2 value of the model between TNL corrected in this paper and GDP was 0.927, the maximum R 2 value was 0.992, and the average R 2 value was 0.9667 (Table 11). According to the reference correction method, the minimum R 2 value of the model between TNL corrected and GDP was only 0.001, the maximum R 2 value was 0.983, and the average R 2 value was 0.8325. Besides, the minimum R 2 value of the model between TNL corrected in this paper and EPC was 0.930, the maximum R 2 value was 0.993, and the average R 2 value was 0.9720. According to the reference correction method, the minimum R 2 value of the model between TNL corrected and EPC was only 0.001, the maximum R 2 value was 0.988 and the average R 2 value was 0.8286. In addition, the TNL-GDP model and TNL-EPC model of 30 provinces obtained in this paper have all passed the F-test at the level of 0.001, while the TNL-GDP model and TNL-EPC model of 28 provinces besides Tianjin and Shandong obtained under the reference correction method passed the F-test at the level of 0.001.  (Table 12) and the predicted EPC (Table 13) in time series and in provincial units under the two correction methods was further analyzed. When modeling based on time series, MARE values of the predicted GDP under the corrected NTL data in this paper are 16.82% (Jiangxi) and 5.12% (Xinjiang), while MARE values of the predicted GDP under that of the reference correction method are 39.12% (Beijing) and 7.55% (Xinjiang) ( Table 12). The correction method in this paper is better than the reference correction method in estimating the GDP of specific provincial units. For the provincial units, the MARE values of the uncorrected nighttime light data in predicting the GDP of Qinghai and Anhui are 210.44% and 13.55%. At the same time, MARE values of the predicted GDP under the corrected nighttime light data in this paper are 114.41% (Hainan) and 5.06% (Guizhou), while MARE values predicted GDP under that of the reference correction method are 115.76% (Xinjiang) and 6.37% (Jilin). Therefore, the performance prediction GDP under nighttime light data processed by the two correction methods is better than that before correction. When forecasting GDP of Hainan, Ningxia, Qinghai, and Xinjiang based on the model on provincial unit, MAREs are all large under two correction methods and before correction because it is difficult to distinguish the real lights in urbanization areas from other artificial structures or some special surface reflected lights. Based on the model of nighttime light data and GDP established by time series, the MARE value under the correction method in this paper is 10.877%, while the MARE value under the reference correction method is 19.249%. It is shown that the correction method in this paper is also better than the reference correction method in the long-term dynamic modeling of GDP on the provincial scale. Based on the model of nighttime light data and GDP established by provincial units, the MARE value under the correction method in this paper is 38.599%, the MARE value under the reference correction method is 43.560%, while the MARE before correction is 52.785%. The MARE value of the model based on provincial units under the two correction methods is far greater than that based on time series. When modeling in provincial units, the nighttime light data of some provincial units are not real data, so the correlation between nighttime light data and GDP will be reduced during modeling, thus reducing the accuracy of GDP prediction. When modeling in time series, it is not necessary to consider the differences in this part. Therefore, when predicting the GDP of specific provincial units, the regression model based on time series is preferred.
The MARE value of the model between the nighttime light data under the correction method in this paper and EPC based on time series is 6.435%, while the MARE value under the reference correction method is 12.178% (Table 13). It is shown that the ability of predicting EPC for dynamic modeling for a long time in the provincial scale under the correction method in this paper is better than that under the reference correction method. Based on the model of nighttime light data and EPC established by provincial units, the MARE value under the calibration method in this paper is 31.487%, the MARE value under the reference correction method is 29.600%, while the MARE value before the calibration is 35.612%. The MARE value based on the modeling of provincial units under the two correction methods is also much greater than that based on the modeling of time series. Therefore, when predicting the EPC of specific provincial units, the regression model based on time series is also optimized. The MARE value is further divided into three categories [61]: 0-25% is high accuracy, 25-50% is medium accuracy, and more than 50% is inaccuracy (Table 14).
In the TNL-GDP and TNL-EPC models established based on time series, the correction method introduced in this paper has extremely significant prediction ability with prediction accuracy of 100%, which is higher than the prediction accuracy (76.67% and 90.00%) of the reference correction method. The prediction inaccuracy in this paper is also less than the prediction inaccuracy (6.67% and 6.67%) of the reference correction method. Compared with the prediction ability of GDP and EPC under the correction method introduced in this paper and the reference correction method, it can be found that the ability of the correction method introduced in this paper in the prediction of GDP and EPC is also better than that of the reference correction method.

Discussion
DMSP/OLS data and NPP/VIIRS data are the most widely used to obtain the nighttime light data for monitoring human socioeconomic activities and natural phenomena [6]. Light saturation in the city center, blooming effect, and lack of airborne calibration are the major obstacles to the application of DMSP/OLS data. In addition, when building a long time series of nighttime light datasets, the inconsistency between NPP/VIIRS data and DMSP/OLS data is more serious than that between the different products of the DMSP/OLS data [50]. Although, scholars have done a lot of research on the processing of DMSP/OLS Data and NPP/VIIRS data, but there are still some limitations in the existing research results.
Aiming at the problems in the current research, this paper proposes the WorldPop and the enhanced vegetation index adjusted nighttime light (WEANTL) using EVI and WorldPop data based on the principle that urban characteristics (including nighttime light) should be inversely proportional to vegetation coverage [38] and directly proportional to population distribution data [53,54], so as to obtain the nighttime light dataset from 2001 to 2018 to represent economic data reliably and accurately. The result shows that the extended nighttime light data set has good quality and reliable time consistency. Considering that the use of nighttime light data to reflect social and economic development is carried out by constructing the statistical relationship between TNL and socioeconomic parameters [7,9,35,74,75], this paper also evaluates the accuracy of the correction method proposed in this paper in predicting socioeconomic parameters by constructing a regression model between TNL and socioeconomic parameters.
In the saturation correction of DMSP/OLS data, the existing research did not distinguish the spatial differences of the saturated pixels and carried out the saturation correction of the unsaturated pixels, thus failed to effectively correct the saturation of the DMSP/OLS data. When fitting DMSP/OLS data and NPP/VIIRS data, some scholars used the same fitting function to fit the data of 30 provinces, without considering the data differences of each province, thus reducing the ability of night light data to represent economic data. Compared with the previous models [63,77], it can effectively distinguish the spatial differences of the nighttime light data because EVI and WorldPop data was used to realize the saturation correction of DMSP/OLS data.
In this paper, when fitting DMSP/OLS data and NPP/VIIRS data, considering the differences of economic development of each province, linear regression models were developed for each province. It can be seen from the above paper, the TNL-GDP model and TNL-EPC model of 30 provinces obtained in this paper all passed the F-test at the level of 0.001, while the TNL-GDP model and TNL-EPC model of 28 provinces besides Tianjin and Shandong obtained under the reference correction method passed the F-test at the level of 0.001. Tianjin and Shandong failed to pass the F-test mainly because when fitting DMSP/OLS and NPP/VIIRS data, the model was fitted on all county-level scales in China in 2012 and 2013. NPP/VIIRS data from 2014 to 2018 were fitted to DMSP/OLS data according to the fitting model, which can ensure that the DMSP/OLS data in later years after fitting on a national scale is almost larger than that in previous years. However, there is no guarantee that DMSP/OLS data in later years after fitting of all the provincial units in the provincial scale is greater than the DMSP/OLS data in previous years. For example, the DMSP/OLS data of Beijing in 2014 or 2015 after fitting is less than that in 2013, leading to the R 2 values of the TNL-GDP model and TNL-EPC model being below the average. The DMSP/OLS data of Tianjin and Shandong after fitting from 2014 to 2018 is all less than that in 2013. As a result, the established regression model could not pass the F test at the 0.001 level. Therefore, it is necessary to establish the fitting equation between the DMSP/OLS and NPP/VIIRS data in 2012 and 2013 for the specific research area when studying the research area on and below the provincial scale.
Compared with the existing research, the main contributions of this paper are as follows: (1) DMSP/OLS image is divided into saturation and unsaturated regions, WEANTL was put forward to correct the image data in the saturation region by using EVI and WorldPop data of the corresponding year, and pixel points in the unsaturated region kept the original value. Therefore, the saturation correction of DMSP/OLS image data from 2001 to 2013 is realized by using EVI and WorldPop data. As a correction method for fitting DMSP/OLS data and NPP/VIIRS data to obtain a long time series of nighttime light data sets for continuous monitoring of society, this method has made some contributions, but there are still some limitations. Firstly, although the prediction accuracy of the proposed correction method in establishing long-term GDP and EPC dynamic models is better than that of the reference correction method, the saturation problem in the DMSP/OLS data and the noise data in NPP/VIIRS data are not removed completely. These problems reduce the accuracy of long time series nighttime light datasets and the prediction accuracy of establishing long-term GDP and EPC dynamic models. Secondly, considering that regression models (quadratic polynomial model and power function model) have been applied in this paper to model socioeconomic activities using nighttime light data, more complex models need to be developed to predict socioeconomic activities more accurately. Therefore, with the updating of NPP/VIIRS data and more sources of nighttime light data, such as luojia1-01 satellite, new methods for multisensor nighttime light data to model social and economic dynamic situations can be developed in future research.

Conclusions
This paper proposes an invariant region method based on EVI and WorldPop data to achieve intercalibration and saturation correction of DMSP/OLS data, and uses a regression model to achieve the fitting of DMSP/OLS data and NPP/VIIRS data, and constructs a long time series of night light images from 2001 to 2018. Comparing the pixel results processed by three data processing methods (uncorrected, calibration method introduced in this paper and the reference correction method), the ability to estimate GDP and EPC based on time series on the national scale and the ability to estimate GDP and EPC based on time series and on provincial units under two correction methods (calibration method introduced in this paper and the reference correction method) were carried out.
In this research, we found that (1) after calibration, the NTL data can guarantee the growth trend on national and regional scales, and the interannual volatility of the corrected NTL data is lower than that of the uncorrected NTL data; are 0.9720 and 6.435%; the established TNL-GDP and TNL-EPC models with 30 provinces data all passed the F-test at the 0.001 level; (5) the prediction accuracy of GDP and EPC on time series was nearly 100%. Therefore, the correction method in this paper can reliably and accurately estimate GDP and EPC on multiple scales.
Author Contributions: P.L. conducted the research. Q.W. revised the paper and guided the research, D.Z. and P.L. were responsible for collection data, creating the figures, and revising the paper. Y.L. revised and improved the paper. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A