A robust gap-filling method for predicting missing observations in daily Black Marble nighttime light data

ABSTRACT

Nighttime light (NTL) remote sensing data plays a crucial role in comprehending changes in human activities. The availability of the daily lunar BRDF-corrected Black Marble NTL product (VNP46A2) enables the use of NTL data to detect and assess the impact of short-term emergencies. However, daily NTL data often experience missing values due to cloud cover and low-quality signals. To address this issue, many studies utilize monthly or annual time-composite NTL products, which restrict the timeliness and potential application scenarios of NTL data usage. Therefore, it is necessary to generate the gap-filled daily NTL product. This study presented a novel NTL gap-filling method comprising rough reconstruction based on spatiotemporal weighting and refined gap-filling using a Bidirectional Long Short-Term Memory (Bi-LSTM) model. We evaluate the accuracy of the proposed method using the “remove-reconstruct-compare” approach, which randomly removes some original data from the complete image, fills the gaps with the proposed gap-filling method, and compares the reconstructed NTL data with the original observations in Beijing, Shanghai, Xi’an and New York. The results reveal that when the rate of missing values in Beijing is 40% and 50%, the proposed gap-filling method achieves accuracy with mean coefficient of determination (R²) values of 0.834 and 0.841, accompanied by corresponding root mean square (RMSE) values of 7.793 and 7.171, respectively. Furthermore, the gap-filling accuracy was evaluated quantitatively, and our proposed gap-filling method performed better than the Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM). Our proposed gap-filling method had R² values of 0.685, 0.781, 0.720 and 0.642, which were higher than those for STARFM (0.430, 0.662, 0.221 and 0.345). The RMSE values of our gap-filling method were 9.628, 12.083, 10.963 and 19.882 for the four sites, while those of STARFM were 12.953, 14.872, 18.280 and 26.990, respectively. The temporal and spatial analysis results demonstrate that this model is robust, capturing city boundaries and NTL high-brightness hotspots with high accuracy and stability. The gap-filling model proposed in this study provides a new technique for expanding the potential applications and reliability of NASA’s daily Black Marble product (VNP46A2) in remote sensing.

KEYWORDS:

• Nighttime light
• gap-filling
• deep learning
• NASA black marble product

1. Introduction

The nighttime light (NTL) data from satellite sensors provide a spatial and temporal record for observing and monitoring human activities. For instance, it offers a new perspective for understanding urban expansion patterns and investigating the influence of human activities on urban natural environments (Levin et al. 2020 leShi et al. 2014). Since the 1970s, a series of satellite sensors have been used for NTL observations, such as the Defense Meteorological Satellite Program-Operational Linescan System (DMSP-OLS) (Elvidge et al. 1999), the Visible Infrared Imaging Radiometer Suite (VIIRS) sensor on the Suomi National Polar-orbiting Partnership (NPP) Satellite (NPP-VIIRS) (Hillger et al. 2013), the Luojia1-01satellite (Wang et al. 2021), the high spatial resolution commercial satellite Jilin1-03B (Cheng et al. 2020), and the commercial satellite EROS-B (Levin et al. 2014). With rapid global urbanization, there is a sharp increase in the use of NTL data to study environmental changes caused by human activities (Tan, Zhou, and Bai 2017; Zheng, Chen Teo, and Pin Koh 2021) and socioeconomic development (Xie and Weng 2015). The release of daily NTL products enables studies to examine events with higher temporal resolution and better insights into human dynamics related to holidays (Ma 2018), conflicts (Zheng et al. 2022), and population movements (Mu et al. 2021, 2022; Stokes and Román).

As the successor to DMSP-OLS NTL remote sensing data, the NTL product obtained from NPP-VIIRS has substantially improved spatial resolution and sensor calibration. Recently high-quality daily VIIRS products have been released and made freely accessible to the public in NASA’s Black Marble NTL product suite (Román et al. 2018). However, NTL data obtained from the VIIRS sensor can be affected by weather, transient light sources, background noise, and atmospheric scattering (Elvidge et al. 2017; Miller and Turner 2009), leading to significant data gaps in certain areas.

Large area gaps dramatically limit the applications of NTL data, such as the assessment of emergency events that require high-quality NTL coverage at a high temporal frequency (Román et al. 2019; Xu and Qiang 2021; Zhao et al. 2018). The potential to monitor abrupt short-term changes can be significantly improved by high-frequency temporal NTL observations. Therefore, the imputation of missing data has important practical implications.

The gap-filling methods of remote sensing data can improve data quality to solve this problem. These methods are divided into three categories (Yao et al. 2021). The first set of methods uses spatial information to fill the missing data and assumes that data are spatially autocorrelated. The missing values are filled by finding the image pixels similar to the missing pixel within the effective spatial window, such as the spline function method (Ding and Li 2013), the regression tree analysis (Wei et al. 2019), and the kriging method (Gotway et al. 1996). These methods do not require additional multi-source information and are simple to implement, but they have many drawbacks. For example, these methods lack better spatial detail and higher gap-filling accuracy, particularly when there are higher rates of missing data in the image. Therefore, these methods are only suitable for cases where the image has a small amount of missing data.

The second set of methods uses temporal information to fill the missing data. The temporal filling method assumes that the feature will change periodically or non-periodically over time. The successive changes constitute the temporal correlation of the feature, and the shorter the interval, the closer the feature properties. The missing values are filled by finding the image pixels similar to the missing pixel within an effective temporal window, such as harmonic analysis (Katznelson 2004), multi-temporal robust regression, singular spectrum analysis (Hassani 2021), and SG filtering method (Jahani et al. 2018). These interpolation methods only learn the temporal features of the time series in the temporal dimension, and smoothing the time series is inevitable, leading to information loss related to abrupt changes. Therefore, the methods described above only apply to images without continuous missing data in the time series.

The third set of methods uses spatiotemporal information to fill the missing data (Militino et al.2019 ; Weiss et al. 2014). The spatiotemporal filling method assumes that the spatial correlation of the feature should not change over time if the feature type remains unchanged under natural conditions. It employs this assumption to find similar pixels within an effective spatiotemporal window to fill the missing values. Compared to spatial or temporal filling methods, the spatiotemporal filling method can utilize more information and achieve higher accuracy, particularly in areas with many missing values. However, the high spatial and temporal heterogeneity of NTL data still dramatically affects the gap-filling accuracy.

In summary, traditional gap-filling methods that use spatial and temporal domains are no longer suitable when numerous missing values exist in the satellite remote sensing data. These methods result in substantial information loss and reduced accuracy, which cannot meet the demand of practical applications. However, with the increase in satellite data and computational power in the last decade, deep learning has surpassed traditional machine learning models in estimating and predicting various remote sensing parameters (Zhu et al. 2017). Traditional spatiotemporal gap-filling methods often fail to adaptively transform the model parameters under complex and different scenarios, thus hindering the improvement of the gap-filling results. Deep learning-based data reconstruction models can maintain high gap-filling accuracy with high robustness under different complex scenes based on the model learning enough high-quality samples (Arslan and Sekertekin 2019; Yang et al. 2018). Moreover, it can capture the spatial details and temporal patterns of NTL data. The long short-term memory (LSTM) network is a type of recurrent neural network, and Bidirectional LSTM (Bi-LSTM) is a variant of the LSTM network, which consists of both forward LSTM and backward LSTM (Cui et al. 2022). It introduces an internal mechanism called “gate” to regulate the flow of information. These gating structures can learn which data in the sequence are important information to keep and which to delete. By doing so, it can pass relevant information along long chains of sequences to perform predictions. Compared with other deep learning algorithms, the Bi-LSTM network can better learn short-term features of the entire time series and prevent abrupt changes from being easily smoothed (Imhoff et al. 1997). Thus, by coupling the spatiotemporal gap-filling method and the deep learning process model, the gap-filling ability is further optimized based on spatiotemporal weighted gap-filling results, leading to more precise prediction results.

This study proposes a novel algorithm that couples the spatiotemporal weighted gap-filling and Bi-LSTM model for NTL data reconstruction. Specifically, the novel algorithm includes two main steps. Firstly, the proposed algorithm uses spatiotemporal weighted gap-filling to extract the spatial details and overall data trends for coarse NTL data reconstruction. Then it uses the Bi-LSTM model to update and complete the refined NTL data reconstruction by learning short-term features of the data.

2. Study area and data

2.1. Study area

The study was conducted in Beijing, the northern region of China (Figure 1), and the center of the research area is located at 116.51°E, 40.01°N. The research area includes most of the urban built-up areas in Beijing. The terrain decreases in turn from northwest to southeast. The study area is surrounded by mountains on three sides, including the west, north, and northeast. The southeast is a plain that slopes gently to the Bohai Sea. The built-up area of the city is mainly concentrated in the plain area. Artificial light sources from urban areas mainly produce the NTL data from satellite sensors.

Figure 1. Image of black marble VNP46A2 in the study area of Beijing. Beijing is the region for modelling research, and three cities including Xi’an, Shanghai, and New York are selected as verification test areas.

As China’s economic and cultural center, Beijing has many commercial and residential areas distributed throughout the city. The NTL changes, often accompanied by economic and cultural activities in these areas, are reflected in the dynamic trends and patterns of economic development in different functional areas of the city. Due to large cities’ highly complex spatial heterogeneity, the NTL data fluctuate significantly within the city, making data gap-filling challenging. Therefore, choosing this area as the research area enables evaluating the accuracy and reliability of the NTL gap-filling model proposed in this paper. In addition, to validate the robustness and generalization of the proposed gap-filling method, we selected Xi’an, Shanghai, and New York as the verification test areas.

2.2. Data and preprocessing

The dataset used in this research comprises of daily measurements of moonlight and atmospheric-corrected NTL, obtained from NASA’s VIIRS Black Marble product (VNP46A2, Collection V001). The VNP46A2 product consists of seven layers, which provide detailed information on various aspects of the BRDF-corrected NTL (500 m), including lunar irradiance, mandatory quality flag, latest high-quality retrieval (number of days), snow flag, and cloud mask flag (2022). A series of corrections including lunar BRDF, clouds, terrain, atmosphere, airglow and stray light corrections, improve the quality of the VNP46A2 product. These minimize the distortion caused by background noise in our subsequent analysis. We downloaded h29v04 and h29v05 tiles overlapping with the research area in 2017 and conducted statistics on the missing NTL data in the study area. As shown in Figure 2, the results indicate that the missing data in the study area is quite severe, with an average missing rate of 42.89%. Therefore, filling the missing data is necessary before it can be used for further analysis. As a data pre-processing process, we sift the BRDF-corrected NTL layer through a mandatory quality flag and cloud mask flag. Then, we remove the outliers to get a cloud-free collection of high-quality imagery for this research.

Figure 2. Overview of the missing data fraction of VNP46A2 data in 2017 in the study area.

3. Method

The detailed process of the NTL gap-filling method is shown in Figure 3, and the NTL gap-filling model is divided into two steps: rough NTL gap-filling method based on spatiotemporal features and refined NTL gap-filling method based on the Bi-LSTM model. In the first step, we use the spatiotemporal information of adjacent images and weight indicators such as distance and variance to identify similar pixels for spatiotemporal weighted gap-filling. The spatiotemporal weighted gap-filling method proposed in this paper differs from the spatiotemporal approach proposed by Weiss (Weiss et al. 2014) in two main aspects: weight indices selection and input data choice. Moreover, we enhance the method proposed by Weiss by incorporating a standard deviation index to reduce errors resulting from spatial heterogeneity variations between the target image and its neighboring counterparts. By considering these distinctions, our method becomes more universally applicable and can efficiently fill gaps in other remote sensing data. The second step establishes the Bi-LSTM model, which learns the samples’ long-term trend and short-term characteristics in the NTL time series. With the Bi-LSTM model, we can iteratively fill and update the spatiotemporal weighted gap-filled NTL data to complete the refined NTL data.

Figure 3. The overall flowchart of the NTL gap-filling method.

3.1 Rough NTL spatiotemporal gap-filling method

3.1.1 Define the spatial and temporal window of image collection

The detailed process of the rough spatiotemporal gap-filling method is shown in Figure 4. First, it is necessary to define a subset of NTL images and determine an appropriate spatial and temporal window to find valid pixel pairs for filling the missing pixel values. A subset of NTL images that need to be filled was selected as the data source of the target image. In addition, larger and longer spatial and temporal window sizes increase the computational load without necessarily improving the accuracy. Therefore, selecting an appropriate spatial and temporal window that balances computational efficiency and accuracy is crucial. In our study, a temporal window size of m days for the target image which includes the previous m/2 days and the following m/2 days, and a spatial window size of n×n pixels centered on the target pixel were model parameters in the gap-filling method that needed to be set. After some experiments and reference to previous research (Gerber et al. 2018; Sun et al. 2017), we chose the temporal window size of 9 days and the spatial window size of 15 × 15 pixels as the parameters in the rough gap-filling method. After this selection and process, the valid NTL pixels in the spatial and temporal windows from the image collection were extracted. These NTL observations were the prediction set for the subsequent gap-filling method.

Figure 4. The framework for rough NTL spatiotemporal gap-filling method.

3.1.2 Calculating weight indicators and predicting the missing pixel value

The pixels in the prediction set were divided into three categories: (a) NTL pixels from different images that are at the same location as the target pixel; (b) NTL pixels from the same image but situated at different locations as the target pixel: (c) other NTL pixels. This prediction set includes different temporal and spatial information. The NTL predictions of the target missing pixel can be predicted using Eq (1)(1) $\mathrm{N}\mathrm{T}{\mathrm{L}}_{P_1}=\mathrm{N}\mathrm{T}{\mathrm{L}}_{{\mathrm{x}}_0,{\mathrm{y}}_0,{\mathrm{i}}_{\mathrm{t}}}+\mathrm{N}\mathrm{T}{\mathrm{L}}_{{\mathrm{x}}_{\mathrm{j}},{\mathrm{y}}_{\mathrm{j}},{\mathrm{i}}_0}-\mathrm{N}\mathrm{T}{\mathrm{L}}_{{\mathrm{x}}_{\mathrm{j}},{\mathrm{y}}_{\mathrm{j}},{\mathrm{i}}_{\mathrm{t}}}$(1) ,which uses two kinds of information.

(1) $\mathrm{N}\mathrm{T}{\mathrm{L}}_{P_1}=\mathrm{N}\mathrm{T}{\mathrm{L}}_{{\mathrm{x}}_0,{\mathrm{y}}_0,{\mathrm{i}}_{\mathrm{t}}}+\mathrm{N}\mathrm{T}{\mathrm{L}}_{{\mathrm{x}}_{\mathrm{j}},{\mathrm{y}}_{\mathrm{j}},{\mathrm{i}}_0}-\mathrm{N}\mathrm{T}{\mathrm{L}}_{{\mathrm{x}}_{\mathrm{j}},{\mathrm{y}}_{\mathrm{j}},{\mathrm{i}}_{\mathrm{t}}}$(1)

where NTL_p1 is the NTL prediction of the target missing pixel; (x₀,y₀) is the location of the missing pixel in the target image; i₀ is the target image; NTLx₀,y₀,i_t is the NTL observation at the location (x₀,y₀) in the image i_t, which is an image in the subset collection. The subset collection belongs to category (a), NTLx_j,y_j,i₀ is the NTL observation at the location (x_j,y_j) in the target image; (x_j,y_j) is a location in the spatial window; NTLx_j,y_j,i₀ belongs to category (b); NTLx_j,y_j,i_t is the NTL observation at the location (x_j,y_j) in image i_t that is adjacent to the target image and belongs to category (c). The spatial explanation of the formula variable for Eq.(1) can be seen from Figure 5.

Figure 5. Spatial diagram for interpreting formula variables.

This spatiotemporal gap-filling method referred to the previous research (Sun et al.) and was based on the hypothesis that the NTL difference between two days of a pixel will be similar to that of a nearby pixel that can be expressed as Eq (2)(2) $\mathrm{N}\mathrm{T}{\mathrm{L}}_{P_1}-\mathrm{N}\mathrm{T}{\mathrm{L}}_{{\mathrm{x}}_0,{\mathrm{y}}_0,{\mathrm{i}}_{\mathrm{t}}}=\mathrm{N}\mathrm{T}{\mathrm{L}}_{{\mathrm{x}}_{\mathrm{j}},{\mathrm{y}}_{\mathrm{j}},{\mathrm{i}}_0}-\mathrm{N}\mathrm{T}{\mathrm{L}}_{{\mathrm{x}}_{\mathrm{j}},{\mathrm{y}}_{\mathrm{j}},{\mathrm{i}}_{\mathrm{t}}}$(2) :

(2) $\mathrm{N}\mathrm{T}{\mathrm{L}}_{P_1}-\mathrm{N}\mathrm{T}{\mathrm{L}}_{{\mathrm{x}}_0,{\mathrm{y}}_0,{\mathrm{i}}_{\mathrm{t}}}=\mathrm{N}\mathrm{T}{\mathrm{L}}_{{\mathrm{x}}_{\mathrm{j}},{\mathrm{y}}_{\mathrm{j}},{\mathrm{i}}_0}-\mathrm{N}\mathrm{T}{\mathrm{L}}_{{\mathrm{x}}_{\mathrm{j}},{\mathrm{y}}_{\mathrm{j}},{\mathrm{i}}_{\mathrm{t}}}$(2)

It should be noted that when making predictions, the variables on the right side of Eq (1)(1) $\mathrm{N}\mathrm{T}{\mathrm{L}}_{P_1}=\mathrm{N}\mathrm{T}{\mathrm{L}}_{{\mathrm{x}}_0,{\mathrm{y}}_0,{\mathrm{i}}_{\mathrm{t}}}+\mathrm{N}\mathrm{T}{\mathrm{L}}_{{\mathrm{x}}_{\mathrm{j}},{\mathrm{y}}_{\mathrm{j}},{\mathrm{i}}_0}-\mathrm{N}\mathrm{T}{\mathrm{L}}_{{\mathrm{x}}_{\mathrm{j}},{\mathrm{y}}_{\mathrm{j}},{\mathrm{i}}_{\mathrm{t}}}$(1) must exist at the same time. The core assumption, as stated in Eq (2)(2) $\mathrm{N}\mathrm{T}{\mathrm{L}}_{P_1}-\mathrm{N}\mathrm{T}{\mathrm{L}}_{{\mathrm{x}}_0,{\mathrm{y}}_0,{\mathrm{i}}_{\mathrm{t}}}=\mathrm{N}\mathrm{T}{\mathrm{L}}_{{\mathrm{x}}_{\mathrm{j}},{\mathrm{y}}_{\mathrm{j}},{\mathrm{i}}_0}-\mathrm{N}\mathrm{T}{\mathrm{L}}_{{\mathrm{x}}_{\mathrm{j}},{\mathrm{y}}_{\mathrm{j}},{\mathrm{i}}_{\mathrm{t}}}$(2) , assumes that the NTL difference between two consecutive days of a pixel is similar to that of an adjacent pixel. The essence of the following rough NTL spatiotemporal gap-filling method is to treat a single missing pixel as a unit. The goal is to identify as many “m” pixel pairs in Eq (2)(2) $\mathrm{N}\mathrm{T}{\mathrm{L}}_{P_1}-\mathrm{N}\mathrm{T}{\mathrm{L}}_{{\mathrm{x}}_0,{\mathrm{y}}_0,{\mathrm{i}}_{\mathrm{t}}}=\mathrm{N}\mathrm{T}{\mathrm{L}}_{{\mathrm{x}}_{\mathrm{j}},{\mathrm{y}}_{\mathrm{j}},{\mathrm{i}}_0}-\mathrm{N}\mathrm{T}{\mathrm{L}}_{{\mathrm{x}}_{\mathrm{j}},{\mathrm{y}}_{\mathrm{j}},{\mathrm{i}}_{\mathrm{t}}}$(2) as possible with missing pixels in adjacent spatiotemporal windows of the image. Each pixel pair will receive a predicted value after the formula calculation. The weight coefficients for various pixel pairs are determined using Eq (3)(3) $\mathrm{N}\mathrm{T}{\mathrm{L}}_{P_2}=\frac{\sum _{n=1}^m\mathrm{N}\mathrm{T}{\mathrm{L}}_{P_1,n}\times w_n}{\sum _{n=1}^m{w}_n}$(3) . The predicted value of the missing pixel is then calculated based on the weighted prediction of these pixel pairs.

(3) $\mathrm{N}\mathrm{T}{\mathrm{L}}_{P_2}=\frac{\sum _{n=1}^m\mathrm{N}\mathrm{T}{\mathrm{L}}_{P_1,n}\times w_n}{\sum _{n=1}^m{w}_n}$(3)

The NTLp₂ is the final NTL prediction for the missing pixel by bringing in spatiotemporal information from adjacent pixels. The m is the total number of valid pixels from the 15 × 15 pixels window in each target image; NTLp₁, n is the mth predicted NTL, and w_n that can be expressed as Eq (4)(4) $w_n=\frac{1}{\mathrm{D}{\mathrm{I}}_n\times \mathrm{S}{\mathrm{I}}_n\times \mathrm{S}\mathrm{D}{\mathrm{I}}_n}$(4) is the weight of the nth predicted NTL; the weight w_n was calculated using three indices, including the distance index, the similarity index, and the standard deviation index.

(4) $w_n=\frac{1}{\mathrm{D}{\mathrm{I}}_n\times \mathrm{S}{\mathrm{I}}_n\times \mathrm{S}\mathrm{D}{\mathrm{I}}_n}$(4)

The location of the target pixel in image t is denoted by (x_0,t, y_0,t), while the kth predicted NTL pixel location in image t is represented by (x_k,t, y_k,t) in the prediction set. The distance index assumes that data are spatially autocorrelated. The closer the distance between the NTL observations, the more similar the NTL observations. Hence, the distance index (DI) was calculated using Eq.(5):

(5) $\mathrm{D}{\mathrm{I}}_n=\sqrt{{(x_{0,t}-x_{k,t})}^2+{(y_{0,t}-y_{k,t})}^2}$(5)

The similarity index is used to filter out the adjacent pixel with the most similar situation to the center pixel that is to be filled. The smaller the difference between pixel pairs, the more similar they are. The minor difference represents the closest conditions of two nearby pixels, which should be set as the most considerable weight in kth weight. The similarity index (SI) was calculated using Eq (6)(6) $\mathrm{S}{\mathrm{I}}_n=\left|\mathrm{N}\mathrm{T}{\mathrm{L}}_{{\mathrm{x}}_0,{\mathrm{y}}_0,{\mathrm{i}}_{\mathrm{t}}}-\mathrm{N}\mathrm{T}{\mathrm{L}}_{{\mathrm{x}}_n,{\mathrm{y}}_n,{\mathrm{i}}_{\mathrm{t}}}\right|+1$(6) :

(6) $\mathrm{S}{\mathrm{I}}_n=\left|\mathrm{N}\mathrm{T}{\mathrm{L}}_{{\mathrm{x}}_0,{\mathrm{y}}_0,{\mathrm{i}}_{\mathrm{t}}}-\mathrm{N}\mathrm{T}{\mathrm{L}}_{{\mathrm{x}}_n,{\mathrm{y}}_n,{\mathrm{i}}_{\mathrm{t}}}\right|+1$(6)

The standard deviation index was used to measure the spatial standard deviation of NTL observations between the target image and the image i_t. The standard deviation of NTL observations between the target image and a given image was slight, and the difference in spatial heterogeneity between the target image and this image was reduced, which ensures the hypothesis of this spatiotemporal gap-filling method that the NTL difference between two images of a target pixel will be similar to that of a nearby pixel. The standard deviation index (SDI) was calculated using Eq (7)(7) $\mathrm{S}\mathrm{D}{\mathrm{I}}_n=\sqrt{\frac{1}{n}\times \sum _{k=1}^n{\left(\mathrm{N}\mathrm{T}{{\mathrm{L}}_{\mathrm{i}}}_{,r}-\overline{\mathrm{N}\mathrm{T}{\mathrm{L}}_{\mathrm{i}}}\right)}^2}$(7) :

(7) $\mathrm{S}\mathrm{D}{\mathrm{I}}_n=\sqrt{\frac{1}{n}\times \sum _{k=1}^n{\left(\mathrm{N}\mathrm{T}{{\mathrm{L}}_{\mathrm{i}}}_{,r}-\overline{\mathrm{N}\mathrm{T}{\mathrm{L}}_{\mathrm{i}}}\right)}^2}$(7)

where n is the total number of valid pixels from the spatial window in the target image, and NTL_i,r is the NTL difference between the target image and image r at the rth pixel and$\overline{\mathrm{N}\mathrm{T}{\mathrm{L}}_{\mathrm{i}}}$is the average NTL difference between the target image and image i.

3.1.3 Fill the remaining gaps

After the above spatiotemporal filling steps, there were still some gaps in the target image. Because the target image may be empty and the gap size was large, no valid pixel was found within the spatial window. To solve this problem, we determine a 3 × 3 spatial window and a temporal window of two adjacent days of the target image with the missing pixel as the center to select a valid pixel. The remaining missing pixels are filled by taking the average value of the collected valid pixels. After the above steps, 80% missing areas of the image have been filled. However, when there are more than 10 consecutive days of missing pixels in the area, some pixels are still not filled, and these pixels will be filled in the subsequent refined gap-filling model.

3.2 Refined gap-filling method based on the Bi-LSTM model

The schematic overview of the refined gap-filling method based on the Bi-LSTM model is illustrated in Figure 6. Firstly, we generate training samples by setting a sliding temporal window and adding temporal prior information features to the samples. Then we train the Bi-LSTM model to learn the temporal characteristics of the NTL time series. Finally, we fill the remaining gaps in NTL images with the Bi-LSTM model and iteratively update the gap-filled NTL values to complete the refined gap-filling process.

Figure 6. The framework for the refined gap-filling method based on the Bi-LSTM model.

3.2.1 Principle of the LSTM model

A recurrent neural network (RNN) is a neural network for processing sequential data. Compared to general neural networks, it can process data with sequence changes. Long short-term memory (LSTM) is a type of RNN developed to solve the problems of gradient disappearance and gradient explosion problems while training long sequences. Unlike the RNN network, LSTM can perform better on longer sequences (Hochester et al. 1997). Figure 7a shows that an LSTM cell consists of the input, output, and forget gates. It filters useful information and captures long-term information in time series.

Figure 7. External and internal structures of the Bi-LSTM neural network. (a) LSTM cell structure; (b) the whole structure of the Bi-LSTM neural network.

The previous cell state C_t-1 information is filtered by the forget gate using the LSTM cell state h_t-1 and x_t and the sigmoid function to output a value between 0 and 1. The current cell state C_t is calculated using the tanh activation value of the input gate and previous information of the forget gate filter. At last, the output of the cell h_t is computed using the sigmoid activation value of the output gate and current cell state.

The Bi-LSTM is another variant of the LSTM model containing two separate LSTM layers: the forward and backward. Figure 7b shows that the Bi-LSTM model can learn both the previous and future information of the time series at each time step, which can use the information in the past and future to capture the features of data and produce more accurate prediction results. It is also worth noting that the Bi-LSTM model requires a relatively complete time series to learn data features before predicting the missing value. The origin data are roughly reconstructed using the spatiotemporal gap-filling method, and most gaps in the image have been filled. Pre-filled NTL images can meet the Bi-LSTM model’s requirements for data integrity.

3.2.2 Construction and training of models

Sample selection and model training mainly include three parts. Firstly, to generate the training sample, we use date information, including day, week, and month which is added to the model as prior information, since these time-related patterns vary widely across the city due to public holidays and weekends. The addition of temporal-prior information allows the model to better capture the trend and pattern of the NTL time series, which is equivalent to a kind of time series decomposition.

Secondly, we calculate the annual average value of pre-filled NTL data and use the natural segmentation method to stratify the NTL observations into five categories using the natural breaks method (Chen et al. 2013), which indicates the different impacts of anthropogenic change disturbance at different levels. In each stratum, the sample points are randomly selected according to the principle of uniform sampling, and a total of 11,837 sample points are selected for model training.

We extract subsequence from the original NTL time series as training and testing samples. The sliding window algorithm anchors the left side in the first time index of the NTL time series and then keeps moving until the right side of the window reaches the end of the NTL time series (Zhang et al. 2019). We filter out samples with null values during the moving and generate the training and testing sample. Different sample NTL time series lengths impact the model’s training time and accuracy. The sliding window length is surveyed to determine the optimal length for the top-performing model. The samples are randomly divided into three parts (70% of them are used for training, 20% for validation, and the rest 10% for test). We use the validation set to select our best model. The final step is to train the model using training samples; the Bi-LSTM model is designed as one input layer, three Bi-LSTM layers with 50 Neurons, two dropout layers, and one dense layer. Furthermore, we use the Adam optimizer, with an initial learning rate of 0.0001, a batch size of 5096, and max epochs of 48.

3.2.3 Gapless image filling and updating

In some places where large areas are continuously missing in spatial and temporal in our study area, there are still 80% missing pixels in NTL images after spatiotemporal weighted gap-filling. So, the Bi-LSTM model is applied to fill the remaining gap in NTL data into a gapless image. To ensure whether the trained Bi-LSTM filling model can effectively fill the continuous missing values, it is crucial to perform experimental verification of the gap-filled results of the Bi-LSTM model under different missing days.

In addition, the Bi-LSTM model is used for the refined gap-filling to iterate the spatiotemporal weighted result. In the sequences, the NTL data of the missing value is predicted at the supervised learning format using existing 25-day observations. When the refined gap-filling NTL data of missing pixels are predicted, the input data of the LSTM gap-filling model are updated, and the rough NTL spatiotemporal gap-filling original result is displaced with the Bi-LSTM model prediction value. If there are n missing values in NTL spatiotemporal weighted sequence, this update process needs to be conducted n times to end the final NTL gap-filling. After a series of fills and updates, we had a gapless 2017 Beijing city daily NTL image collection.

3.3. Accuracy assessment

We conducted a comparative analysis method to assess the performance of our proposed algorithm at different missing rates. This method involves randomly removing some original data from the complete image, filling the missing data with the proposed gap-filling method, and comparing the original observed NTL values of the missing pixels with the filled NTL values. These artificial gaps were then filled with the method proposed in this study. Therefore, we select a day with relatively complete data as the evaluation image. The NTL data from the day of year (DOY) 325 of 2017 in the study area were removed for this comparative analysis to achieve different missing rates from 10% to 50% (10% intervals). The root mean square error (RMSE), mean absolute error (MAE), Bias, and coefficient of determination (R²) were calculated to evaluate the accuracy.

4 Results

4.1 Optimization of Bi-LSTM gap-filling model parameters

The selection of the optimal temporal length in the NTL time series is crucial for the performance of the Bi-LSTM model. During the training of the Bi-LSTM model, the temporal length of the NTL time series can impact both the training time and accuracy. A sliding window length test was conducted in study area of Beijing to determine the optimal temporal length by evaluating the training data with 10 different temporal lengths, ranging from 5 to 50 days with 5-day intervals.

During the Bi-LSTM model training stage, we only utilized training samples from study area in Beijing. In the process of optimizing model parameters, only data from study area in Beijing were involved in the analysis and presentation of results. The results in Figure 8 illustrate the Bi-LSTM model’s performance under different temporal lengths in study area of Beijing. Notably, the performance of R² across different temporal lengths for the training, validation, and test datasets were all above 0.75, indicating high accuracy. Furthermore, the RMSE for the same datasets was consistently lower than 12.5, indicating strong predictive capabilities regardless of the temporal length. The curves show that a longer temporal length leads to a lower RMSE for the training, validation, and test datasets. However, when generating time series training samples, for data quality, it is necessary to meet the condition that a single sample cannot contain missing values to filter samples. A temporal length that is too long may result in the inability to produce uniformly distributed samples in terms of space and quantity in regions where the original data is missing, leading to reduced generalization and robustness of the gap-filling model. After balancing the trade-off between accuracy and time consumption, a temporal length of 25 days was chosen as the optimal length.

Figure 8. Evaluation of temporal length in Bi-LSTM model in study area of Beijing. The upper and bottom curves represent R² and RMSE, respectively.

4.2 Accuracy of the gap-filling method

To evaluate the accuracy of our gap-filling method for NTL images, we first visually compared multiple pairs of gap-filled and actual NTL images from various months. We analyzed the trends and verified the continuity of typical landmarks in urban functional areas. Then, we quantitatively evaluated the gap-filled data using scatterplots and regression analysis based on all available precise NTL observations.

4.2.1 Visual comparison

The most straightforward method to determine the precision of the gap-filled images is to compare them with actual NTL images obtained at a similar time. To evaluate the effectiveness of the gap-filling approach regarding visual quality, a visual comparison was performed between images acquired on weekend, holiday and weekday before and after gap-filling in our study area. We selected four sets of before and after gap-filling comparison images in four cities, Beijing, Shanghai, Xi’an and New York, to assess the visual effect of our proposed gap-filling model. As shown in Figure 9, the results demonstrated that the proposed gap-filling method successfully filled many missing pixels, leading to gapless NTL images that retain their spatial patterns and details.

Figure 9. Comparison of the gap-filling effect before and after NTL data reconstruction in the research area. The blank in the original image shows the missing areas, and 57, 121, 212, 305 represents the day of year (DOY).

To emphasize the significance of using the proposed data, which differs from the Gap Filled BRDF corrected DNB NTL data, it is essential to highlight the advantages of the Bi-LSTM technique. This can be achieved by conducting a comprehensive comparison with various other gap-filling methods and existing products to ensure a robust and credible validation process. We compared the Gap Filled BRDF corrected DNB NTL with gap-filled data obtained using our proposed method, respectively, in time and space. We selected three locations in Shanghai area: airport, commercial area, and residential area to analyze the temporal sequence of the two types of gap-filled data. As we can see from Figure 10 that the Gap Filled BRDF corrected DNB NTL product lost crucial temporal information, which was a key feature for evaluating urban development status, urban function and urban spatial structure. The gap-filled data using our proposed method effectively captured the regional characteristics of the NTL time series in three regions. So it is necessary to fill the missing data based on the BRDF Corrected DNB NTL data.

Figure 10. A comparison of our result with existing product gap filled BRDF corrected DNB NTL is analyzed spatially and temporally in three typical urban functional areas of Shanghai: the residential area, the commercial area, and the airport. The blue curve represents our gap-filling results, and the red curve represents gap filled BRDF corrected DNB values.

To demonstrate the consistency and trend of the gap-filled NTL time series, it is crucial to examine and verify if the trend of the gap-filled NTL time series correctly reflects reality to meet the demands of further applications. Therefore, we selected representative locations in different urban functional areas to visualize the NTL time series, including the Beijing Capital International Airport, the Beijing Business District in Chaoyang, the Olympic Park, and the Tiantongyuan Community largest residential community in Asia, with a population of 400,000 to 500,000.

These locations were selected as they all feature frequent commercial and cultural activities, and the NTL data in these areas exhibit complex heterogeneity. By using the model to fill the complete NTL time series in these areas, the actual progression of the light, trend characteristics, and changing patterns can be demonstrated, proving the rationality and robustness of the gap-filling method. As shown in Figure 11, the pixel NTL time series after gap-filling captures the trend features retained in the source data. It displays excellent temporal continuity, thereby confirming the superiority of the proposed gap-filling method in this paper. In addition, the results of the time series analysis of four representative landmark buildings confirm the Bi-LSTM model’s capability in learning and predicting short-term trend features. These results emphasize the robustness and reliability of our proposed method based on Bi-LSTM in handling missing data and accurately forecasting time series patterns.

Figure 11. The gap-filling results for temporal analysis in iconic city locations: (a) Beijing Capital International airport, (b) Beijing Business District, and (c) Olympic Park. The images on the left represent the high-resolution satellite image in the target area. The plots on the right represent the time series of NTL pixels after gap-filling in the target locations in 2017. The blue curve represents the actual NTL pixel time series values, and the red curve represents the gap-filled NTL pixel time series values for the missing observations.

4.2.2 Quantitative assessment

Figure 12 presents the results of the proposed gap-filling method with varying missing rates. The original gap-free NTL images can be seen in Figure 12a. The black frames in Figure 12a were artificially created gaps in the original NTL images with different missing rates. Figure 12b shows the gap-filled images, and it is apparent that all the gaps well preserved the details of the original NTL images. The relationship between the NTL observations and the gap-filled NTL predictions with varying missing rates is shown in Figure 12c. The scatterplots of the original NTL values and gap-filled NTL pixel values demonstrated good pixel-wise performance. The R² values were between 0.770 and 0.850 with different missing rates, indicating that different missing rates caused minor effects on gap-filling performance. The gap-filled NTL image had the highest accuracy at a missing rate of 50% (R² = 0.841,Bias = −0.732, RMSE = 7.171) and the lowest accuracy at a missing rate of 10% (R² = 0.771,Bias = −2.921, RMSE = 12.512). Despite the high missing rate of 50%, we observed that the scatterplot of the NTL observations and the gap-filled NTL predictions was concentrated around the 1:1 line (Figure 12c), which further indicated the good performance of our proposed method.

Figure 12. Accuracy comparison of the gap-filled NTL data using the proposed gap-filling method ranging from 10% to 50% missing rates. The black frame represents artificially created gaps in the original NTL data, (a) represents the observed NTL image, (b) represents the gap-filled NTL image, and (c) represents the corresponding scatterplot and accuracy evaluation. N represents the total number of pixels participating in the scatterplot evaluation.

In this section, our main goal is to validate and demonstrate the rationality and superiority of our proposed algorithm. To achieve this, we conducted a comprehensive comparison with the Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM) (Gao et al. 2006), a widely recognized spatiotemporal fusion model. This model takes into account both spatial and temporal variability in the fusion process and is one of the most commonly used spatiotemporal fusion approaches. Unlike the application of the STARFM algorithm, which is mainly used to fuse MODIS and Landsat images with different resolutions, in this scenario, the spatiotemporal fusion module of the STARFM algorithm is employed to predict missing NTL values. The STARFM algorithm is the original version of the rough NTL spatiotemporal gap-filling method. It uses spatiotemporal information and fills gaps by calculating weight indicators. The weight indicators in the algorithm structures of these two algorithms are different. Due to the setting of the original parameters of the algorithm, the STARFM algorithm selects a pair of NTL images with fewer missing values from all eight-day NTL images for the gap-filling process, whereas the rough NTL spatiotemporal gap-filling method utilizes all eight-day NTL images to fill missing values.

Figure 13 presents the comparison results for four evaluated cities, showing that both the STARFM method and our proposed method exhibit improved spatial gap-filling details compared to the original images. However, our proposed gap-filling method demonstrates superior performance, particularly in areas with high luminance within the city’s built-up areas. Furthermore, Figure 13 displays the scatter plot evaluation, which reveals that our proposed gap-filling method achieves significantly higher R² values compared to those of the STARFM method. Specifically, the R² values for Shanghai and New York are 0.781 and 0.642, respectively, providing compelling evidence for the superiority of our gap-filling method.

Figure 13. Comparison of the accuracy of the gap-filling method proposed in this paper with the STARFM method. Four cities, Beijing, Shanghai, Xi’an and New York, were selected for comparison and verification. The blank in the original image shows the missing areas, and 205,47,103,38 represents the day of year (DOY).

4.3 Model robustness of sequence missing

To evaluate the accuracy of the Bi-LSTM model in predicting continuous missing NTL values and to assess the performance of the refined gap-filling model based on Bi-LSTM under varying levels of missingness, we tested it in study area of Beijing. As one of the largest cities in China, the spatial structure of Beijing exhibits a complex urban configuration characterized by multiple levels, multiple centers, and multiple functions. By conducting experiments in this region, we are able to examine the robustness of the Bi-LSTM model in extreme data missing situations. In specific, we selected sample points with four different levels of missing data ranging from 3 to 12 days with a 3-day interval from the annual average NTL images in study area of Beijing. To ensure that the sample points selected are representative of the entire study area and are evenly distributed across different levels of urban development, we have used the natural segmentation method to stratify the annual average NTL image value into five NTL brightness categories, namely class1 (0–5.5), class2 (5.5–18.3), class3 (18.3–34.8), class4 (34.8–71.4), and class5 (71.4–233.6). The classification of class1 to class5 represents the categorization of different luminance regions. Then we selected sample points from each category for subsequent evaluation.

According to Table 1, the accuracy of the Bi-LSTM model was affected by the increase in successive missing days and NTL brightness in study area of Beijing. To quantitatively compare the performance of the Bi-LSTM model, we computed the MAE, Bias, and RMSE for sample points from four different levels of missing days. As shown in Table 1, the model’s MAE, Bias, and RMSE increase accordingly with the length of the missing days to be filled and the NTL brightness values of the different classes. Specifically, the Bi-LSTM model predicted the best in three continuous missing days belonging to class 1 with an RMSE of 0.035 and MAE of 0.005. However, this method predicted the worst in 12 continuous missing days belonging to class 5 with an RMSE of 38.597 and MAE of 24.041.

Table 1. Evaluation of gap-filling accuracy based on the Bi-LSTM model at different missing lengths in study area of Beijing.

Missing length	Index	Class1	Class2	Class3	Class4	Class5
3 days	Bias	−0.004	0.042	0.028	0.307	3.332
	MAE	0.005	0.144	0.262	0.863	11.980
	RMSE	0.035	1.108	1.756	4.105	21.672
6 days	Bias	0.063	0.865	2.392	5.199	11.675
	MAE	0.465	2.696	5.478	8.922	16.738
	RMSE	0.958	4.297	7.489	13.163	24.401
9 days	Bias	−0.282	0.766	1.438	4.298	11.799
	MAE	0.748	4.381	6.824	10.727	24.654
	RMSE	1.493	6.396	9.655	16.417	43.677
12 days	Bias	−0.253	0.744	1.460	4.514	11.754
	MAE	0.744	3.952	6.772	10.952	24.041
	RMSE	1.528	5.745	9.15	18.611	38.597

From Table 1, we can see that as the length of the missing days increases, the gap-filling error becomes larger, and this is because as the successive missing days in NTL data increase, errors from the previous time step gap-filling result are transferred to the next time step, lowering the accuracy of the Bi-LSTM model. Additionally, high-brightness areas often indicate bustling commercial areas with complex human activities, leading to high spatial heterogeneity in NTL images. Spatial heterogeneity and angle effects of tall buildings in commercial areas increase the uncertainty of the NTL inversion and make daily NTL data gap-filling more challenging.

To demonstrate the effectiveness of the Bi-LSTM model in dealing with continuously missing time series data in study area of Beijing, we conducted experiments to validate its capacity to capture time series patterns and changes when 3, 6, 9, and 12 consecutive days of data are missing, and some sample points with varying lengths of missing data were selected from the table and visualized. Figure 14 demonstrates that the gap-filled values based on the Bi-LSTM model deviate within a reasonable range from the actual data and accurately capture the trends and patterns of the NTL data even in the case of continuous data gaps. This further supports the robustness of the proposed gap-filling model and its potential application in areas with continuous missing data.

Figure 14. Comparison of the gap-filled NTL value and actual NTL value, and rows (a) to (d) represent four specific sample points, indicating the gap-filling results by the model at different missing lengths in study area of Beijing.

5. Discussion

5.1 Why add a Bi-LSTM network to the gap-filling model and distinguish between refined and rough gap-fillings?

After NTL data has been gap-filled with the rough spatiotemporal gap-filling method, there are still two problems. Firstly, when a large fraction of the NTL image is missing for more than seven consecutive days, the rough spatiotemporal gap-filling method cannot be performed due to the unavailability of pixels in the spatial and temporal window. Secondly, in months where NTL data are significantly missing, the gap-filling performance for NTL time series with the rough spatiotemporal gap-filling method can be negatively affected. To address these issues, the Bi-LSTM model is proposed in the refined gap-filling method, which leverages the advantage of the Bi-LSTM model in handling time series data to learn the long-term trend and short-term characteristics of NTL time series. The crucial aspect of effectively utilizing deep learning network to learn the short-term characteristics of NTL time series data lies in the selection strategy for training samples and the incorporation of temporal prior information. NTL time series data shows a strong correlation with human activities, particularly in urban areas where human activities exhibit distinct periodic patterns. By including the daily, weekly and monthly temporal information of NTL time series data in the training samples, the model can adeptly learn the short-term change patterns in NTL time series data. Additionally, the training sample selection strategy involves using a sliding window approach to filter and select samples from NTL time series. This method ensures that the time series samples genuinely utilize the temporal information from the 365-day images in the test area. Regarding spatial considerations, the annual mean image value of the test area is initially divided into five categories using the natural breakpoint method. Subsequently, a hierarchical sampling technique is applied to uniformly select sample points from different layers, resulting in a final set of 11,837 samples for model training. After applying a 25-day sliding window to the 365-day time series images in the test area, 4,019,523 time series samples are obtained for participation in model training. These two approaches ensure that the sample selection considers areas with varying intensities of human activity across time and space. Consequently, the model can achieve superior gap-filling results across regions with different brightness levels.

The Bi-LSTM model in the refined gap-filling method is used to fill the remaining large gaps iteratively while updating the previous rough spatiotemporal gap-filling method results. We also compared the accuracy of the spatiotemporally weighted gap-filling method with the Bi-LSTM method. Figure 15 shows an average accuracy improvement of 0.1 in R² for the Bi-LSTM gap-filling model when the images are missing between 10% and 50%, in comparison to the spatiotemporally weighted gap-filling method. This result demonstrates the improved gap-filling quality after applying the Bi-LSTM network. The refined gap-filling method based on the Bi-LSTM model is crucial for addressing persistent data gaps and improving the accuracy of the rough spatiotemporal gap-filling method. In addition, incorporating the Bi-LSTM model strengthens the resilience of the gap-filling model and enhances its ability to handle extreme data absence.

Figure 15. Accuracy comparison of the gap-filled NTL data using between the rough NTL spatiotemporal gap-filling method and refined gap-filling method based on the Bi-LSTM model ranging from 10% to 50% missing rates.

5.2 Potential for monitoring sudden urban short-term dynamic changes

Worldwide urbanization has led to diverse changes in urban land use and land cover (LULC), and the temporal features of NTL time series can be used to distinguish various urban land changes and to record patterns of human activity and socioeconomic features (Li et al. 2022; Zheng, Weng, and Wang 2021). Coupling the Bi-LSTM model allows us to take advantage of Bi-LSTM in temporal prediction, making it possible to predict spatially discontinuous scenes with temporal gaps in remote sensing data. By including daily, weekly and monthly temporal information of images in the Bi-LSTM model’s training sample, the model can effectively capture short-term changes in the NTL time series. As a result, the gap-filled image time series maintains reliable trend characteristics and ensures that it accurately reflects the temporal pattern of the NTL data. Comparing our proposed gap-filling data with existing gap-filling products, it stands out in preserving essential time series characteristics. Figure 10 clearly demonstrates its capability to capture short-term fluctuations effectively. This is particularly valuable for applications with high timeliness requirements, such as disaster assessment, as well as for evaluating urban development status, urban function, and urban spatial structure. The ability to preserve temporal information makes our method highly suitable for studies requiring real-time analysis and monitoring weekly and monthly change patterns of urban land changes and human activities. After filling the gaps in NTL data, the daily NTL time series can better capture the change characteristics of the NTL time series within a year and even detect the monthly and weekly change patterns of urban land changes and human activities. To verify the ability of daily NTL sequences to monitor short-term dynamic changes, we used daily seamless NTL data after gap-filling to evaluate the impact of the short-term migration of large populations during the 2017 Chinese New Year holiday.

We selected four NTL images showing the research area before and after the Chinese New Year. These images compare Figure 16(a) before the Chinese New Year holiday and Figure 16(d) after the Chinese New Year holiday. Notably, the NTL brightness value of commercial and residential areas in the Beijing Chaoyang district exhibited much higher after the holiday than before. The reason can be attributed to the migration of many non-local workers who returned to their hometowns three days before the Chinese New Year holiday, leading to a decline in commercial and economic activities. Once the holiday ended and people resumed work, there was an increase in passenger flow intensity at transportation hubs. As observed in Figure 11(d), the NTL brightness at location B of the capital airport suddenly experienced an extreme increase.

Figure 16. The impact of the short-term migration of a huge population during the Chinese New Year holiday on NTL images. Site A and B correspond to the areas of Beijing Capital International airport and Beijing Business District, respectively.

Secondly, Figures 16b,c precisely depict two moments: the scene of people celebrating throughout the night on Chinese New Year’s Eve and the sudden decrease in urban population due to the holiday break and population migration, resulting in a decrease in the level of activity intensity in the city’s business districts. This short-term change can be reflected in the trend pattern of the daily NTL time series of the target pixel and the spatial change of NTL brightness. Such monitoring capability can be used for subsequent analysis to assess urban development status, city functions, and urban spatial structure, providing effective spatiotemporal features.

It is worth noting that the quality of NTL data is susceptible to the influence of other factors, such as cloud cover, lunar irradiance, and stray light. However, after filling gaps in NTL data with our proposed method, researchers can synthesize different temporal lengths of composite data based on the final daily gapless NTL data generated by the proposed method to meet the needs of various research applications.

In the past, filling gaps for remote sensing images using traditional empirical statistical methods often heavily relied heavily on a large amount of high-quality training data and auxiliary data (Li et al. 2013; Shen et al. 2011). However, since NTL data is relatively new and lacks auxiliary information, mining the spatiotemporal characteristics of the missing data itself and utilizing emerging deep-learning methods can significantly reduce the model’s dependence on multi-source data. The inclusion of Bi-LSTM model has two main advantages. First, using the temporal prediction ability of Bi-LSTM, we can temporally fill remote sensing images that are completely missing in space, which ensures the robustness and efficiency of the neural network in accurately filling the target area. Second, by incorporating daily, weekly and monthly temporal information in the training samples, the Bi-LSTM model ensures that the short-term changes in the NTL time series are efficiently captured. Therefore, the gap-filled image time series retains reliable trend features. In summary, the Bi-LSTM model proved to be a good fit for our specific gap-filling prediction scenario, providing accurate and reliable results. The parameters of the network model can be adjusted according to the test area and specific gap-filling requirements. For larger-scale regions, further experiments can help determine the optimal number of layers in a deep learning model, taking into account available training samples and region sizes. This relatively universal gap-filling method can be applied to NTL data and provides new insights and methods for the inversion and reconstruction of other remote sensing parameters.

5.3 Limitation of the proposed method

Although the proposed gap-filling method can achieve great accuracy in filling gaps for NTL data when there is a large fraction of missing observations, it also has some limitations. First, our proposed gap-filling method performs poorly in capturing precise variation patterns in regions with high spatial heterogeneity. High brightness observations in NTL data indicate frequent economic and commercial activities and typically correspond to the core functional area of the city, which often includes various urban functional elements. Therefore, capturing the precise variation pattern within regions of high spatial heterogeneity is challenging when only using the temporal characteristics of NTL data and a limited number of samples with high NTL brightness.

Secondly, the Bi-LSTM model in the refined gap-filling method only considers the feature patterns in the temporal domain and not the spatial domain, which may lead to the gap-filled results lacking fine spatial details. Therefore, future research can investigate the potential benefits of employing recurrent neural network models with spatial convolution information to capture more relevant and precise spatiotemporal features for gap-filling models. The lack of support from multi-source remote sensing data related to NTL is a significant issue. The availability and reliability of existing remote sensing data products, such as Land Surface Temperature (LST) and PM2.5 data, pose significant challenges due to limitations at specific sites, sensor capabilities on satellites, cloud cover issues, and data gaps. Large-scale application of these auxiliary data sources becomes less feasible due to these constraints. However, with rapid advancements in remote sensing technology and data gap-filling techniques, there is hope for future release of relevant gap-filling products by scholars. These products could provide valuable multi-source data support, helping us address current limitations and enhance the accuracy and efficiency of the NTL gap-filling process. In addition, incorporating other datasets related to NTL data can offer more relevant features and a comprehensive understanding of urban development, city functions, and urban spatial structure. The proposed method can be further improved and applied to a wider range of applications by addressing these limitations.

6. Conclusions

Our study presents an innovative solution to the problem of many missing observations in NASA’s Black Marble VNP46A2 product by introducing a gap-filling technique to obtain continuous NTL data. This approach combines spatiotemporal weighting and deep learning model to exploit the spatiotemporal features of the data, which can achieve good performance and model generalization without redundant auxiliary data. To evaluate the performance of our proposed gap-filling method, four cities with different regions in Beijing, Shanghai, Xi‘an, and New York, were tested.

The accuracy of the proposed method is evaluated using the “remove- reconstruct -compare” approach, using metrics such as RMSE, R² and MAE. Results indicate that the proposed method can fill up to 50% of image gaps with accuracy comparable to lower gap percentages and capture city boundaries and highly lit urban areas even with most of the data missing.

The proposed gap-filling method performed better than STARFM, and it offers several advantages, including better gap-filling accuracy and spatial details. These results indicate that this gap-filling method can provide a novel technique for filling gaps in the NASA VNP46A2 product.

Highlights

• A novel gap-filling method to fill gaps of daily Black Marble nighttime light (NTL) data was proposed.

• The spatiotemporal weighted features and Bi-LSTM model were utilized.

• Seamless daily Black Marble NTL data was produced.

• The method demonstrated good performance in four testing sites.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study are available upon request by contact with the corresponding author, or accessed through https://earthexplorer.usgs.gov.

Additional information

Funding

The work was supported by the Key Laboratory of Territorial Spatial Planning and Development-Protection of the Ministry of Natural Resources of PRC and CAUPD Beijing Planning & Design Consultants LTD (grant number TSPDP23/03), and the Fundamental Research Funds for the Central Universities (grant number 2652018077).