Intercalibration of Luojia1-01 and Suomi-NPP-VIIRS Monthly Nighttime Light Composite Using a Spatial-Temporal Residuals Correction Random Forest Model

Satellite-derived nighttime light (NTL) data from the Suomi National Polar-orbiting Partnership Visible Infrared Imaging Radiometer Suite (VIIRS) and Luojia1-01 (LJ1-01) have been widely used in multiple urban aspects. However, the inconsistency in NTL intensity (the same spatial coverage and month) between the two sensors impedes their integration on finer scales. To alleviate this issue, a novel model was developed to simulate LJ1-01-like NTL in March 2019, considering Guangzhou main city as the study area, Shenzhen, Dongguan, and Huizhou as validation cities to confirm the spatial simulation ability of the model. First, the spatial features originated from VIIRS NTL and geographical data to enhance the prior spatial information for NTL. Second, the temporal features were derived from the VIIRS NTL time series between 2018 and 2020, such as interannual and intra-annual trends, to describe the temporal fluctuation of NTL. Third, a spatial-temporal residuals correction random forest (STRCRF) regression model was proposed to generate LJ1-01-like NTL validation with ten-fold cross-validation (CV). It established the nonlinear relationship between LJ1-01 NTL intensity and spatial-temporal features at the pixel level. It generated LJ1-01-like NTL data with a CV-R2 of 0.982 and a CV-RMSE of 3469. The validation accuracy with R2 ranged from 0.918 to 0.945, and RMSE varied from 4139 to 7298 in Shenzhen, Dongguan, and Huizhou. The evaluations of LJ1-01-like NTL based on the statistic analyses, spatial patterns, and profile analyses indicated the effectiveness of the STRCRF model in improving NTL intensity consistency at the pixel and landscape level. Moreover, the STRCRF model outperformed the previous related studies and improved the CV-R2 by 0.4. Our approach offers the potential for simulating LJ1-01-like NTL by taking advantage of the prior spatial and temporal information and the ensemble regression with residuals correction. This method could contribute to a higher quality NTL dataset for mapping human activities and monitoring the urban environment.


I. INTRODUCTION
T HE nighttime light (NTL) has been used as a proxy indicator of human activities because it can detect the intensity and spatial extent of lights at night [1]. The NTL were collected from the Operational Line-scan System (OLS) of Defense Meteorology Satellite Program (DMSP), subsequent Visible Infrared Imaging Radiometer Suite (VIIRS) of the Suomi National Polar-orbiting Partnership (NPP) satellite, and the newly Luojia1-01 (LJ1-01) have triggered many analyses in multiple urban aspects. The NTL provides a unique view to map the urbanization process [2], [3], [4], estimate and spatialize the socio-economic factors [5], [6], monitor the light pollution [7], [8], disasters [9], [10], and conflict [11], [12].
Most of the studies based on NTL (i.e., DMSP and VIIRS) were done to a large spatial extent, which was not conducive to examining the intraurban variation in NTL intensity [13]. Fortunately, some studies have proved that LJ1-01 NTL has a more excellent performance in extracting urban areas [14], [15], investigating light pollution [7], and estimating socio-economic activities [16] at a finer scale than VIIRS NTL data. However, it is challenging for LJ1-01 NTL to explore dynamic urban changes since this dataset was obtained in June 2018 [7]. Some studies affirmed the usefulness of VIIRS NTL in urban spatial-temporal changes [17], [18]. Therefore, it is crucial to combine the advantages of VIIRS and LJ1-01 NTL exploring the urban light characteristics at a finer scale.
Tackling the inconsistency between VIIRS and LJ1-01 NTL is a critical challenge in their integration into urban application. This issue can be ascribed to two reasons: 1) The sensor performance differences between VIIRS and LJ1-01. LJ1-01 present a significant advance over VIIRS in spatial and radiometric resolution [7], [14]. The local overpass time of VIIRS is about 1:30 AM, [19] when there is generally less illuminated activity This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ occurring on the ground than during the 22:00 overpass time of LJ1-01 [20]. Meanwhile, the low light detection limits of the two sensors are also different [14]. 2) The inherent flaws of both satellite sensors. For example, Version 1 (V1) VIIRS monthly datasets are commonly used; however, they are easily affected by background noise and seasonal effects [21]. LJ1-01 has some geographical position errors and abnormal values that need to be corrected. Thus, there are some inconsistencies between the two sensors at the same spatial coverage, i.e., 1) NTL intensity for the same month; 2) the trend and variability of NTL time series; and 3) application results, such as extraction of urban areas [14], [15]. This study attempts to solve the first issue by generating the high-quality LJ1-01-like NTL based on VIIRS and LJ1-01 NTL for the same month.
To address the above issues between two sets of NTL data, many useful approaches have been proposed. The intercalibration is one of the mainstream methods. This method assumes that the observation of light intensity of the same extent has a slight difference owing to the difference in overpass time. Then it establishes a robust regression relationship to generate consistent NTL, including characteristics of two NTLs. There are three types of intercalibration, i.e., intercalibration at the same scaling, upscaling, and downscaling.
First, the same scaling intercalibration was used in DMSP annual composite NTL data with six sensors [22]. The pseudoinvariant features (PIF) paradigm was primarily adopted. These studies assumed that the digital number (DN) value at the pixel level would hardly change in stable and highly developing areas over several years. Then, the remaining pixels were calibrated to the reference pixels with the empirical model. Pandey et al. [26] reviewed six common intercalibration methods [4], [12], [22], [23], [24], [25] and systematically analyzed their calibration performance at global, regional, and pixel level [26]. However, this paradigm cannot reveal the actual development in rapid urbanization areas, i.e., China and India.
Second, the upscaling intercalibration mainly oriented DMSP and VIIRS NTL to generate DMSP-like NTL. There are two methods summarized. 1) Image preprocessing methods of VIIRS to narrow the gap between VIIRS and DMSP NTL, including Gaussian low pass [11], logarithmic transformation [27], and breaks for additive seasonal and trend time series decomposition algorithm [28]. Then the regression functions, i.e., power function, sigmoid function, and geographically weighted regression, were selected to simulate the DMSP-like data. 2) Based on pixels with stable relationships and the PIF paradigm to simulate DMSP-like NTL [29].
Third, the downscaling intercalibration was adopted between VIIRS and DMSP NTL to generate VIIRS-like NTL. Some researchers adopted power function [30] and machine learning (ML) algorithms, i.e., the auto-encoder [31], random forest (RF), and multi-layer perception [32], to achieve this goal. Compared with traditional statistical models, ML has vastly improved calibration performance due to exploring the complex relationships between different NTLs.
Despite these creative efforts, there are some limitations in current intercalibration methods. The spatial information of the land surface and temporal information of NTL is less considered the independent variables in the regression model. This may lead to unsatisfactory results [31], [33]. In addition, the residuals originating from the regression model are neglected mostly. Some researchers attempted to correct residuals of simulated products and had excellent performance [28]. Finally, most researchers were committed to calibrating DMSP and VIIRS NTL, lacking studies aimed at intercalibration between LJ1-01 and VIIRS NTL.
Due to the gap in existing studies, a novel method for generating LJ1-01-like NTL with high spatial and radiometric resolution integrates VIIRS and LJ1-01 NTL using the ML approach. The proposed approach intends to fully use geographical data related to NTL and temporal variability of VIIRS time-series data to obtain significant spatial-temporal features and to establish the spatial-temporal residuals correction random forest (STRCRF), intercalibration model, to produce high-quality LJ1-01-like NTL.

A. Study Area
Guangzhou is the capital city of Guangdong province. It is located centrally in Guangdong-Hong Kong-Macao Greater Bay Area (GBA) within the Pearl River Delta (see Fig. 1). According to the Statistical Yearbook of Guangzhou and previous studies, there was high development from 2010 to 2019 [34], [35]. The Guangzhou main city (GZMC) is the central area of Guangzhou, which consists of six districts, i.e., Liwan(LW), Yuexiu(YX), Haizhu(HZ), Tianhe(TH), Baiyun(BY), and Huangpu(HP), with a land area of 1472.93 km 2 and more than nine million people in 2019. Its elevation ranges from 0 to 255 m. It is dominated by a southern subtropical humid monsoon climate and complex land cover physical condition. The GZMC also represents the most active area of Guangzhou in terms of human activities. Thus, its NTL intensity shows apparent spatial characteristics, which provides an ideal study area to test the proposed method.
Moreover, three cities were selected to examine the robustness of the proposed method. The city network of the three cities shows different levels in the GBA; Shenzhen, Dongguan, and Huizhou are regarded as the core city, the semi-peripheral city, and the peripheral city, respectively [36].  [37]. It can obtain high-quality NTL imagery with a spatial resolution of 130 m and a swath of 250 km. The data are free to download at the High-Resolution Earth Observation System of the Hubei Data and Application Center. 1 In general, the quality of LJ1-01 NTL will improve after the radiometric calibration, which enlarges the gap with VIIRS NTL. Thus, LJ1-01 NTL with DN value covering the study area on March 11, 2019, was used in this experiment for modeling and analysis.
2) VIIRS: The VIIRS NTL datasets (Version 1) with a spatial resolution of 500 m and a panchromatic band ranging between 505 and 890 nm have been available in month composites since April 2012 from NOAA/NCEI website. Two datasets, including "vcm" and "vcmsl" are available for users. Contaminated by lightning, lunar illumination, and cloud cover before averaging are already excluded in both datasets [21], [38]. The main difference between the two datasets is that the "vcm" files remove stray-light-contaminated data, whereas the "vcmsl" files include the stray-light corrected data. Therefore, the difference is pretty slight in most regions except for the high latitudes [39]. The "vcm" version of the DNB monthly composite imagery with "avg_rad" band was used for observations and analysis. These datasets with pixel values represent radiance in nanoWatts (nW)/(cm 2 ·sr). VIIRS monthly composite products from 2018 to 2020 (a total of 36 months) were used in this study. The comparison of parameters for VIIRS and LJ1-01 NTL data is shown in Table I.
3) Auxiliary Data: Related auxiliary data were used in this study, including raster imageries, vector data, and statistical data. The raster imageries include the 30 m resolution Landsat 8 operational land imager (OLI) multispectral bands on January 30, 2019, the globeland30 with ten classes in 2020,  TABLE II  GEOGRAPHICAL AUXILIARY DATA INFORMATION TABLE   and Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (DEM) data in 2019, to represent the physical environment within the study area. The vector data are road networks from Open Street Map to show the socio-economic condition indirectly. The statistical data are the permanent population (POP) in six districts of GZMC in 2019 used to validate the performance of LJ1-01-like NTL. The above data source is summarized in Table II.

C. Data Preprocessing
First, all data were clipped with the boundary of GZMC and projected into Albers Conical Equal Area to keep the area unchanged. Due to different spatial resolutions, all data used in this study were uniformly resampled to the same spatial size (130 m) using the cubic convolution interpolation method.
According to the related study [7], the positioning accuracy of LJ1-01 NTL ranges from 490 to 930 m. Consequently, precise geometric correction is necessary for LJ1-01 NTL to be accurately modeled. Owing to the high spatial resolution of LJ1-01 NTL, the road networks can be observed. Thus, accurate ground control points (GCPs) can be collected using road intersections. Seven GCPs were selected to achieve the geometric correction [root mean square error (RMSE) = 0.48 pixel] using the second-degree polynomial regression model [see Fig. 2(a)]. The geolocating performance of LJ1-01 NTL imagery in GZMC is shown in [see Fig. 2 Besides, the LJ1-01 NTL needs to remove high-value noise owing to high sensitivity and extensive dynamic range [7]. We set the highest DN value of airports in Guangdong Province as the threshold; the process method can be expressed as in (1). In the end, the LJ1-01-corrected NTL on March 11, 2019, in GZMC was obtained.
where DN i and DN thr denote DN value in the ith pixel and the maximum value, respectively.

III. METHODOLOGY
In general, there are apparent tradeoffs between spatial and temporal resolution in the above NTL datasets. The intercalibration model should take these factors into account. First, NTL intensity varies dramatically in space and time due to complex geographical conditions [21] and diverse human activities   [17]. Second, the intercalibrated NTL has residuals deducing some unsatisfactory results [28]. An STRCRF model is, thus, developed to consider the influence of spatial and temporal features with residuals correction improving performance based on the RF algorithm [see Fig. 3(a)]. Based on the experience of previous studies [40], three models, i.e., the spatial random forest (SRF) only with spatial features, the temporal random forest (TRF) only with temporal features, and spatial-TRF (STRF) with spatial-temporal features, are used to prove the superiority of the STRCRF model.
The proposed methodology is comprised of five parts [see Fig. 3(b)]. First, two spatial analysis methods, i.e., kernel density and slope analysis, were adopted to generate spatial features related to NTL. Second, harmonic analysis (HA) was used to generate harmonic coefficients. Those coefficients were regarded as temporal features. Third, correlation and collinearity analyses were used to select significant and noncollinear features. Fourth, the STRCRF model was developed as an intercalibration method to simulate high-quality LJ1-01-like NTL in March 2019. At last, three evaluation methods, i.e., statistic analyses, spatial patterns, and profile analyses, were adopted. Below are details of the methodology.

A. Generating Spatial Features Based on Geographical Data
The geographical data related to NTL were processed to obtain the spatial features. According to Landsat 8, OLI surface reflectivity bands computed the Normalized Difference Vegetation Index (NDVI). The DEM data were adopted with slope analysis to generate slope data. Besides, five land cover raster datasets, i.e., cropland, grassland, forestland, impervious surface percentage (ISP), and water, were transformed into polygons. The 130 m GZMC fishnet points being inside the landcover polygons were selected to obtain the land cover points. Then, those points data and road networks with rank selection (Primary, Secondary, Trunk, and Trunk_link) were performed on the kernel density. In the end, the NDVI, slope, DEM, road density (RD), cropland density (CD), forestland density (FD), grassland density (GD), water density (WD), impervious surface percentage density (ISPD), and 130 m VIIRS radiance intensity were combined to generate spatial features where Band NIR , Band Red denote the near infrared and the red band of the Landsat 8 OLI multispectral imagery, respectively where f (x, y) is the kernel density in the pixel (x, y); h >0, h denotes bandwidth, which is the key factor of smoothing raster surface. We set the value of h to 130 m in correspondence with the raster data resolution; n is the number observation; K h denotes kernel function; d i denotes the distance between pixel (x, y) and ith observation location.

B. Generating Temporal Features Based on Monthly VIIRS NTL Data
HA decomposed a periodic function into a series of basic sinusoids, a set of sine or cosine terms having unique coefficients and intercept [28], [41]. The function f(t) represents a general cyclic model with harmonic components [42] f (t) = a 0 + 3 n=1 a n cos 2πnt where f t : radiance intensity for the VIIRS NTL at t month from model prediction.
t: the tth month from 2018 to 2020, a total of 36 months (t = 1, 2, 3, …, 35, and 36) n: temporal frequency of the harmonic component (n = 1, 2, and 3) L: number of months per year (L = 12) a 0 : coefficient for overall value for the VIIRS NTL radiance intensity value c 1 : coefficient for interannual change (slope) for the VIIRS NTL radiance intensity value a n , b n : coefficients for intra-annual change, intra-annal bimodal changes, and intra-annual trimodal change for the VIIRS NTL radiance intensity value, n = 1, 2, and 3, respectively. Higher order harmonic components have a higher frequency. We set the highest frequency at 3 to remove higher frequency variations resulting from noise, which could obtain the robust temporal features (a 0 , c 1 , a 1 , b 1 , a 2 , b 2 , a 3 , and b 3 named VI-IRS_b1, VIIRS_b2, …, and VIIRS_b8, respectively) of VIIRS NTL.
Some studies confirmed that the first-order and second-order harmonic terms of time series NDVI allow ecological interpretation [43], [44] as they represent the annual and intra-annual bimodal signals, respectively. Similarly, the first-order harmonic terms of time series VIIRS NTL, i.e., VIIRS_b3 and VIIRS_b4, describe the annual signal of land cover phenology (LCP) and comprehensive human activities (CHA). While the second-order harmonic terms, i.e., VIIRS_b5 and VIIRS_b6, were more likely to represent the LCP and CHA with bimodal (i.e., two peaks per year) phenological characteristics.

C. Selecting Spatial-Temporal Features
To select significant and noncollinear variables from all spatial-temporal features, the bivariate correlation analysis and variable inflation factor (VIF) were applied. Only statistically significant independent variables have explanatory power for the dependent variable. First, bivariate correlation analysis was employed to analyze the correlation of two variables. Then, the VIF method was used to escape the multicollinearity of independent variables. When multicollinearity exists, the variance of the independent variable inflates. It is recommended that if the VIF value of any independent variable is greater than 5, and then there is high collinearity with the other independent variables, which will lead to significant standard errors while performing regression [45]. Therefore, we adopted the VIF method to select noncollinear features. In each VIF iteration, the feature with the highest VIF value was removed until the VIF values of all spatial-temporal features were less than 5.

D. Developing STRCRF Model
The basic idea for intercalibrating LJ1-01 DN intensity with VIIRS radiance intensity is to learn the relationships between LJ1-01 DN intensity and the spatial-temporal features at 130 m resolution using ML algorithms. In this study, we developed the STRCRF model to quantify the relationships. As shown in Fig. 4, the STRCRF model comprises two RF models using spatial-temporal features. RF is a relatively modern and highly flexible ML approach that integrates multiple decision trees. It can evaluate the importance of each feature during the regression [32]. Overall, when the STRF model fitting and simulating (5) is over, the residuals between the LJ1-01-like and LJ1-01 NTL could be obtained. To correct the residuals and make them closer to zero, we further input the LJ1-01-like data from the STRF model into the STRCRF model and train the model (6) to make the LJ1-01-like closer to LJ1-01 NTL.
In the end, we regarded the coefficient determination (R 2 ), and RMSE of ten-fold CV as the accuracy indicators to validate the performance of the STRCRF model. Besides, the mean-residuals (MEAN RES ) of all pixels within the study area were employed to compare the effect of residuals correction where (LJ1 − 01) i and LJ1 − 01 − like) i are the values of the ith pixel in the LJ1-01-corrected and simulated LJ1-01-like NTL, respectively; (LJ1 − 01) is the mean values of the LJ1-01-corrected data; and k is the total number of pixels.

E. Evaluating of LJ1-01-Like NTL
We evaluated the LJ1-01-like NTL from three aspects, i.e., statistic analyses, spatial patterns, and profile analyses at GZMC. Statistic analyses included linear relationships between different NTLs, i.e., VIIRS, LJ1-01, and LJ1-01-like, at the pixel level. Meanwhile, we summed the pixel value of the ISP area and POP statistic data at the district level. The linear relationships between them and the sum of total light (SOTL) of VIIRS, LJ1-01, and LJ1-01-like NTL were built. On the other hand, we used the spatial error ratio (SER) to evaluate the spatial error of LJ1-01-like NTL. To further evaluate intuitively the simulated NTL of the STRCRF model, we adopted spatial pattern analyses and profile analyses. Then, we could find the difference and similarity of the above NTLs in spatial textures and trends, respectively

A. Descriptive Statistics of Spatial and Temporal Features
Before statistics of spatial and temporal features, it is essential to consider the correlation and collinearity of these features and to select significant and nonlinear information to avoid duplication of information between the seemly features. For this, the bivariate correlation analysis and VIF method were chosen to select features (i.e., VIIRS, VIIRS_b1-VIIRS_b8, DEM, Slope, NDVI, RD, CD, WD, GD, ISPD, and FD). A total of 18 features were evaluated based on the above methods. Through correlation analysis, all features have a significant Pearson correlation coefficient (CC) at the 0.01 level (see Table VI). Then, in each VIF iteration, the feature with the highest VIF was removed. Finally, three features, including VIIRS_b2, ISPD, and FD, were excluded, and 15 features with VIF <5 were then used to fit the model (see Table VII). Fig. 5 plots the histograms of all pixel-scale features used to fit the LJ1-01-like estimation models from 124 166 samples over GZMC. In general, the histogram distribution of LJ1-01 intensity is similar to which of VIIRS intensity and VIIRS_b5. Besides, the VIIRS_b4 shows negative distributions with the LJ1-01 intensity. The bivariate correlation analysis illustrates that all features were significantly related to LJ1-01 intensity, especially VIIRS (CC = 0.629, p < 0.01), VIIRS_b5 (CC = 0.623, p < 0.01), and VIIRS_b4 (CC = −0.575, p < 0.01, Table VI).   To test whether the proposed method is applicable for regions outside GZMC, Shenzhen, Dongguan, and Huizhou were validated by this method. The simulated imageries in these three cities in March 2019 were also validated with the observed imagery based on R 2 and RMSE (see Table IV). Owing to the hierarchical structure of Guangzhou, Shenzhen, Dongguan, and Huizhou in the GBA [36], the performance of the STRCRF model showed an apparent decreasing trend from Shenzhen [R 2 = 0.945, RMSE = 7298 (30.91% of the average DN value)] to Huizhou [R 2 = 0.918, RMSE = 4139 (46.20% of the average DN value)]. These results confirmed that the performance of the STRCRF model has close relationships with the differences  To evaluate the application potentiality of the simulated LJ1-01-like NTL, the correlation between SOTL of NTL (VIIRS, LJ1-01, and LJ1-01-like) and the sum of the POP statistics and ISP area were assessed at the district level. As shown in Fig. 7, LJ1-01-like and POP have a significant linear relationship with an R 2 value of 0.561 [see Fig. 7(c)], which is higher than that of VIIRS and POP (R 2 = 0.436) in Fig. 7(a). Additionally, the trend of NTL is in line with the ISP area in six districts, and LJ1-01-like NTL has a higher R 2 value of 0.795 [see Fig. 7(f)] than that of VIIRS (R 2 = 0.769) in Fig. 7(d). Moreover, the relationship between LJ1-01-like and statistics data (POP and ISP area) [see Fig. 7(c) and (f)] is close to that of LJ1-01 NTL [see Fig. 7(b) and (e)]. Therefore, the linear relationship between SOTL in LJ1-01-like and POP, ISP area can promote the application of LJ1-01-like NTL in urban populations and human footprints at a finer scale.

B. Model Fitting and Validation
2) Spatial Patterns: Fig. 8 shows the NTLs of different data sources and SER of LJ1-01-like in GZMC in March 2019. Compared with other NTL images, VIIRS NTL lost more details because of the limitation of coarse spatial resolution (500 m). Additionally, VIIRS NTL imagery was blooming in the suburban area [see Fig. 8(b)]. After intercalibration between VIIRS and LJ1-01 NTL, LJ1-01-like NTL was obtained. Compared  . The SER was also calculated to demonstrate the error distribution in space [see Fig. 8(d)]. The area with a high value of SER (in purple) only accounts for 1.89%, mainly distributed in the southeastern area, which covers high-density vegetation. Additionally, four typical areas were selected to compare the above three NTLs. We found that the LJ1-01 and LJ1-01-like NTL have apparent advantages compared with VIIRS NTL [see Fig. 8 Besides, the simulated results in three validation cities, i.e., Shenzhen, Dongguan, and Huizhou are displayed in Fig. 9. Overall, the STRCRF model still obtained high-quality LJ1-01-like NTL with detailed information. On the one hand, compared with VIIRS NTL, the overall spatial distribution of LJ1-01-like NTL was not changed; on the other hand, compared with LJ1-01 NTL, some similar details were built in LJ1-01-like NTL. These demonstrated the effectiveness of the STRCRF model.
3) Profile Analyses: The fluctuation of LJ1-01-like NTL (the gray part) agrees well with that of the LJ1-01 NTL (solid lines) from the two profiles across GZMC in Fig. 10. In GZMC, the ranges and distribution patterns of LJ1-01-like NTL were similar Fig. 11. Importance ranking of the STRF modeling features. The blue and orange histograms denote the importance of spatial and temporal features, respectively.
to that of LJ1-01 NTL, especially in the medium-intensity NTL area [see Fig. 10(b)], even though LJ1-01-like profiles showed slight overestimation in general. Meanwhile, parts of the LJ1-01like NTL profiles had a more serious overestimation within the urban center area [e.g., ID = 244 in Fig. 10(a)] compared with LJ1-01 NTL. This issue is mainly caused by the abnormally high values in the urban center area, while the surrounding pixels with abrupt high intensity. Under such a situation, the surrounding pixels might improve the center pixel's NTL intensity owing to the accumulated effect.

A. Comparison of Spatial-Temporal Features Importance
The spatial and temporal features are both input into the RF model to provide helpful information to enhance the consistency between LJ1-01 and VIIRS NTL. In the STRF regression process, the model obtained the importance scores by calculating the variance of each feature. In Fig. 11, the most important feature was VIIRS NTL intensity accounting for 30% nearly. In other words, the importance of other features was summed up to 70%, which means that spatial heterogeneity of land cover and temporal fluctuation of VIIRS NTL time series indeed affects the spatial intensity of LJ1-01 NTL. Moreover, the temporal features contributed more than the spatial features expected for VIIRS. The VIIRS_b5 (the intra-annual bimodal signal) and the VIIRS_b4 (the intra-annual signal) were the top factors among temporal features. This result affirmed that the intra-annual and intra-annual bimodal signals of NTL affect the spatial intensity of LJ1-01 to some degree. RD was much more important than the remaining spatial features, which means a close relationship between the road networks and NTL. Therefore, some research combined road vector data to correct the saturation issue of NTL [33].

B. Comparison With Related Studies
Due to the serious inconsistency of different kinds of NTL datasets hindering their application in the long term, there are many studies trying to enhance the consistency of NTLs in spatial and radiation resolution. In earlier studies, the intercalibration object was DMSP NTL [12], [22], [24]. Owing to sensor degradation and lack of on-board calibration, the DMSP NTL time series could not be integrated. These calibration methods corrected the DN values in the same spatial resolution (1000 m), thus achieving tangible results at the global and regional levels. With the development of NTL sensors, VIIRS NTL data with 500 m spatial resolution and alleviating saturated and blooming issues were developed [19]. Then, some studies have occurred to resolve the inconsistency between DMSP and VIIRS NTL. Based on simulated NTL, DMSP-like NTL with 1000 m spatial resolution [11], [27], [28], [29] or VIIRS-like NTL with 500 m spatial resolution [30], [31], [32] could be generated.
In this study, we developed an STRCRF model to intercalibrate the LJ1-01 and VIIRS NTL and simulated LJ1-01-like NTL. To be comparable with our study, only these studies intercalibrating of different NTLs were selected (see Table V). Considering the differences in fitting relationships, the linear and nonlinear regression models were chosen. The STRCRF model has CV-R 2 of 0.982 and CV-RMSE of 3469 at the pixel scale, and outperforms most previous models used for simulated 130 m LJ1-01-like NTL product, e.g., the LR model (CV-R 2 = 0.579, CV-RMSE = 15825; [29]) and the RF model (CV-R 2 = 0.765, CV-RMSE = 12713; [32]). Meanwhile, the application validation of LJ1-01-like from the R 2 _pop of the STRCRF model (0.561) is better than which of the other two models, i.e., the LR model (0.541) and RF model (0.549). The comparison results affirm that the STRCRF model performs better than most models in intercalibrating NTLs.

C. Limitation
In this study, we did not consider the differences in the overpass time of LJ1-01 and VIIRS NTL. However, the pixels of NTL images may change entirely because of this difference [46], which may lead to errors in training samples. To further improve intercalibration performance, more studies can focus on reducing the impact of the overpass time differences between LJ1-01 and VIIRS. Additionally, removing the abnormally high values of LJ1-01 NTL data based on the threshold method still has some drawbacks. We still found some abnormal high pixels in the urban center area. Therefore, more experiments can be operated to denoise the LJ1-01 NTL in future studies [47]. Moreover, VIIRS NTL seasonality and the characteristics of features are stable except for VIIRS radiance value with timing changes; those factors affect the time prediction performance of the STRCRF model (see Table VIII). In future studies, it is necessary to remove the seasonality of NTL and consider model features before building the intercalibration model.

VI. CONCLUSION
To generate high-quality LJ1-01-like NTL data, in this study, we proposed a novel intercalibration model for VIIRS and LJ1-01 NTL in March 2019 by considering the GZMC as the study area, Shenzhen, Dongguan, and Huizhou as the validation cities to confirm the spatial simulation ability of the model. First, the related spatial-temporal features were generated based on spatial analysis and harmonic analysis. Second, the significant and noncollinear spatial-temporal features were selected by bivariate correlation analysis and the VIF method. Third, by adopting the Consistency between the LJ1-01-like and LJ1-01 NTL was validated with evaluations on statistic analyses, spatial patterns, and profile analyses. As results indicated that the temporal features contribute more information than spatial features during the intercalibration process. Compared with most models presented in the previous studies, the STRCRF outperformed model fitting, prediction, and application. The newly simulated LJ1-01-like NTL data with more spatial details and fewer blooming issues than VIIRS NTL. More urban high-quality NTL information at the finer scale is, thus, obtained. These results illustrate that the STRCRF model is effective in intercalibrating the spatial inconsistency between the LJ1-01 and VIIRS NTL, and LJ1-01-like NTL data can be helpful in NTL-based studies such as mapping human activities and monitoring the urban environment.

ACKNOWLEDGMENT
We thank Wuhan University for providing Luojia1-01 nighttime light data. We also thank the anonymous reviewers for helping us to improve the manuscript.