Correcting the Pixel Blooming Effect (PiBE) of DMSP-OLS nighttime light imagery

Remote Sensing of Environment 240 (2020) 111707

Contents lists available at ScienceDirect

Remote Sensing of Environment journal homepage: www.elsevier.com/locate/rse

Correcting the Pixel Blooming Effect (PiBE) of DMSP-OLS nighttime light imagery

T

⁎

Qiming Zhenga,b, Qihao Wenga, , Ke Wangb a

Center for Urban and Environmental Change, Department of Earth and Environmental Systems, Indiana State University, Terre Haute, IN 47809, USA Institute of Applied Remote Sensing and Information Technology, College of Environmental and Resource Sciences, Zhejiang University, Hangzhou, Zhejiang 310058, China

b

A R T I C LE I N FO

A B S T R A C T

Edited by Emilio Chuvieco

In the last two decades, the advance in nighttime light (NTL) remote sensing has fueled a surge in extensive research towards mapping human footprints. Nevertheless, the full potential of NTL data is largely constrained by the blooming effect. In this study, we propose a new concept, the Pixel Blooming Effect (PiBE), to delineate the mutual influence of lights from a pixel and its neighbors, and an integrated framework to eliminate the PiBE in radiance calibrated DMSP-OLS datasets (DMSPgrc). First, lights from isolated gas flaring sources and a Gaussian model were used to model how the PiBE functions on each pixel through point spread function (PSF). Second, a two-stage deblurring approach (TSDA) was developed to deconvolve DMSPgrc images with Tikhonov regularization to correct the PiBE and reconstruct PiBE-free images. Third, the proposed framework was assessed by synthetic data and VIIRS imagery and by testing the resulting image with two applications. We found that high impervious surface fraction pixels (ISF > 0.6) were impacted by the highest absolute magnitude of PiBE, whereas NTL pattern of low ISF pixels (ISF < 0.2) was more sensitive to the PiBE. By using TSDA the PiBE in DMSPgrc images was effectively corrected which enhanced data variation and suppressed pseudo lights from non-built-up pixels in urban areas. The reconstructed image had the highest similarity to reference data from synthetic image (SSIM = 0.759) and VIIRS image (r = 0.79). TSDA showed an acceptable performance for linear objects (width > 1.5 km) and circular objects (radius > 0.5 km), and for NTL data with different noise levels (< 0.6σ). In summary, the proposed framework offers a new opportunity to improve the quality of DMSP-OLS images and subsequently will be conducive to NTL-based applications, such as mapping urban extent, estimating socioeconomic variables, and exploring eco-impact of artificial lights.

Keywords: Nighttime light imagery DMSP-OLS Pixel Blooming Effect Point spread function Urban areas

1. Introduction Nighttime lights (NTL), as a unique indicator of human footprints, have been extensively utilized in diverse fields of researches, such as mapping urban growth (Cao et al., 2009; Imhoff et al., 1997), disaggregating environmental and socioeconomic metrics (Bennett and Smith, 2017; Zhao et al., 2018) and assessing the ecological consequences of artificial light encroachment in natural ecosystems (Gaston et al., 2015). Defense Meteorological Satellite Program-Operational Linescan System (DMSP-OLS) has collected NTL images from 1992 to 2013 and released them as an annually composite product, providing a valuable long-term observation portraying artificial lights from earth surface (Elvidge et al., 1997; Small et al., 2005). Yet, the full potential of DMSPOLS image data has not been realized because of the inherent

limitations of DMSP-OLS sensors, especially the blooming effect (Henderson et al., 2003). The blooming effect is a unique sensor effect of NTL data which causes actual lit areas to be enlarged in NTL images. The blooming effect has given rise to a succession of adverse impacts on NTL-based applications. For example, one of the most striking problems in using DMSP-OLS data is that the urban area derived from DMSP-OLS image tends to be overestimated (Xie and Weng, 2017). Various methods have been developed to correct the blooming effect. The most widely-used thresholding method attempts to correct the blooming effect by isolating pixels with a lower detection frequency or NTL intensity than a specific threshold value (Henderson et al., 2003; Zhou et al., 2014). Another popular method is image classification that often incorporates additional features, e.g., Vegetation Adjusted NTL Urban Index (VANUI) (Zhang et al., 2013) into NTL data to identify and remove the blooming effect (Cao et al., 2009).

⁎ Corresponding author at: Center for Urban and Environmental Change, Department of Earth and Environmental Systems, Indiana State University, Terre Haute, IN 47809, USA. E-mail addresses: [email protected] (Q. Zheng), [email protected] (Q. Weng).

https://doi.org/10.1016/j.rse.2020.111707 Received 29 May 2019; Received in revised form 18 December 2019; Accepted 4 February 2020 0034-4257/ © 2020 Elsevier Inc. All rights reserved.

Remote Sensing of Environment 240 (2020) 111707

Q. Zheng, et al.

Fig. 1. Transects of different nighttime light data (a). Nighttime light images of Berlin Tegel Airport, Germany, derived from DMSPstl (b), VIIRS (c), Aerial photo (Kuechly et al., 2012) (d) and photo from International Space Station (ISS) (e).

Effect (PiBE)”. PiBE emphasizes the fact that the blooming effect exists for all DMSP-OLS pixels rather than only in city outskirts. It is necessary to propose this term to distinguish from the traditional concept of the blooming effect, because the proposed methodology for eliminating PiBE is completely different from traditional methods, e.g., Henderson et al. (2003). Due to the PiBE, the detected light of a center pixel appears to “spill over” to its neighboring pixels, while lights from its neighboring pixels also scatter back to the center pixel. Consequently, the PiBE smooths data variation of DMSP-OLS and gives a pseudo indication of lights in non-BU pixels both inside and outside urban areas. Abrahams et al. (2018) discussed the causes of the blooming effect, which is briefly summarized as follows:

However, most of the methods are more effective for correcting the blooming effect outside the urban boundary than inside urban areas (Zhou et al., 2018). This limitation is mainly because the objective of these methods is to identify and remove dimly-lit pixels as blooming effect. The bulk of dimly-lit pixels are located at the edge of a city, so previous methods could only deal with the blooming effect around city outskirts. Nevertheless, there are evidences indicating that pixels around city outskirts, as well as those inside urban areas, are subject to the blooming effect. First, bright but pseudo lights exist in non-built-up (non-BU) pixels inside urban areas, such as inner-city lakes, green spaces, and rivers (Small et al., 2005). In fact, these pixels were dark as revealed by NTL images obtained from a more advanced nighttime light sensor—Visible Infrared Imaging Radiometer Suite (VIIRS) of Suomi NPP satellite (Fig. 1a). In VIIRS image pixels with the blooming effect, both inside and outside the city, are dimly-lit so they can be easily identified and eliminated. In contrast, many non-BU pixels and their surrounding BU pixels are equally bright in DMSP-OLS image, making it difficult to separate pseudo lights from true lights in BU areas. Second, only little data variation could be observed in DMSP-OLS data partly due to the saturation of DMSP-OLS stable light product (DMSPstl). However, even with gain-adjusted global radiance calibrated DMSPOLS (DMSPgrc) data which is an unsaturated dataset (Hsu et al., 2015), it still fails to reveal the spatial variation of NTL (Fig. 1a). For example, at Tegel Airport in Berlin, Germany, the NTL variation is clearly shown in VIIRS image (Fig. 1c), aerial photo (Fig. 1d), and an image captured from International Space Station (ISS) (Fig. 1e). Nevertheless, the image from DMSP-OLS exhibited a distinct pattern with limited spatial variation (Fig. 1b).

(1) Field-of-view variation. When the scanner swung from nadir point to the further end, the circular ground scanning area (GSA) was stretched into an elliptical shape. Each elliptical GSA would thereby capture the light from multiple neighboring light sources (Abrahams et al., 2018). (2) Geolocation error. The DMSP-OLS sensors were developed based on sensor technology in the 1970s. Because only satellite navigation data were used, raw DMSP-OLS images were geometrically mismatched among all daily scenes (Imhoff et al., 1997). When these mismatched daily scenes were stacked into composite data, the actual size of the light footprint was enlarged (Elvidge et al., 2004). (3) On-board data management issues. The DMSP-OLS systems were only able to store a limited amount of information with integer values, so the data was top-coded to 6-bit. Meanwhile, the five by five blocks of fine spatial resolution mode data (0.56 km) were averaged to become “smoothed” data with coarse spatial resolution (2.8 km). The data was later processed and released by NOAA as composite data product at 1-km grid, offering a uniform grid cell size (Elvidge et al., 1997). This procedure brought about a shift in the nighttime light footprints. (4) Atmospheric scattering resulted in a slight displacement of ground

2. The Pixel Blooming Effect (PiBE) The blooming effect has only been investigated by a limited number of studies (Abrahams et al., 2018; Cao et al., 2019). Herein, we propose a term to explain this phenomenon, and we name it “Pixel Blooming 2

Remote Sensing of Environment 240 (2020) 111707

Q. Zheng, et al.

B = PSF ∗ X + E → b = Ax + e

scanning area (Abrahams et al., 2018; Henderson et al., 2003). (5) The annual composite DMSP-OLS data, accessible to the public, were created by stacking all cloud-free daily observations, and the problems identified in (1)–(4) were accumulated through the stacking operation and exacerbated the PiBE problem.

(1)

where B is a blurred image, i.e., DMSPgrc image; PSF is a point spread function, which is also termed a blurring kernel; X is the assumed “original” PiBE-free image; E is a noise term; ∗ denotes the convolution operation; A presents the corresponding blurring matrix derived from PSF. It denotes the convolution operation with PSF is equivalent to matrix multiplication (Ax) (see Supplementary file S1); and b, x, and e are the column-wise vector form of B, X and E, respectively. Hence, to correct the PiBE one needs to inverse the process expressed by Eq. (1), which is simply to deconvolve the blurred image. The overarching goal of this study is to develop an integrated framework to model and correct the PiBE in DMSP-OLS image which will improve the data quality and benefit NTL-based applications. To achieve this goal our study (1) models how PiBE functions on DMSPgrc image by fitting a Gaussian model, and (2) develops a two-stage deblurring approach (TSDA) to reconstruct the PiBE-free DMSPgrc image by recovering the data variation and suppressing pseudo lights in nonBU pixels.

Thus far, only a few studies have provided direct solutions for modeling or correcting the blooming effect at pixel level. Small et al. (2005) found a linear relationship between the diameter of an island and the width of the blooming effect. Abrahams et al. (2018) suggested the blooming effect as an image blurring problem and deblurred stable light DMSP-OLS images with Wiener deconvolution. Cao et al. (2019) proposed a self-adjusting model (SEAM) with a spatial response function to correct the blooming effect at pixel scale. Nevertheless, several limitations remain unsolved: (1) Previous studies do not model the blooming effect at pixel scale; (2) The impact of data saturation has not been taken into consideration when addressing the blooming effect; (3) The pseudo lights of non-BU pixels in urban areas are not fully suppressed; (4) The corrected images are generated by making some simple assumptions. For instance, the SEAM model assumed that only neighboring pixels brighter than the corresponding center pixel would cast the blooming effect (Cao et al., 2019); and (5) Some previous algorithms possess a heavy computation load because they iteratively operated over the whole NTL image. Therefore, an improved method for correcting the PiBE becomes imperative. In this study, we intend to correct the PiBE in DMSP-OLS data by modeling the PiBE as an image blurring problem (Fig. 2a). That is, the “original” image is blurred by the PiBE and the blurring process is regarded as a convolution operation with a point spread function (PSF) expressed by Eq. (1):

3. Materials and methods 3.1. Study area and datasets The Yangtze River Delta (YRD), located in the east coast of China, was selected as our experimental site. With over 150 million residents, the YRD accommodates 25 prefectural cities with different population levels and landscape layout patterns. Cities in the YRD also share a similar climate so that local differences can be minimized (Fig. 3). To deal with the data saturation problem of DMSP-OLS, a number of NTL-based indices have been proposed, such as VANUI (Zhang et al.,

Fig. 2. The conceptual framework (a), analytical procedures (b), and two models of the proposed framework (c) for modeling and correcting the PiBE. 3

Remote Sensing of Environment 240 (2020) 111707

Q. Zheng, et al.

Fig. 3. DMSP-OLS stable light product (2013) for the Yangtze River Delta, China. The red circles denote eight sample prefectural-level cities. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

3.2. Methods

2013) and HSI (Lu et al., 2008). However, the assumption behind these indices is that if a pixel has low NDVI (or, high impervious surface fraction) it should have high NTL value. Studies using high resolution NTL data reveal that this assumption might not be true in reality and NTL intensity varies largely across urban land use classes (Levin et al., 2014). Besides, the variation of NDVI in urban areas is limited compared with the dominantly high NTL values so they could not sufficiently increase NTL variation. In addition, different NTL-based indices provide different, sometimes even contradictory, NTL patterns (Ma et al., 2018). Thus, global radiance calibrated DMSP-OLS (DMSPgrc) data were employed in this study, instead of commonly used DMSP-OLS stable-light product (DMSPstl) or NTL-based indices. DMSPgrc data have a 30-arc spatial resolution and were produced by blending gain-adjusted scenes with DMSPstl. The DMSPgrc data are superior to the DMSPstl because the former one have a larger dynamic range (i.e., without data saturation) and are capable of retrieving dimer lights and revealing more NTL variations (Hsu et al., 2015) (refer to the blue line in Fig. 1a). The Version 1 monthly averaged VIIRS NTL images product was used to assess the results. Compared with DMSP-OLS, an annually averaged composite product, VIIRS has a more advanced sensor with a higher spatial resolution (15-arc sec) and temporal resolution (monthly), a larger dynamic range (14-bit) and a better low-light detection sensitivity (~2E-11 W/cm2/sr) (Elvidge et al., 2013). DMSPgrc, VIIRS, and global gas flaring record were downloaded from the Earth Observation Group (EOG) of NOAA (https://www.ngdc.noaa.gov/eog). Globeland30 land cover map (30 m, 2010) (Chen et al., 2014), MODIS land cover map product (MCD12Q1) were obtained to evaluate the accuracy of urban mapping. MCD12Q1 was also used as an auxiliary input data in the proposed method.

Fig. 2b shows the operational flowchart of our method. In short, we utilized DMSPgrc and the estimated PSF (Section 3.2.2) and a two-stage deblurring approach (TSDA) (Section 3.2.3) to remove the PiBE and reconstruct the “original” PiBE-free NTL image. Lastly, the resulting images were assessed in Section 4. 3.2.1. Preprocessing DMSPgrc and VIIRS images were both projected into the Albers Equal Area Conic projection system. They were resampled to 500 m resolution with the cubic convolution method to ensure sufficient data points available for PSF estimation and to retain the high data quality of VIIRS. The geometric correction and inter-calibration were applied to DMSPgrc data (Zhao et al., 2015). The DMSPgrc data are generated from unique fixed-gain DMSP-OLS scenes, so DMSPgrc are only available in eight discrete years: 1996, 1999, 2000, 2003, 2004, 2006, 2010, and 2011 (Hsu et al., 2015). A logistic regression model was employed to interpolate DMSPgrc into a consecutive annual data from 1996 to 2013 to ensure a sufficient number of valid gas flare sites can be found and to support long-term applications of the resulting datasets (Zheng et al., 2019). The monthly VIIRS images were aggregated into a yearly averaged dataset. The DMSPgrc and VIIRS in 2013 were used to test and evaluate our method, respectively. 3.2.2. Modeling the point spread function (PSF) The PSF manifests how the blurring process, i.e., the PiBE, functions on each pixel. It is difficult to estimate the PSF directly from urban areas because the NTL footprints are already intertwined with each other. Light from isolated gas flare has a great potential to depict the shape of underlying PSF. Gas flare sites are usually located in remote and offshore areas, making them less likely to be affected by neighboring lighting sources (Elvidge et al., 2004). Two distinct merits of gas 4

Remote Sensing of Environment 240 (2020) 111707

Q. Zheng, et al.

Fig. 4. The selected single-point lighting sources and an example of gas flare identified in Google Earth historical image (a); scatter plot of the estimated ax and ay (b), and the valid points (black) and outlier points (gray) beyond mean ± 2σ (the blue dash-line square); the corresponding DMSPgrc image of the zoom-in image (c); and the estimated PSF (d). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

flares were found from our benchmark tests: (1) the shape of estimated PSF became stable when at least 25 gas flare sites were used; and (2) the shape of estimated PSF was not correlated to the brightness, spatial location, extent and time of gas flare in statistical significance (see Supplementary file S2). It is worth noting that what PSF modeled (i.e., PiBE) is a combinational result of multiple sensor effects instead of the process of a certain cause. From the perspective of PiBE's causes, PSF should be irrelevant to the type of detected light emission. Thus, it is reasonable to assume that gas flares are an ideal source to estimate PSF and the resulting PSF can be used for modeling PiBE of other light sources in addition to gas flares. To select isolated gas flares, we first employed the gas flare map and DMSPgrc image in 2004 to identify candidates. Then, Google Earth (GE) historical images in the corresponding year were used to double-check whether a candidate gas flare was isolated and whether the light was emitted by mining facilities (Fig. 4). The estimated shape of PSF would be inevitably affected by geometric distortion of the projection system, but it could be minimized by selecting the gas flare sites within two standard parallels of Albers projection systems (ESRI, 2018). The DMSPgrc images of the identified gas flare sites were fitted by a Gaussian model (Eq. (2)). The rationale of using Gaussian model was that the PiBE might not be spread equally along the longitude direction and latitude direction, while the Gaussian model was able to capture this discrepancy with ax and ay. In our case, 29 valid single-pixel lighting sources (black points) were selected, which showed a good fitting performance with an R2 of 0.981 ± 0.012 (Fig. 4b). Given that the Gaussian model had a double-symmetric structure, NTL intensity was not relevant to the spatial location of the center of Gaussian surface. Thus, x0 and y0 were set as the coordinate origin, and the shape of

the Gaussian surface (i.e., PSF) was only determined by ax and ay. Lighting sources, whose ax and ay were beyond mean ± 2σ, were considered as outliers and discarded. A universal PSF was determined by the average ax (3.26 pixels) and ay (3.65 pixels) of the 29 valid lighting sources and was then normalized to ensure the sum of PSF to equal to 1.

PSF = a0 ×

⎧ (x − x 0 )2 (y − y0 )2 ⎫ − + 2 ⎨ 2a y2 ⎬ ⎭ e ⎩ 2ax

(2)

where x and y are the coordinates of the pixel; a0 is a magnitude factor; x0 and y0 are the center location of Gaussian surface; ax and ay are the half-long and half-short axes. The size of the PSF was set to 15 × 15 pixels (7.5 × 7.5 km) after trial and error tests. In a nutshell, we examined the corresponding deblurring performance with a series of PSF window sizes (5–25 pixels). We found that it was necessary to avoid small PSF size, which resulted in an undersmoothed reconstructed image with a fragmented NTL pattern. And when the window size was sufficiently large (> 11 pixels), the impact of PSF size on the deblurring performance was negligible. 3.2.3. Correcting the PiBE by TSDA Now that the PSF (i.e., matrix A) was known, the next step was to reconstruct the original PiBE-free image. However, it was impossible to reconstruct the exact original image (Hansen et al., 2006), because the blurred image contained a noise component (A-1e) and a noise-free blurred image (A-1bexact) (Eq. (4)). If we directly calculated the inverse of Eq. (1) through Eq. (3), the reconstructed image would be overwhelmed by noises and lose all useful information. 5

Remote Sensing of Environment 240 (2020) 111707

Q. Zheng, et al.

(3)

x = A−1 b x=

A−1 bexact

+

A−1 e

(oversmoothed). On the other hand, if α is small, the magnified noises would be undersmoothed.

(4)

GCV(α) =

Stage 1: Deconvolution with Tikhonov regularization (simple model) Singular value decomposition of matrix A is used to further explain the noise-overwhelming problem (Eq. (5)). As the singular value σi became smaller, |uiTe| remained approximately in the same magnitude order. The high-frequency noise component was substantially magnified when divided by small singular value (σi), so the reconstructed image was dominated by noises. Thus, eliminating the expansion of inverted noise component played a determinant role in eliminating the PiBE. N

x = VΣ −1U T b =

∑ i=1

uiT b vi = σi

N

∑ i=1

uiT bexact vi + σi

N

∑ i=1

uiT e vi σi

N

∑ φi i=1

σi2/(σi2

uiT b vi σi

(7)

Stage 2: Suppressing pseudo lights in non-BU pixels Although the NTL variation could be recovered by optimal regularization parameter α1, non-BU pixels in urban areas were still as bright as the surrounding BU pixels. It was difficult to discriminate effectively these pseudo lights from the true lights of BU pixels. If these pseudo lights were not properly corrected, it would lead to errors when NTL data were used, such as commission error of urban mapping (i.e., mislabeling non-BU pixels as BU pixels) or spurious indication of high population density in non-BU areas. To suppress pseudo lights of nonBU pixels, a second-stage deconvolution was implemented. We aimed to find the second optimal Tikhonov regularization parameter α2, where α2 was a bit smaller than α1. Deconvolving with α2 would lead to an undersmoothed reconstructed image, partially contaminated by noises, but the NTL intensity difference between non-BU and BU pixels could be sufficiently enlarged in the reconstructed image. Then, a single predefined threshold could identify and mask out pseudo lights of nonBU pixels easily. Detailed steps of Stage 2 deblurring process are as follows:

(5)

where U and V refer to the left and right singular matrix, respectively; Σ denotes the diagonal matrix of singular value σi in a non-increasing order (i = 1,2,…,N), so σ1 ≥ σ2 ≥ … ≥σN; and N is the number of singular values; and u and v denote vector form of U and V, respectively. The Tikhonov regularization filter was employed to solve this problem (Tikhonov et al., 2013) (Eq. (6)).

x reg = VΦΣ −1U T b =

‖b − Axreg ‖22 ‖(IN − AVΦΣ−1U T)b‖22 = 2 N (trace(IN − AVΦΣ−1U T))2 (N − ∑i = 1 φi )

(6) Step 1: Creating a non-BU buffer zone at urban boundary

+ α ); Φ is the diagonal matrix of filter factors φi; α where φi = denotes the regularization parameter and α ∈ [σN, σ1]; Xreg stands for the reconstructed image. The resultant φi ≈ 1 for large σi, denoting large singular values with useful information would be assigned a high weight close to 1. When φi decreased sharply (φi ≈ σi2/α2) as σi became smaller, small singular values would be assigned a small weight close to 0. This small weight suppressed the noise expansion and could reserve part of useful information for small singular values instead of dumping them all. An ideal α was expected to recover most of the lost data variation and constrain the expansion of noise component simultaneously. The Generalized Cross Validation (GCV) method was adopted to find the optimal value of α (α1) that gave the minimum GCV value for the GCV- α curve (Eq. (7)) (Golub et al., 1979). As Fig. 5 shows if α is too large the reconstructed image would be still blurred 2

A non-BU buffer zone at urban boundary was created. MCD12Q1 was utilized to establish a rough non-BU buffer outside the BU areas. MCD12Q1 land cover map is one of the most widely-used land cover map products, whose user's accuracy and producer's accuracy of urban and built-up land (IGBP-scheme) are 99.2% and 92.7%, respectively (Sulla-Menashe et al., 2019). Therefore, the impact of classification uncertainty on TSDA would be negligible. When MCD12Q1 data were not available, we suggest using other land cover map products (e.g., CCI Land cover map from ESA: https://www.esa-landcover-cci.org) or manually interpret the non-BU buffer zone. Since our testing cities are all prefectural level cities in China that shared a similar development stage, the universal size of 5-pixel was used for the non-BU buffer.

Fig. 5. An illustration of using the GCV-α curve to determine the optimal α (α1). 6

Remote Sensing of Environment 240 (2020) 111707

Q. Zheng, et al.

PSF rather than with a huge-size matrix A. The computation process can thereby be greatly accelerated (Hansen et al., 2006). 4. Experiments 4.1. Experiment with synthetic data Two simulation tests were designed to evaluate the performance of the proposed deblurring method with a known PSF. First, the PiBE is associated with the spatial configuration between BU and non-BU pixels, while urban morphology is a key determining factor for this configuration relationship. Therefore, we simulated the reference PiBEfree image and the PiBE-contaminated image with two urban morphology prototype scenarios, i.e., concentric and polycentric model (Lynch, 1984). Two typical cities, Chicago (US) and Wuhan (China) were selected to represent each scenario (Liu and Wang, 2016). Ground truth BU and non-BU areas were extracted from the MCD12Q1 land cover map product. The synthetic reference (PiBE-free) image was generated by randomly assigning 1–3 urban centers with the highest NTL intensity, where NTL intensity decreased from urban center to urban outskirt. The NTL intensity of non-BU pixels was set to zero. The synthetic blurred image was created by adding normally distributed noises and then using estimated PSF to convolve the synthetic reference image. Second, because the noise magnitude of DMSPgrc was unknown, the effectiveness of TSDA was further examined under different noise levels ~N(μ = 0, σ = k ∙ σref), where k is the noise level factor and σref denotes the standard deviation of the reference image. Moreover, to further explore the effectiveness of our method in dealing with small non-BU objects inside urban areas, a sensitivity analysis was conducted by evaluating the performance upon different shapes of non-BU objects (e.g., circle or linear), sizes of non-BU objects, and distance to urban center (Fig. 6). Linear and circular objects are commonly selected to represent small objects in digital image analysis (Zhu et al., 2010). We created linear objects and circular objects with different sizes (i.e., radius or width) and at different distance to urban center, and then evaluated the effectiveness of TSDA for these objects. Each simulation process was repeated 500 times to reduce the impact of simulated noise with extreme values and then assessed by the mean value of the following two indicators. Structural SIMilarity index (SSIM) was designed to measure the quality of the deblurred (reconstructed) image from the perspective of luminance, contrast and structure (Wang et al., 2004). A larger SSIM implies a higher similarity to the reference image of the reconstructed image. Meanwhile, a Coefficient of Relative Variation (CRV) was established to quantify the relative data variation between the reconstructed image and the reference image (Eq. (8)).

Fig. 6. A diagram for sensitivity analysis of TSDA when being applied to circular and linear non-BU objects at different sizes and distances to urban center (red point) which was set with the highest NTL intensity in the simulated image. The dashed square denotes BU area and the area encompassed by solid lines are non-BU areas. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Step 2: Deconvolving with another regularization parameter α2 The DMSPgrc image was then deconvolved by a certain α2 where α2 ∈ [σN, α1] to obtain the reconstructed image. Step 3: Finding an optimal α2 to effectively identify pseudo lights The average NTL intensity of the buffer zone in the reconstructed image was set as a threshold. Pixels with NTL intensity lower than the threshold were labeled as non-BU areas while other BU areas. The optimized α2 was determined by repeating step 2 and step 3 and lowering the value of α2 and the α2 that maximized the correlation coefficient between the BU/non-BU binary map derived from step 3 and MCD12Q1. Step 4: Generating the PiBE-free image Finally, non-BU pixels identified with the optimized α2 in step 3 were masked out in the reconstructed image obtained from Stage 1 deblurring. After this two-stage deblurring process, the resulting image was regarded as a PiBE-free DMSPgrc image. To distinguish between them, we named the Stage 1 deblurring process as a simple deblurring model, as the resulting image was only partially deblurred; and the combination of Stage 1 and Stage 2 deblurring process (i.e., TSDA) as a full deblurring model. Considering that the size of blurring matrix A is huge when applied to a large study area, it would require a large computational load to implement singular value decomposition. In this study, the PSF estimated by the Gaussian model had a double-symmetric structure and matrix A was constructed by a reflexive boundary condition. Consequently, the blurring matrix A was a normal matrix, which can be uniquely represented by a small number of entries through discrete cosine transformation (DCT), so that efficient approaches can take the place of computationally expensive Singular Value Decomposition (SVD). Using DCT bases allowed us to operate with a 15-by-15 pixel's

CRV =

σntl σref / μntl μref

(8)

where μ and σ denote the mean and standard deviation; and ntl and ref refer to the testing image and reference image, respectively. 4.2. Experiments and application tests Evaluating reconstructed PiBE-free NTL imagery was not a straightforward task because no ground-truth PiBE-free DMSP-OLS imagery was available to serve as reference data. Although there are many high spatial resolution NTL data, VIIRS was the most suitable one to serve as an alternative “reference” image due to its data availability, spectral range and low light detection sensitivity. It should be further noted that the image quality of VIIRS is better than DMSP-OLS, suggesting VIIRS can only be used for comparison rather than as groundtruth data. Considering the ultimate goal of correcting the blooming effect is to support NTL-based applications, two experiments were carried out to test the applicability of the deblurred data: (1) Urban 7

Remote Sensing of Environment 240 (2020) 111707

Q. Zheng, et al.

Fig. 7. An example of the simulation test for concentric city (a–d) and polycentric city (e–h) (noise level = 0.4∙σref). S: SSIM; C: CRV. The white squares in (a) are four sample non-BU areas.

areas inside the city. The image reconstructed by TSDA produces a higher correlation coefficient to reference image than simple model because TSDA successfully removed a majority of pseudo lights of nonBU pixels that failed to be suppressed by the simple deburring model. The simulation test indicates that our method could enhance NTL variation and recover a considerable portion of lost information from the PiBE-contaminated image. The transects and frequency histograms further explain why deconvolving with a small regularization parameter (α2) could suppress pseudo lights in non-BU pixels (Fig. 9). Although deconvolving with an optimal α (α1) was helpful to balance noise expansion and lost information recovery, the NTL intensity difference between BU and nonBU pixels was inadequately enlarged (the gray line in Fig. 9a). A large number of non-BU pixels, especially those in urban areas, were still as bright as BU pixels in the blurred image and image reconstructed by simple model (purple bars in Fig. 9b & c). This explained why pseudo lights of those non-BU pixels were inclined to be confused with BU pixels in previous methods. When being deconvolved with a proper α2, NTL difference between BU and non-BU pixels was evidently increased (green line in Fig. 9a). Most BU pixels were brighter than non-BU pixels in the image reconstructed in Stage 2 deblurring process (Fig. 9d). Our method failed in dealing with a small proportion (5.42%) of non-BU areas in urban centers, which were extremely bright or had a narrow size. The brightness of these non-BU and their surrounding BU pixels were nearly identical, thereby making the deblurring algorithm less effective. Sensitivity analysis illustrates the effectiveness of TSDA against different types of non-BU objects at different noise levels (Fig. 10). The average correlation coefficients of all testing linear objects were higher than 0.5, and reached as high as 0.87 ± 0.04 for the ones wider than 1.5 km (Fig. 10a). For the circular object with a radius larger than one pixel (0.5 km), the correlation coefficients yielded 0.88 ± 0.09 (Fig. 10b). TSDA shows a stable and acceptable deblurring performance (SSIM > 0.75) when the noise level is lower than 0.6σref. SSIM decreased rapidly as noise level became larger than 0.6σref (concentric scenario) and 0.7σref (polycentric scenario) (Fig. 10c) where all useful information of the reconstructed image was lost, except for the noises. This was not the case for DMSPgrc image as a general spatial variation pattern of NTL intensity could at least be observed in DMSPgrc. It further implies that the noise level of DMSPgrc, although undetermined, was

mapping. We used the reconstructed image to map urban extent by using the local optimal threshold method (Liu et al., 2012). The ground truth urban extent was derived from GlobeLand30 and MCD12Q1 land cover map product and visual interpretation on Landsat images. Kappa coefficient and the exaggeration rate were computed for accuracy assessment. (2) Modeling the Gross Domestic Product (GDP) dynamics. We adopted linear regression to model the relationship between the sum of total light and the GDP of each prefectural city from 1996 to 2013 (Shi et al., 2014). The GDP data of each city were obtained from Provincial/City Statistical Yearbook of China. The model performance was evaluated by the determination coefficient (R2). 5. Results 5.1. Result with the synthetic data Fig. 7 demonstrates the result of the simulation test. Subject to the PiBE, the spatial variation of DMSPgrc data was noticeably smoothed in the blurred image (Fig. 7b & f). Taking synthetic polycentric city scenario as an example (Fig. 7e–h), BU pixels appeared to be equally bright and the data variation in the blurred image was nearly half of the synthetic reference image (CRV = 0.563). In contrast, the reconstructed images (Fig. 7g & h) showed an improved data variation, and a NTL gradient from urban center to the periphery could be clearly observed. Compared with the simple deblurring model (CRV = 0.640), the full model was able to recover more data variation (CRV = 0.889). A similar case was also found in concentric city scenario (Fig. 7a–d). More importantly, our method performed well in suppressing pseudo lights in non-BU pixels. For the synthetic blurred image and image reconstructed by the simple model, non-BU pixels in urban areas were as bright as surrounding BU pixels. For instance, the linear object (see the white square in Fig. 7g) across urban centers is even brighter than most of the BU pixels in urban outskirts. Attempting to remove these pixels in Fig. 7g would lead to an overcorrection of the PiBE by mislabeling 38.28% of the BU pixels as non-BU pixels. Most of the pseudo lights in non-BU pixels were successfully identified and removed through the full deblurring model (TSDA). The sharp increase of SSIM (0.735 and 0.759) implies a higher resemblance between reference images and images reconstructed by TSDA (Fig. 7d & h). Fig. 8 provides a close look at the performance of the deblurring model over four sample non-BU 8

Remote Sensing of Environment 240 (2020) 111707

Q. Zheng, et al.

Fig. 8. Correlation coefficients of four sample areas in Fig. 7a between the synthetic reference image and the testing image, including the blurred image (B), the image reconstructed by the simple model and TSDA (full model).

more likely to be lower than 0.7 σref. Thus, TSDA was expected to work well for DMSPgrc data. The performance of TSDA was slightly degraded for objects far from urban centers because these synthetic objects had a low NTL intensity and were apt to be overcorrected during Stage 2 deblurring process. Besides, the TSDA was less effective for narrow linear objects (width ≤ 3 pixels, 1.5 km) and small circular objects (radius ≤ 1 pixel, 500 m). The NTL intensity of these tiny objects was extremely high, so they were difficult to be separated from surrounding BU pixels, especially for the ones surrounded by bright pixels, like in urban centers (e.g., sample area 4 in Fig. 7a).

5.2. Result with the VIIRS data Fig. 11a demonstrates the relationship between impervious surface fraction (ISF) and PiBE intensity where DMSPPiBE = |DMSPgrc, TSDA − DMSPgrc|, and DMSPgrc,TSDA and DMSPgrc denote the normalized NTL intensity of the image reconstructed by TSDA and DMSPgrc, respectively. The relative PiBE intensity (DMSPPiBE, rel) was determined by the ratio of DMSPPiBE and DMSPgrc. The ISF of each 500 m pixel was computed from GlobeLand30 land cover data (30 m). It was found that pixels with high ISF (ISF > 0.6) were subject to intensive light scattering from neighboring pixels, whereas the PiBE only contributed to a

Fig. 9. NTL transects from different datasets of sample area 1 in Fig. 7a (a). Frequency histogram for non-BU and BU pixels: blurred image (b), image deblurred by simple deblurring model (c) and image deblurred by Stage 2 deblurring (d). 9

Remote Sensing of Environment 240 (2020) 111707

Q. Zheng, et al.

Fig. 10. The correlation coefficient of different linear objects (a) and circular objects (b). SSIM values for different noise levels (k·σref) (c).

Fig. 11. The relative PiBE intensity (primary y-axis) and the absolute PiBE intensity (secondary y-axis) against different categories of the impervious surface fraction (ISF) (a). The correlation coefficient (r) and coefficient of relative variation (CRV) between VIIRS (as reference data) and different testing datasets (b). VANUI (Zhang et al., 2013), SEAM (Cao et al., 2019) and Wiener filter (Abrahams et al., 2018) are employed with DMSPstl and DMSPgrc as input data for comparison.

variation. By removing the PiBE the reconstructed image recovered a considerable amount of the variation inside urban areas. As a result, local NTL peaks, not discernible in either DMSPstl or DMSPgrc images, could now be identified, such as commercial areas in Fig. 12g & h and transportation hubs in Fig. 12e. The TSDA (full model) achieved a better deblurring quality (r = 0.79, CRV = 1.49) than the simple deblurring model (r = 0.72, CRV = 1.20) by suppressing pseudo lights of non-BU pixels in urban areas. For instance in Fig. 12a, because of the existence of PiBE, lights from bright arterial traffic lanes and commercial buildings smeared into the neighboring wetland. Consequently, the wetland turned out to be equally bright as the neighboring BU pixels in DMSPstl and DMSPgrc images, whereas the wetland was supposed to be dark as shown in the corresponding VIIRS image. After applying TSDA, the illumination discrepancy between the dark wetland and the surrounding bright traffic lanes was clearly presented, which were largely in line with the NTL pattern in the VIIRS image. Similar cases could also be found, such as the small hill in Fig. 12f and river in Fig. 12d. In addition, we downloaded the source code of SEAM (Cao et al., 2019) and the Wiener filter (Abrahams et al., 2018) from https:// xiaolinzhu.weebly.com and https://github.com/alexeiabrahams, respectively, and performed the same experiment with DMSPstl and DMSPgrc as input data (Fig. 11b). It was found that compared with the blurred DMSPstl and DMSPgrc image the Wiener filter and SEAM were able to reduce PiBE, but their performance was inferior to TSDA. Besides, using DMSPgrc (r = 0.79, CRV = 1.49), instead of DMSPstl (r = 0.60, CRV = 0.91), as input data brought a more noticeable improvement in TSDA performance than other compared methods. Our further investigation revealed the limitations of these two methods. First, both methods failed to suppress pseudo lights of non-BU pixels in

minor NTL intensity increase for medium ISF pixels (0.2 < ISF < 0.6). Although high ISF pixels (ISF > 0.6) had the highest PiBE magnitude (DMSPPiBE), the PiBE only accounted for 31% of the total NTL intensity changes. Because the baseline NTL intensity of high ISF pixels was high, it offset the NTL changes resulted from the PiBE. Interestingly, the PiBE rendered a notable influence on low ISF pixels (ISF < 0.2), which appeared to have a moderate DMSPPiBE but the highest relative PiBE intensity (DMSPPiBE,rel) among other ISF groups. This finding does not only underline that the blooming effect occurs in all pixels of DMSP-OLS images but suggests that NTL intensity pattern of low ISF pixels was sensitive to the impact of the PiBE. Therefore, special precaution should be paid to the PiBE in non-BU pixels (ISF < 0.2 is a common indication of non-BU pixels), otherwise the PiBE would bring about dramatic changes in their NTL intensity pattern. Fig. 11b illustrates the CRVs and correlation coefficients between the reference data VIIRS and different testing datasets, and Fig. 12 gives a detailed comparison of different datasets of eight sample cities in YRD. First, the image quality of DMSPstl-based datasets was found poorer than that of DMSPgrc-based datasets. Due to PiBE and data saturation problem, DMSPstl showed the lowest data variation (r = 0.34, CRV = 0.27) and the entire urban and peri-urban area appeared to be dominated by high NTL values (Fig. 12). It was observed that substantial pseudo lights caused by the PiBE “scattered” into the surrounding non-BU areas, such as Wuxi, Changzhou and Suzhou. Even with the DMSPstl image enhanced by NDVI (e.g., VANUI), the data variation remained low (r = 0.42, CRV = 0.32). This implies that NTLbased indices may not be suitable to be used as input data, as they could worsen the deblurring quality. Second, the image reconstructed by both the simple model and the full model resembled with VIIRS image, and had relatively high data

10

Remote Sensing of Environment 240 (2020) 111707

Q. Zheng, et al.

Fig. 12. Comparison between different NTL images and reconstructed PiBE-free images. JL1-3B is a high-resolution NTL image with 0.9 m spatial resolution (b). For each city, all the nighttime light images are in the same scale, spatial resolution and stretch method. 11

Remote Sensing of Environment 240 (2020) 111707

Q. Zheng, et al.

suppressed pseudo lights. Therefore, the resulting urban extent was more accurate and more comparable to the urban extent derived from VIIRS which possessed the lowest exaggeration rate (7.96% ± 1.83%) and the highest Kappa coefficient (0.72 ± 0.05). Another test for modeling GDP dynamics with NTL data also corroborated that the deblurred image was beneficial to NTL based applications (Table 2). The deblurred image shows the highest R2 (0.822 ± 0.079), compared to DMSPstl (0.504 ± 0.133), DMSPgrc (0.676 ± 0.091), and the image reconstructed by the Wiener filter (0.721 ± 0.035) and SEAM (0.767 ± 0.072).

urban areas. Second, the deblurring quality was found sensitive to the PSF setting and related parameters, such as the signal-noise-ratio used in the Wiener filter. The above findings highlight the effectiveness of TSDA for eliminating the PiBE and the necessity for using DMSPgrc to avoid data saturation in the deconvolution process. Nevertheless, as the simulation test indicated pseudo lights from narrow objects were not adequately suppressed, particularly for features close to bright urban centers. Fig. 12b gives an example and a probable explanation in a real situation. The reason for this ineffectiveness was that there were ten bridges connecting urban business centers at both sides of the Qiantang River across Hangzhou city. Previous studies using high spatial resolution NTL images had reported the illumination facilities for main traffic arteries and car headlights contributed to a considerable proportion of total upward light emission (Kuechly et al., 2012), and the percentage was over 58% for Hangzhou (Zheng et al., 2018). Bright lights scattered from bridges and surrounding business centers together made it challenging to remove the pseudo lights of the PiBE over a narrow river. However, being away from the influence of luminosity from urban center, another part of the river lit-up by the PiBE in the same city, appeared to be successfully discriminated from the bright bridges (Fig. 12c). An over-correction of PiBE was observed in Changzhou (Fig. 12i). The overcorrection problem occurs in the built-up area which is usually close to urban outskirts and has a lower NTL intensity than its surrounding. In the case of Changzhou city two holes were a combination of non-built-up areas (bare land and green space) and less developed built-up areas alongside a canal (i.e., the Datong Canal) crossing through the city. Thus, when we implemented the Stage 2 deblurring process to suppress pseudo lights in the non-built-up area, it overcorrected the PiBE by masking out the pseudo light and dim light in adjacent built-up areas (see Supplementary file S3). Moreover, although the reconstructed DMSPgrc image recovered part of lost spatial variation, it did not match exactly the VIIRS image. VIIRS image was able to provide more details than what the reconstructed DMSPgrc image could reveal. The differences between the reconstructed DMSPgrc image and VIIRS image can be attributed to two reasons: (1) VIIRS was only employed for comparison rather than as reference data because VIIRS and DMSP-OLS data are essentially different. VIIRS is more advanced than DMSP-OLS regarding its spatial and temporal resolution, dynamic range, and low light detection limit. The overpass time of VIIRS (~1:30 am) and DMSP-OLS (~7:30 pm) is also different (Elvidge et al., 2013). Although the general quality of the reconstructed PiBE-free DMSP-OLS image has been improved, it is still lower than VIIRS image. Even under the most ideal condition the reconstructed DMSPgrc and VIIRS are not likely to be identical; and (2) DMSPgrc data, generated by blending DMSPstl images and the fix-gain scenes, are a product of DMSP-OLS. It is inevitable for DMSPgrc to inherit part of the drawbacks from DMSP-OLS sensors, thereby limiting the amount of information that could be recovered through TSDA.

6. Discussion and conclusions This study probed into the PiBE of DMSP-OLS data and proposed an integrated framework to eliminate the impact of pseudo lights from neighboring pixels. The point spread function (PSF) was estimated by 29 isolated single-pixel lighting sources. We then developed a two-stage deblurring approach (TSDA) to reconstructed PiBE-free DMSPgrc. NTL pattern of low ISF pixels was found noticeably vulnerable to the impact of PiBE. The evaluation results from synthetic datasets and VIIRS image both showed that the image reconstructed by the TSDA had the highest similarity to the synthetic reference image (SSIM = 0.759/0.735) and the highest correlation coefficients with VIIRS image (r = 0.79). TSDA was effective in recovering the lost data variation (CRV = 1.49) and suppressing lights in inner-city non-BU pixels. The sensitivity analysis revealed that the TSDA had an acceptable performance for linear objects wider than 1.5 km and circular objects with a radius larger than 0.5 km, and when the noise level was lower than 0.6 σref. Our findings indicate that the proposed method was a feasible solution to model and correct the PiBE of DMSPgrc imagery. The advantages of our method can be summarized as follows. First, removing PiBE is of critical necessity for NTL based studies. On one hand, the enhanced NTL variation makes it possible to uncover more spatial variation of NTL data. On the other hand, the bright pseudo lights of inner-city non-BU pixels are masked out through Stage 2 deblurring process. It significantly lowers the possibility of mislabeling non-BU pixels as BU pixels, when applied to map urban extent. It also has great potential to benefit other NTL-based applications, such as estimating socioeconomic variables. For example, if the PiBE of non-BU pixels (e.g., a lake closed to the urban center) is not properly suppressed, the estimated population would turn out to be as high as that of the densely populated residential areas. Methods, like Wiener deconvolution employed in Abrahams et al. (2018) can only recover the lost information in the blurred image. Those methods are unable to suppress PiBE in non-BU pixels inside urban areas, which, as shown in Fig. 11a, are the most vulnerable to PiBE. Second, we provide a series of top-down approaches to estimate the PSF (i.e., Gaussian model) and to obtain the optimal parameters for Tikhonov regularization (i.e., GCV curve). This avoids a subjective choice of modeling parameters and is good for the quality of the deblurring process. Third, DMSPgrc is used as the input data instead of DMSPstl. The data saturation would strongly impair the quality of the reconstructed image by limiting the total variation that could be recovered and bringing about a ringing effect in the reconstructed image (Ji and Wang, 2012). Furthermore, our method is computationally efficient. Accelerating by the Discrete Cosine Transformation, it took < 10 min to implement the entire framework for all the cities in YRD (2123×2391 pixels) using a quad-core laptop (2.8 GHz, Intel® Core™ i7-4810MQ, 16GM RAM). At the same time, several limitations should also be noted. First of all, in this study, it is assumed that all pixels were subject to the same kind of PiBE, which can be modeled by a spatial invariant model. To adopt the spatial variant PSF model is a new trend in image deblurring studies. However, the image deblurring problem is not exactly the same as the PiBE of NTL data from the point of the origin of blurring. The image deblurring approaches are mostly designed for blurred images caused by motion blurring from ordinary cameras (Ji and Wang, 2012).

5.3. Applications The accuracy assessment of urban mapping with different testing datasets is listed in Table 1. The urban extent derived from DMSPstl and DMSPgrc resulted in a much higher exaggeration rate, 38.70%, and 24.13%, respectively. About 36% of the exaggeration rate of DMSPgrc image was incurred by the commission error of non-BU pixels in urban areas. A city with a large percentage of non-BU areas (e.g., lakes and large-scale greenspaces) tended to suffer from a higher exaggeration rate, such as Hangzhou (43.16%), Nanjing (32.61%), and Suzhou (29.50%). Regardless of the thresholding-based or a classification-based method, it was difficult to discriminate BU pixels from non-BU with pseudo lights. The urban extent obtained from the image reconstructed by TSDA showed a much lower exaggeration rate (8.96% ± 2.08%) and a higher accuracy for urban mapping (Kappa = 0.68 ± 0.05). This improvement was mainly due to the enlarged data variation and 12

Remote Sensing of Environment 240 (2020) 111707

Q. Zheng, et al.

Table 1 Exaggeration rate and Kappa coefficient of urban extent derived by different datasets. City

Hangzhou Shanghai Nanjing Yancheng Ningbo Suzhou Wuxi Changzhou YRD⁎

City

Hangzhou Shanghai Nanjing Yancheng Ningbo Suzhou Wuxi Changzhou YRD ⁎

Exaggeration rate DMSPstl

VANUI

DMSPgrc

TSDA

VIIRS

43.16% 31.01% 32.61% 29.46% 28.48% 29.50% 31.51% 30.40% 38.70% ± 4.42%

32.98% 15.81% 28.02% 24.82% 24.51% 25.15% 25.07% 23.35% 27.23% ± 4.78%

25.14% 19.91% 16.76% 20.45% 26.00% 20.15% 26.67% 19.83% 24.13% ± 4.12%

8.35% 9.10% 9.23% 7.01% 7.59% 6.05% 5.53% 7.68% 8.96% ± 2.08%

6.99% 5.88% 5.05% 6.42% 7.53% 4.12% 5.82% 7.14% 7.96% ± 1.83%

Kappa coefficient DMSPstl

VANUI

DMSPgrc

TSDA

VIIRS

0.37 0.48 0.55 0.53 0.44 0.50 0.47 0.53 0.48 ± 0.06

0.58 0.63 0.60 0.59 0.63 0.65 0.66 0.65 0.62 ± 0.04

0.56 0.63 0.56 0.64 0.56 0.61 0.56 0.59 0.59 ± 0.04

0.60 0.66 0.68 0.69 0.60 0.59 0.65 0.70 0.68 ± 0.05

0.68 0.70 0.76 0.76 0.61 0.71 0.68 0.75 0.72 ± 0.05

The statistics for YRD are calculated for all prefectural-level cities (n = 25) in YRD.

Table 2 The determination coefficient (R2) of the regression between the sum of total night-time light (TNL) and GDP in 1996–2013. City

Hangzhou Shanghai Nanjing Yancheng Ningbo Suzhou Wuxi Changzhou YRD⁎ ⁎

R2 (TNL-GDP) DMSPstl

DMSPgrc

Wiener

SEAM

TSDA

0.570 0.363 0.562 0.660 0.508 0.369 0.591 0.522 0.504 ± 0.133

0.680 0.621 0.769 0.587 0.667 0.399 0.586 0.726 0.676 ± 0.091

0.707 0.723 0.771 0.731 0.639 0.713 0.724 0.739 0.721 ± 0.035

0.834 0.766 0.672 0.833 0.799 0.752 0.808 0.622 0.767 ± 0.072

0.769 0.875 0.721 0.690 0.835 0.744 0.784 0.859 0.822 ± 0.079

The statistics for YRD are calculated for all prefectural-level cities (n = 25) in YRD.

process is established to identify and to mask out the pixels with pseudo lights. We set the NTL intensity of these pixels to zero, so TSDA, in some cases, overcorrected the PiBE. This operation should be helpful for most NTL-based applications, such as urban mapping, GDP or population estimation, because removing lights from non-BU pixels is an essential step for these applications (Bennett and Smith, 2017). However, we suggest to employ the simple model only, rather than the full model, for other studies that include assessing light pollution of natural landscapes and disaggregating energy consumption. In our future work, we plan to establish a predictive model to describe the relationship between urban land use and nighttime light emission, and use this model to further solve the overcorrection problem. Third, although TSDA significantly reduces the overcorrection problem, a small percentage of BU pixels are still masked out as non-BU pixels, which in this case was 2.98% for synthetic data and was a little bit higher for DMSPgrc data. The mislabeled BU pixels were mainly dimly-lit highways and rural settlements around cities. In summary, the proposed framework offers a new approach to deal with the blooming effect of DMSP-OLS data at pixel scale and opens an opportunity to improve the quality of DMSP-OLS data, as well as subsequent applications. The reconstructed PiBE-free DMSPgrc image could be conducive to use in conjunction with VIIRS images to explore the

The PiBE of NTL imagery is induced by the combination of multiple sensor effects in each scene and the image stacking process. Meanwhile, NTL imagery is taken under a complicated imaging condition with low luminosity and a substantial amount of noise, such as lunar contamination and atmospheric scattering (Elvidge et al., 2017). The deblurring approach used in this study is not among the most advanced ones in light of how we model the PiBE and how we reconstruct the PiBE-free image. We compared the Tikhonov filter used in TSDA with other widely used deconvolution methods—Truncated Singular Value deconvolution and Wiener filter deconvolution. The proposed method showed a better performance in dealing with the PiBE. The superior of TSDA lies in that we employ a top-down method to obtain the optimal regularization parameter and Tikhonov regularization can retain part of useful information when suppressing noise expansion (see Supplementary file S4). A further in-depth investigation is needed to explore the effectiveness of more advanced deblurring methods in addressing the PiBE. Second, two stages of deblurring model are proposed to eliminate the PiBE. An appropriate method should be selected according to the objective of an application. Specifically, the simple model was found effective in recovering lost data variation but would under-correct pseudo lights of narrow non-BU areas in urban areas. Stage 2 deblurring 13

Remote Sensing of Environment 240 (2020) 111707

Q. Zheng, et al.

long-term dynamics of nighttime lights. In future works, we intend to (1) explore the possibility of establishing a spatial variant PSF model by examining the contributing factors of PiBE and developing an algorithm to eliminate the impact of noises and outliers; and (2) improve the deblurring method for small non-BU objects and to avoid overcorrecting the light from BU pixels.

Filtering. 3 Siam. Henderson, M., Yeh, E.T., Gong, P., Elvidge, C., Baugh, K., 2003. Validation of urban boundaries derived from global night-time satellite imagery. Int. J. Remote Sens. 24 (3), 595–609. Hsu, F.-C., Baugh, K.E., Ghosh, T., Zhizhin, M., Elvidge, C.D., 2015. DMSP-OLS radiance calibrated nighttime lights time series with intercalibration. Remote Sens. 7 (2), 1855–1876. https://doi.org/10.3390/rs70201855. Imhoff, M.L., Lawrence, W.T., Stutzer, D.C., Elvidge, C.D., 1997. A technique for using composite DMSP/OLS “city lights” satellite data to map urban area. Remote Sens. Environ. 61 (3), 361–370. Ji, H., Wang, K., 2012. A two-stage approach to blind spatially-varying motion deblurring. In: Paper Presented at the Computer Vision and Pattern Recognition (CVPR), 2012 IEEE Conference on. Kuechly, H.U., Kyba, C.C.M., Ruhtz, T., Lindemann, C., Wolter, C., Fischer, J., Hölker, F., 2012. Aerial survey and spatial analysis of sources of light pollution in Berlin, Germany. Remote Sens. Environ. 126, 39–50. https://doi.org/10.1016/j.rse.2012.08. 008. Levin, N., Johansen, K., Hacker, J.M., Phinn, S., 2014. A new source for high spatial resolution night time images — the EROS-B commercial satellite. Remote Sens. Environ. 149, 1–12. https://doi.org/10.1016/j.rse.2014.03.019. Liu, X., Wang, M., 2016. How polycentric is urban China and why? A case study of 318 cities. Landsc. Urban Plan. 151, 10–20. https://doi.org/10.1016/j.landurbplan.2016. 03.007. Liu, Z., He, C., Zhang, Q., Huang, Q., Yang, Y., 2012. Extracting the dynamics of urban expansion in China using DMSP-OLS nighttime light data from 1992 to 2008. Landsc. Urban Plan. 106 (1), 62–72. https://doi.org/10.1016/j.landurbplan.2012.02.013. Lu, D., Tian, H., Zhou, G., Ge, H., 2008. Regional mapping of human settlements in southeastern China with multisensor remotely sensed data. Remote Sens. Environ. 112 (9), 3668–3679. https://doi.org/10.1016/j.rse.2008.05.009. Lynch, K., 1984. Good City Form. MIT press. Ma, T., Xu, T., Huang, L., Zhou, A., 2018. A human settlement composite index (HSCI) derived from nighttime luminosity associated with imperviousness and vegetation indexes. Remote Sens. 10 (3). https://doi.org/10.3390/rs10030455. Shi, K., Yu, B., Huang, Y., Hu, Y., Yin, B., Chen, Z., ... Wu, J., 2014. Evaluating the Ability of NPP-VIIRS Nighttime Light Data to Estimate the Gross Domestic Product and the Electric Power Consumption of China at Multiple Scales: A Comparison with DMSPOLS Data. Remote Sensing 6 (2), 1705–1724. https://doi.org/10.3390/rs6021705. Small, C., Pozzi, F., Elvidge, C.D., 2005. Spatial analysis of global urban extent from DMSP-OLS night lights. Remote Sens. Environ. 96 (3–4), 277–291. https://doi.org/ 10.1016/j.rse.2005.02.002. Sulla-Menashe, D., Gray, J.M., Abercrombie, S.P., Friedl, M.A., 2019. Hierarchical mapping of annual global land cover 2001 to present: the MODIS collection 6 land cover product. Remote Sens. Environ. 222, 183–194. https://doi.org/10.1016/j.rse.2018. 12.013. Tikhonov, A.N., Goncharsky, A., Stepanov, V., Yagola, A.G., 2013. Numerical Methods for the Solution of Ill-posed Problems. 328 Springer Science & Business Media. Wang, Z., Bovik, A.C., Sheikh, H.R., Simoncelli, E.P., 2004. Image quality assessment: from error visibility to structural similarity. IEEE Trans. Image Process. 13 (4), 600–612. Xie, Y., Weng, Q., 2017. Spatiotemporally enhancing time-series DMSP/OLS nighttime light imagery for assessing large-scale urban dynamics. ISPRS J. Photogramm. Remote Sens. 128, 1–15. https://doi.org/10.1016/j.isprsjprs.2017.03.003. Zhang, Q., Schaaf, C., Seto, K.C., 2013. The vegetation adjusted NTL urban index: a new approach to reduce saturation and increase variation in nighttime luminosity. Remote Sens. Environ. 129, 32–41. https://doi.org/10.1016/j.rse.2012.10.022. Zhao, N., Zhou, Y., Samson, E.L., 2015. Correcting incompatible DN values and geometric errors in nighttime lights time-series images. IEEE Trans. Geosci. Remote Sens. 53 (4), 2039–2049. https://doi.org/10.1109/tgrs.2014.2352598. Zhao, N., Zhang, W., Liu, Y., Samson, E.L., Chen, Y., Cao, G., 2018. Improving nighttime light imagery with location-based social media data. In: IEEE Transactions on Geoscience and Remote Sensing, pp. 1–12. https://doi.org/10.1109/tgrs.2018. 2871788. Zheng, Q., Weng, Q., Huang, L., Wang, K., Deng, J., ... Jiang, R., Gan, M., 2018. A new source of multi-spectral high spatial resolution night-time light imagery—JL1-3B. Remote Sensing Of Environment 215, 300–312. https://doi.org/10.1016/j.rse.2018. 06.016. Zheng, Q., Weng, Q., Wang, K., 2019. Developing a new cross-sensor calibration model for DMSP-OLS and Suomi-NPP VIIRS night-light imageries. ISPRS J. Photogramm. Remote Sens. 153, 36–47. https://doi.org/10.1016/j.isprsjprs.2019.04.019. Zhou, Y., Smith, S.J., Elvidge, C.D., Zhao, K., Thomson, A., Imhoff, M., 2014. A clusterbased method to map urban area from DMSP/OLS nightlights. Remote Sens. Environ. 147, 173–185. https://doi.org/10.1016/j.rse.2014.03.004. Zhou, Y., Li, X., Asrar, G.R., Smith, S.J., Imhoff, M., 2018. A global record of annual urban dynamics (1992–2013) from nighttime lights. Remote Sens. Environ. 219, 206–220. https://doi.org/10.1016/j.rse.2018.10.015. Zhu, X., Chen, J., Gao, F., Chen, X., Masek, J.G., 2010. An enhanced spatial and temporal adaptive reflectance fusion model for complex heterogeneous regions. Remote Sens. Environ. 114 (11), 2610–2623. https://doi.org/10.1016/j.rse.2010.05.032.

Author contributions Qiming Zheng, and Qihao Weng designed research; Qiming Zheng performed data analysis and wrote the draft; Qihao Weng revised and edited the paper; Ke Wang provided comments on the paper. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgment We would like to thank anonymous reviewers and the editor for their constructive comments and suggestions. This study is supported in part by China Scholarship Council (Grant No. 201806320144), Zhejiang University (Grant No. 2018092), and USGS IndianaView Program via Purdue University. Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.rse.2020.111707. References Abrahams, A., Oram, C., Lozano-Gracia, N., 2018. Deblurring DMSP nighttime lights: a new method using Gaussian filters and frequencies of illumination. Remote Sens. Environ. 210, 242–258. https://doi.org/10.1016/j.rse.2018.03.018. Bennett, M.M., Smith, L.C., 2017. Advances in using multitemporal night-time lights satellite imagery to detect, estimate, and monitor socioeconomic dynamics. Remote Sens. Environ. 192, 176–197. https://doi.org/10.1016/j.rse.2017.01.005. Cao, X., Chen, J., Imura, H., Higashi, O., 2009. A SVM-based method to extract urban areas from DMSP-OLS and SPOT VGT data. Remote Sens. Environ. 113 (10), 2205–2209. https://doi.org/10.1016/j.rse.2009.06.001. Cao, X., Hu, Y., Zhu, X., Shi, F., Zhuo, L., Chen, J., 2019. A simple self-adjusting model for correcting the blooming effects in DMSP-OLS nighttime light images. Remote Sens. Environ. 224, 401–411. https://doi.org/10.1016/j.rse.2019.02.019. Chen, J., Ban, Y., Li, S., 2014. China: open access to Earth land-cover map. Nature 514 (7523), 434. Elvidge, C., Baugh, K., Kihn, E.A., Kroehl, H.W., Davis, E.R., 1997. Mapping city lights with nighttime data from the DMSP operational linescan system. Photogramm. Eng. Remote. Sens. 63 (6), 727–734. Elvidge, C.D., Safran, J., Nelson, I., Tuttle, B., Hobson, V.R., ... Baugh, K.E., Lyon, J., 2004. Area and positional accuracy of DMSP nighttime lights data. In Remote sensing and GIS accuracy assessment 281–292. Elvidge, C., Baugh, K., Zhizhin, M., Hsu, F.-C., 2013. Why VIIRS data are superior to DMSP for mapping nighttime lights. Proc. Asia Pac. Adv. Netw. 35 (0), 62. https:// doi.org/10.7125/apan.35.7. Elvidge, C., Baugh, K., Zhizhin, M., Hsu, F.C., Ghosh, T., 2017. VIIRS night-time lights. Int. J. Remote Sens. 38 (21), 5860–5879. https://doi.org/10.1080/01431161.2017. 1342050. ESRI, 2018. Albers equal area conic. Retrieved from. http://desktop.arcgis.com/en/ arcmap/latest/map/projections/albers-equal-area-conic.htm. Gaston, K.J., Gaston, S., Bennie, J., Hopkins, J., 2015. Benefits and costs of artificial nighttime lighting of the environment. Environ. Rev. 23 (1), 14–23. https://doi.org/ 10.1139/er-2014-0041. Golub, G.H., Heath, M., Wahba, G., 1979. Generalized cross-validation as a method for choosing a good ridge parameter. Technometrics 21 (2), 215–223. Hansen, P.C., Nagy, J.G., O’leary, D.P., 2006. Deblurring Images: Matrices, Spectra, and

14