Assessment of Night-Time Lighting for Global Terrestrial Protected and Wilderness Areas

Protected areas (PAs) play an important role in biodiversity conservation and ecosystem integrity. However, human development has threatened and affected the function and effectiveness of PAs. The Defense Meteorological Satellite Program/Operational Linescan System (DMSP/OLS) night-time stable light (NTL) data have proven to be an effective indicator of the intensity and change of human-induced urban development over a long time span and at a larger spatial scale. We used the NTL data from 1992 to 2013 to characterize the human-induced urban development and studied the spatial and temporal variation of the NTL of global terrestrial PAs. We selected seven types of PAs defined by the International Union for Conversation of Nature (IUCN), including strict nature reserve (Ia), wilderness area (Ib), national park (II), natural monument or feature (III), habitat/species management area (IV), protected landscape/seascape (V), and protected area with sustainable use of natural resources (VI). We evaluated the NTL digital number (DN) in PAs and their surrounding buffer zones, i.e., 0–1 km, 1–5 km, 5–10 km, 10–25 km, 25–50 km, and 50–100 km. The results revealed the level, growth rate, trend, and distribution pattern of NTL in PAs. Within PAs, areas of types V and Ib had the highest and lowest NTL levels, respectively. In the surrounding 1–100 km buffer zones, type V PAs also had the highest NTL level, but type VI PAs had the lowest NTL level. The NTL level in the areas surrounding PAs was higher than that within PAs. Types Ia and III PAs showed the highest and lowest NTL growth rate from 1992 to 2013, respectively, both inside and outside of PAs. The NTL distributions surrounding the Ib and VI PAs were different from other types. The areas close to Ib and VI boundaries, i.e., in the 0–25 km buffer zones, showed lower NTL levels, for which the highest NTL level was observed within the 25–100 km buffer zone. However, other types of PAs showed the opposite NTL patterns. The NTL level was lower in the distant buffer zones, and the lowest night light was within the 1–25 km buffer zones. Globally, 6.9% of PAs are being affected by NTL. Conditions of wilderness areas, e.g., high latitude regions, Tibetan Plateau, Amazon, and Caribbean, are the least affected by NTL. The PAs in Europe, Asia, and North America are more affected by NTL than South America, Africa, and Oceania.

surrounding areas [29]. Studies suggested that the NTL level of the boundary of PAs is particularly high [39]. In addition, some studies have combined NTL with other data to conduct research, for example, Geldmann et al. used the NTL and population data to construct a temporal human pressure index (THPI) on the time series, which quantifies the human pressure on the protected area. However, the spatial resolution of the study was 10 km 2 , and only attempted comparisons between 1995 and 2010 [40].
In this study, we used NTL data as an indicator of the intensity of human development to portray the characteristics of NTL in various types of PAs defined by the International Union for Conversation of Nature (IUCN). The objectives of this study were to: (1) reveal human development surrounding different types of terrestrial PAs, growth rate, and trend characteristics measured using the NTL level; and (2) make an assessment about the NTL effects on global terrestrial PAs.

Protected Areas and Data
We used the 2019 World Database on Protected Areas (WDPA, https://protectedplanet.net/), as the most updated data source for obtaining the coverage of global protected areas. We selected seven types of terrestrial PAs defined by the IUCN, including Ia, Ib, and II-VI, as defined in Table 1. In order to explore the spatial distribution and change of NTL within PAs and their vicinity, we examined buffer zones with distances of 0-1 km, 1-5 km, 5-10 km, 10-25 km, 25-50 km, and 50-100 km from PA boundaries. This range of buffer distances was chosen to capture different types of NTL effects. Human activities within PA boundaries and within a 1 km distance exert a direct and significant influence on protected areas, e.g., habitat loss, noise, and light pollution [41]. At further distances, artificial surfaces contribute to the landscape disturbances and effects, such as the isolation of protected areas, disruption of connectivity, and introduction and spread of invasive species. Even as far as 50 km from the protected area, human development can impose effects on PAs. PAs may be impacted even if the source of lighting lies kilometers away owing to skyglow [24]. PAs that are strictly set aside to protect biodiversity and also possibly geological/geomorphological features, where human visitation, use, and impacts are strictly controlled and limited to ensure the protection of the conservation values.
Wilderness Area Ib 3114 1,152,403 PAs that are usually large unmodified or slightly modified areas, retaining their natural character and influence, without permanent or significant human habitation.
National Park II 5523 6,178,074 Large natural or near natural areas set aside to protect large-scale ecological processes, along with the complement of species and ecosystems characteristic of the area, which also provide a foundation for environmentally and culturally compatible spiritual, scientific, educational, recreational, and visitor opportunities.
Natural Monument or Feature III 14,317 434,460 PAs set aside to protect specific natural monuments, which can be a landform; sea mount; submarine cavern; geological feature, such as a cave; or even a living feature, such as an ancient grove. They are generally quite small protected areas and often have a high visitor value.

DMSP/OLS NTL Data
The Operational Linescan System (OLS) is one of the sensors that is carried by the Defense Meteorological Satellite Program (DMSP). The DMSP/OLS NTL data came from the National Geophysical Data Center (NGDC) under the National Oceanic and Atmospheric Administration (NOAA), which eliminates accidental noise effects, such as clouds and flaring, with a 0.008333 degrees spatial resolution. The number of each pixel is not the actual light level on the ground, but rather the relative brightness level recorded as a digital number (DN) from 0 to 63. Currently, NTL data from 1992 to 2013 are available through online access. DMSP/OLS is different from Landsat, Satellite Pour l'Observation de la Terre (SPOT), and the Advanced Very High Resolution Radiometer (AVHRR) sensors that use the ground object to monitor the reflected radiation characteristics of sunlight. The DMSP/OLS sensor can work at night and capture the low intensity of urban lights and even small-scale residential areas and traffic lights, providing powerful data support for monitoring human development research on a large scale [42].

Data Processing
For origin data, discrepancies appeared between the DN values and the number of lit pixels from different satellites in the same year, and abnormal fluctuations appeared in the DN values for different years derived from the same satellite. It was necessary to calibrate the original NTL between different years and satellites. In this study, we assumed that all regions of the world that developed positive reflection in NTL would keep the DN value and the number of lit pixels either consequently increased or remained the same. The NTL were corrected using four main steps as described below.

Inter-Annual Correction
As the economy and population continue to develop over time, the assumption was that all PA regions would be positively developing and increasing in terms of the amount of NTL. Therefore, the DN values of each pixel in the time series would either increase or remain unchanged. With this assumption, the light pixels detected in the early image would remain in the latter image, and the DN value should not be larger than the DN value detected in their subsequent images. If the pixel that detected the light in the earlier image disappeared in the later year, the pixel with light would be Remote Sens. 2019, 11, 2699 5 of 20 considered as an unstable light pixel, and its DN value would be replaced with a value of 0. Second, the DN value of each stable light pixel was corrected to ensure that the DN value in the early image was smaller than the DN value of the corresponding position pixel in the later image. Due to the replacement of satellites between 1992 and 2013, the light sensitivity of each sensor was different, leading to the difference in the number of lit pixels and the DN values of the corresponding position pixels from different sensors. Thus, each satellite was interannually corrected using the following step [43]: where DN (n−1, i) , DN (n, i) , and DN (n+1, i) are the DN values of the ith lit pixel of the NTL data in the (n − 1)th, the nth, and (n + 1)th years, respectively.

Inter-Satellite Correction
The NTL from 1992 to 2013 were derived from multiple sensors. The data were not continuous and the data between different sensors could not be directly compared. It was necessary to correct the sensors to improve the continuity and comparability of the NTL. Although the image of satellite F18 did not intersect with those of the other satellites, there was an intersection between F10 and F12, F12 and F14, F14 and F5, and F15 and F16. First, based on F10, using the data of the intersection year, 100,000 pixels were randomly selected on each continent to establish a minimum linear-square regression model between two satellites; thus, the data of F12 could be corrected, which could then be used to correct the data of F14. The same process was applied to data from other satellite pairs.
The regression formulas are as follows: Correct F12 based on F10: Correct F14 based on F12: Correct F15 based on F14: Correct F16 based on F15:

Intra-Annual Correction
The intra-annual composition correction used images from two satellites in the same year to remove any intra-annually unstable lit pixels. There were two images from two different satellites in 1994 and 1997-2008. We used the average DN values of the two NTL images to calculate the annual NTL data for these years. First, we examined all the lit pixels to determine whether the pixel was an unstable light pixel during the year. If only one satellite was detected, the light pixel was defined as an unstable light pixel during the year. Second, in intra-annual composites, the DN values of intra-annually unstable lit pixels were replaced with values of 0, and the DN value of each intra-annually stable lit pixel was replaced by the average DN value of two NTL images from the same year. This produced one intra-annual composite for each year [43].
where DN a (n,i) and DN b (n,i) are the DN values of the ith lit pixel from two NTL data in the nth year, and DN (n,i) is the DN value of the ith lit pixel of the intra-annual composite in the nth year.

Correction for Different Sensors
We performed another first-step correction on this series of data after the first three steps, which resulted in a dataset that grew positively and with continuity.

Evaluation of the Calibrated NTL Time Series
A number of studies have shown that the DMSP/OLS NTL data are correlated with economic activities [44][45][46][47][48] and have used gross domestic product (GDP) and electricity consumption (EC) to evaluate the performance of intercalibration techniques [28,[49][50][51][52][53]. We applied this assessment method in this study to evaluate the calibrated NTL time series. We obtained country-level GDP data (1992-2013) from the World Bank (http://data.worldbank.org/) and EC records derived from International Energy Statistics (http://www.eia.gov/beta/international/data/browser/#). We compared the raw NTL and the calibrated NTL with GDP and EC data ( Figure 1) to reduce the inconsistency of the data and enhanced correlation with GDP and EC. We calculated the Pearson correlation between the raw NTL with GDP and EC, and between the calibrated NTL and the two ancillary datasets ( Figure 2). The calibrated method improved the correlations between the NTL time series and GDP and EC.  [44][45][46][47][48] and have used gross domestic product (GDP) and electricity consumption (EC) to evaluate the performance of intercalibration techniques [28,[49][50][51][52][53]. We applied this assessment method in this study to evaluate the calibrated NTL time series. We obtained country-level GDP data (1992-2013) from the World Bank (http://data.worldbank.org/) and EC records derived from International Energy Statistics (http://www.eia.gov/beta/international/data/browser/#). We compared the raw NTL and the calibrated NTL with GDP and EC data ( Figure 1) to reduce the inconsistency of the data and enhanced correlation with GDP and EC. We calculated the Pearson correlation between the raw NTL with GDP and EC, and between the calibrated NTL and the two ancillary datasets ( Figure 2). The calibrated method improved the correlations between the NTL time series and GDP and EC.

Calculation of the Trend Using Sen's Slope
The Sen's slope method can be used to avoid the loss of time series data and the influence of data distribution on the analysis results to eliminate the interference of outliers in a time series [54]. Sen's slope has been used in the analysis of long-term sequence data sets to detect the magnitude of the change trend. We used Sen's slope to calculate the trend of the NTL from 1992 to 2013. The calculation formula is: where Q is the Sen's slope; and Xj and Xi are the average DN value corresponding to j and i year, respectively. If the length of the time series is n, the number of Qi is N = n(n − 1)/2; and the Q is determined by N.
where Q > 0, Q < 0 and Q = 0 indicate that there is a rising trend, a downward trend, and no obvious trend, respectively.

Spatial Distribution and Trend of NTL for Global PAs
We superimposed the global PAs and the global NTL, except for Antarctica ( Figure 3). The result revealed that the DN values for most of the PAs were lower than that of the surrounding areas. Further assessment ranked the average NTL DN values of all PAs on each continent as: Europe > Asia > North America > South America > Africa > Oceania. The DN values in Europe was much higher than those in the other continents ( Figure 4). The ranking of the average change trend on all continents was consistent with the average DN value, i.e., the highest was in Europe and the lowest was in Oceania ( Figure 5).
We superimposed the global PAs and the global NTL, except for Antarctica ( Figure 3). The result revealed that the DN values for most of the PAs were lower than that of the surrounding areas. Further assessment ranked the average NTL DN values of all PAs on each continent as: Europe > Asia > North America > South America > Africa > Oceania. The DN values in Europe was much higher than those in the other continents ( Figure 4). The ranking of the average change trend on all continents was consistent with the average DN value, i.e., the highest was in Europe and the lowest was in Oceania ( Figure 5).    For (a), we used the natural breaks (Jens) method to divide the PAs into three levels according to the average trend of NTL. The blue-colored PAs had an average trend of 0, the yellowcolored PAs had a lower than average trend, and the red-colored PAs had a higher than average trend. For (b), the ranking of the average NTL trend of PAs in each continent gave the order as Europe > Asia > North America > South America > Africa > Oceania. Figure 6 illustrates that NTL DN values in types III and V PAs were significantly higher than that of the other five types of PAs in the interior and buffer zones. With the buffer radius increased, the average NTL DN values of all PAs increased first and then decreased (Figure 7). The highest DN values of different types of PAs appeared in the range of 1-10 km.

Average NTL Level in Different PAs and Buffers
From the time series, the NTL level within and outside PAs increased without exception, but the NTL distributions among seven types of PAs had different characteristics ( Figure 8). For (a), we used the natural breaks (Jens) method to divide the PAs into three levels according to the average trend of NTL. The blue-colored PAs had an average trend of 0, the yellow-colored PAs had a lower than average trend, and the red-colored PAs had a higher than average trend. For (b), the ranking of the average NTL trend of PAs in each continent gave the order as Europe > Asia > North America > South America > Africa > Oceania.       Type Ia PAs represent strict nature reserves. As shown in Figure 8a, the average NTL DN value in buffer zones of type Ia PAs increased in the 0-10 km buffer zone and then decreased with over the 10-100 km buffer zone, reaching the maximum DN value in the 5-10 km buffer zone. The average DN value in the 0-100 km buffer zone of type Ia was much higher than the average DN value within the boundaries of the PAs.

Average NTL Level in Different PAs and Buffers
Type Ib PAs represent wilderness area with large unmodified or slightly modified areas ( Figure  8b). The average DN value within boundaries of type Ib was lower than that of type Ia PAs with the lowest NTL value among all types of PAs. The average DN value of the type Ib buffer increased from the 0-1 km to the 50-100 km buffer zones.
Type II PAs represent national parks (Figure 8c). The change trend of the DN value of the type II PAs buffer was similar to that of Ia (Figure 8a), suggesting a trend of increasing in the close buffer zones (0-10 km) and decreasing in distant buffer zones (10-100 km), with the maximum DN value in the 5-10 km buffer areas. The fluctuation of the NTL DN value of type II was the lowest among all types of PAs.
Type III PAs represent natural monuments or features (Figure 8d). The difference between DN values within PAs and 0-100 km buffer zones were the highest among all types of PAs. The NTL outside of type III PAs increased in the 1-5 km buffer zone and then decreased, indicating more human development in this range.
Type IV PAs represent habitat/species management areas (Figure 8e). The nearest (0-1 km) buffer zone had the lowest DN value. The area with the highest NTL level was concentrated in the 1-10 km buffer zone.
Type V PAs represent protected landscapes/seascapes (Figure 8f). The DN value of type V PAs was much higher than that of the other types of PAs (Figure 9). The average DN value and the range Type Ia PAs represent strict nature reserves. As shown in Figure 8a, the average NTL DN value in buffer zones of type Ia PAs increased in the 0-10 km buffer zone and then decreased with over the 10-100 km buffer zone, reaching the maximum DN value in the 5-10 km buffer zone. The average DN value in the 0-100 km buffer zone of type Ia was much higher than the average DN value within the boundaries of the PAs.
Type Ib PAs represent wilderness area with large unmodified or slightly modified areas (Figure 8b). The average DN value within boundaries of type Ib was lower than that of type Ia PAs with the lowest NTL value among all types of PAs. The average DN value of the type Ib buffer increased from the 0-1 km to the 50-100 km buffer zones.
Type II PAs represent national parks (Figure 8c). The change trend of the DN value of the type II PAs buffer was similar to that of Ia (Figure 8a), suggesting a trend of increasing in the close buffer zones (0-10 km) and decreasing in distant buffer zones (10-100 km), with the maximum DN value in the 5-10 km buffer areas. The fluctuation of the NTL DN value of type II was the lowest among all types of PAs.
Type III PAs represent natural monuments or features (Figure 8d). The difference between DN values within PAs and 0-100 km buffer zones were the highest among all types of PAs. The NTL outside of type III PAs increased in the 1-5 km buffer zone and then decreased, indicating more human development in this range.
Type IV PAs represent habitat/species management areas (Figure 8e). The nearest (0-1 km) buffer zone had the lowest DN value. The area with the highest NTL level was concentrated in the 1-10 km buffer zone.
Type V PAs represent protected landscapes/seascapes (Figure 8f). The DN value of type V PAs was much higher than that of the other types of PAs (Figure 9). The average DN value and the range of the buffer zones was also the highest among all PA types (Figure 9). The brighter areas around the PAs were concentrated within 1-25 km, and the 0-1 km and 25-100 km areas were darker.
Remote Sens. 2019, 11, x FOR PEER REVIEW 13 of 20 of the buffer zones was also the highest among all PA types (Figure 9). The brighter areas around the PAs were concentrated within 1-25 km, and the 0-1 km and 25-100 km areas were darker. The average DN values in type VI PAs buffers were the lowest among all PAs (Figure 8g). The DN values of the 0-1 km buffer zone were the lowest among all types of PAs. As distance from the boundaries of the PAs increased, the DN value in the buffer zones increased, and reached a maximum in the 25-50 km buffer zone. The results showed that the NTL levels within PAs were lower than that of surrounding areas, indicating that human development was limited and controlled by the boundary of the PAs. The NTL level within and outside the PAs increased from 1992 to 2013 (Figure 8). The areas with the highest NTL level in the buffer zone constantly approached the boundary of the type Ia PAs (Figure 8a).
In general, the ranking of the NTL level among the seven types of PAs was V > IV > III > VI > II > Ia > Ib ( Table 2). The ranking of NTL level in the 0-100 km buffer zones outside the boundaries of types of PAs was V > III > IV > Ia > II > Ib > VI (Figure 7). The ranking of the NTL level of buffer zones and the interior was 5-10 km > 1-5 km > 10-25 km > 25-50 km > 50-100 km > 0-1 km. The NTL level within the PAs was significantly lower than that observed outside the PAs (Table 2). For most types of PAs, e.g., Ia, II, III, IV, and V, the brightest areas around the PAs were concentrated in the 1-25 km buffer zone. However, for types Ib and VI, the brightest areas around PAs were concentrated at further distances ( Figure 9).

NTL Growth Rate in Different PAs and Buffers
The 1992-2013 NTL growth rates of every type of PAs and their buffer zones were doubled, except within the 0-1 km for type III PAs (Figure 10), indicating enhanced human development The results showed that the NTL levels within PAs were lower than that of surrounding areas, indicating that human development was limited and controlled by the boundary of the PAs. The NTL level within and outside the PAs increased from 1992 to 2013 (Figure 8). The areas with the highest NTL level in the buffer zone constantly approached the boundary of the type Ia PAs (Figure 8a).
In general, the ranking of the NTL level among the seven types of PAs was V > IV > III > VI > II > Ia > Ib ( Table 2). The ranking of NTL level in the 0-100 km buffer zones outside the boundaries of types of PAs was V > III > IV > Ia > II > Ib > VI (Figure 7). The ranking of the NTL level of buffer zones and the interior was 5-10 km > 1-5 km > 10-25 km > 25-50 km > 50-100 km > 0-1 km. The NTL level within the PAs was significantly lower than that observed outside the PAs (Table 2). For most types of PAs, e.g., Ia, II, III, IV, and V, the brightest areas around the PAs were concentrated in the 1-25 km buffer zone. However, for types Ib and VI, the brightest areas around PAs were concentrated at further distances ( Figure 9).

NTL Growth Rate in Different PAs and Buffers
The 1992-2013 NTL growth rates of every type of PAs and their buffer zones were doubled, except within the 0-1 km for type III PAs (Figure 10), indicating enhanced human development within and around PAs. The NTL growth rate for type Ia PAs was the highest both within and the surrounding outside among all types. Especially, NTL growth rates within the PAs and in the 0-10 km buffer zones were much higher than that of other types. The DN value of the interior Ia increased from 0.04 in 1992 to 0.19 in 2013, representing an increase of 378%. The 0-1 km buffer zone increased by 441% from 0.54 in 1992 to 2.93 in 2013. The growth rates of the 1-5 km and 5-10 km buffer zones were 372% and 274%, respectively. The NTL growth rate from 1992 to 2013 for type Ib PAs decreased from 217% for the internal area to 135% for the furthest buffer zone. The growth rate in the 0-1 km buffer zone was 215%, which was much greater than that of other buffer zones. NTL growth rates for type Ib PAs were much lower than that of Ia. The NTL growth rate for type II PAs ranged from 160% in the 5-10 km buffer zone to 226% in the outermost 50-100 km buffer zone. The growth rate within the type II PAs was 222%, ranking second. For type III PAs, the NTL growth rate was 114% within the protected boundary area. The average growth rate of the buffer zones was the lowest among all types of PAs. The growth rate increased with distance, from 90% in the 0-1 km buffer zone to 170% in the 50-100 km buffer zone. The NTL growth rate within type IV PAs (178%) was slightly higher than that of its surrounding areas (157-175%). The difference in growth rates between buffer zones was small, with a minimum of 157% in the 1-5 km buffer zone and a maximum of 175% in the 5-25 km buffer zone. The NTL growth rate within the type V PAs was 157%. The growth rate of the 0-1 km buffer zone was 164%, which was the highest among all type V PA buffers. The growth rate of the remaining buffer zones decreased with distance from the boundary, i.e., from 133% in the 1-5 km buffer zone to 151% in the 50-100 km buffer zone. From 1992 to 2013, the NTL growth rate of type VI PAs was 245%. The buffers near the boundary had a higher growth rate than that of the buffers farther away. The highest growth rate (263%) occurred in the 0-1 km buffer zone, and the minimum growth rate was 142% in the 50-100 km buffer zone. within and around PAs. The NTL growth rate for type Ia PAs was the highest both within and the surrounding outside among all types. Especially, NTL growth rates within the PAs and in the 0-10 km buffer zones were much higher than that of other types. The DN value of the interior Ia increased from 0.04 in 1992 to 0.19 in 2013, representing an increase of 378%. The 0-1 km buffer zone increased by 441% from 0.54 in 1992 to 2.93 in 2013. The growth rates of the 1-5 km and 5-10 km buffer zones were 372% and 274%, respectively. The NTL growth rate from 1992 to 2013 for type Ib PAs decreased from 217% for the internal area to 135% for the furthest buffer zone. The growth rate in the 0-1 km buffer zone was 215%, which was much greater than that of other buffer zones. NTL growth rates for type Ib PAs were much lower than that of Ia. The NTL growth rate for type II PAs ranged from 160% in the 5-10 km buffer zone to 226% in the outermost 50-100 km buffer zone. The growth rate within the type II PAs was 222%, ranking second. For type III PAs, the NTL growth rate was 114% within the protected boundary area. The average growth rate of the buffer zones was the lowest among all types of PAs. The growth rate increased with distance, from 90% in the 0-1 km buffer zone to 170% in the 50-100 km buffer zone. The NTL growth rate within type IV PAs (178%) was slightly higher than that of its surrounding areas (157-175%). The difference in growth rates between buffer zones was small, with a minimum of 157% in the 1-5 km buffer zone and a maximum of 175% in the 5-25 km buffer zone. The NTL growth rate within the type V PAs was 157%. The growth rate of the 0-1 km buffer zone was 164%, which was the highest among all type V PA buffers. The growth rate of the remaining buffer zones decreased with distance from the boundary, i.e., from 133% in the 1-5 km buffer zone to 151% in the 50-100 km buffer zone. From 1992 to 2013, the NTL growth rate of type VI PAs was 245%. The buffers near the boundary had a higher growth rate than that of the buffers farther away. The highest growth rate (263%) occurred in the 0-1 km buffer zone, and the minimum growth rate was 142% in the 50-100 km buffer zone.

Trends in Different PAs and Buffer Zones
NTL DN values of all PAs and their buffers showed a significant increasing trend from 1992 to 2013 ( Figure 11). Except for the type V PAs, the NTL trends within the other six types of PAs were lower than that in any other buffer zone at the 0-100 km. The change trend within type Ib PAs was the lowest of all the PAs, with a value close to 0. The change trend of type V (0.13 DN/year) was the highest, which was greater than the trend in 50-100 km buffer zone. The average trend of type V PAs buffers was the lowest among all types of PAs.
The trend for type Ia PAs was significantly higher than that of type Ib, indicating that Ia was affected more by human development. The change trends of different types of PAs showed characteristics. The rank of mean values of change trends in buffer zones were V > III > Ia > IV > II >

Trends in Different PAs and Buffer Zones
NTL DN values of all PAs and their buffers showed a significant increasing trend from 1992 to 2013 ( Figure 11). Except for the type V PAs, the NTL trends within the other six types of PAs were lower than that in any other buffer zone at the 0-100 km. The change trend within type Ib PAs was the lowest of all the PAs, with a value close to 0. The change trend of type V (0.13 DN/year) was the highest, which was greater than the trend in 50-100 km buffer zone. The average trend of type V PAs buffers was the lowest among all types of PAs. Ib > VI. The change trend of the type V PAs was the highest, e.g., 0.26/year in 1-5 km and 0.24/year in 5-10 km buffers. The trends in type V, III, Ia, and IV PAs were concentrated in the 1-10 km buffer zone, indicating active urban development. The type VI, Ib, and II PAs had the lower change trends inside PAs, but higher change trends in the peripheral areas far from the PA boundary. The change trends in the area near the PA boundary between 0-10 km were relatively low.

The NTL Distribution Pattern in Different Types of PAs
The primary objective of type Ia PAs is to conserve regionally, nationally, and globally outstanding ecosystems, species (occurrences or aggregations), and/or geodiversity features: these attributes will have been formed mostly or entirely by non-human forces [55]. However, the global NTL level inside type Ia PAs was not the lowest. The NTL value in type Ia PAs was 2.5 times higher than that of the type Ib PAs. Among seven types of PAs, the NTL growth rate inside type Ia was very high from 1992 to 2013, reaching 378%, which was far greater than the others. The NTL for 0-1 km and 1-5 km buffer zones had growth rates that were ranked highly, with values of 441% and 372%, respectively. The area with the highest NTL level around Ia consistently approached the boundary of the PAs. The brightest NTL areas between 1992 and 1999 was in the 25-50 km buffer zone, which shifted to the 5-10 km buffer zone between 2001 and 2003. According to the growth rate and trend of each buffer zone, the brightest NTL area would be encroaching in the 0-5 km area quickly. Human development in the 0-1 km outside PAs could have direct impacts [56,57]. Because the baselines of the nightlight level of different types of protected areas are very different, the actual change values corresponding to the same growth rate are very different even though the growth rates are the same. For example, the NTL growth rate of type Ia PAs exceeded the others, but the absolute increase of the NTL DN value inside and around the Ia PAs was still low in comparison to other types because its baseline of the nightlight level was far less than other types. It is noteworthy that, according to the guidelines, Ia will be degraded or destroyed when subjected to all but very limited human impact [55]. In addition, the area of category Ia was generally small, and human-induced impacts contributed more in small PAs than to stochastic processes [58].
For type Ib PAs, the average NTL DN value was 0.04, which was the lowest among all types. Also, the NTL levels in and out of the PAs boundaries was lower than type Ia PAs. The 50-100 km buffer zone was the area with the most concentrated human population around Ib PAs. The NTL The trend for type Ia PAs was significantly higher than that of type Ib, indicating that Ia was affected more by human development. The change trends of different types of PAs showed characteristics. The rank of mean values of change trends in buffer zones were V > III > Ia > IV > II > Ib > VI. The change trend of the type V PAs was the highest, e.g., 0.26/year in 1-5 km and 0.24/year in 5-10 km buffers.
The trends in type V, III, Ia, and IV PAs were concentrated in the 1-10 km buffer zone, indicating active urban development. The type VI, Ib, and II PAs had the lower change trends inside PAs, but higher change trends in the peripheral areas far from the PA boundary. The change trends in the area near the PA boundary between 0-10 km were relatively low.

The NTL Distribution Pattern in Different Types of PAs
The primary objective of type Ia PAs is to conserve regionally, nationally, and globally outstanding ecosystems, species (occurrences or aggregations), and/or geodiversity features: these attributes will have been formed mostly or entirely by non-human forces [55]. However, the global NTL level inside type Ia PAs was not the lowest. The NTL value in type Ia PAs was 2.5 times higher than that of the type Ib PAs. Among seven types of PAs, the NTL growth rate inside type Ia was very high from 1992 to 2013, reaching 378%, which was far greater than the others. The NTL for 0-1 km and 1-5 km buffer zones had growth rates that were ranked highly, with values of 441% and 372%, respectively. The area with the highest NTL level around Ia consistently approached the boundary of the PAs. The brightest NTL areas between 1992 and 1999 was in the 25-50 km buffer zone, which shifted to the 5-10 km buffer zone between 2001 and 2003. According to the growth rate and trend of each buffer zone, the brightest NTL area would be encroaching in the 0-5 km area quickly. Human development in the 0-1 km outside PAs could have direct impacts [56,57]. Because the baselines of the nightlight level of different types of protected areas are very different, the actual change values corresponding to the same growth rate are very different even though the growth rates are the same. For example, the NTL growth rate of type Ia PAs exceeded the others, but the absolute increase of the NTL DN value inside and around the Ia PAs was still low in comparison to other types because its baseline of the nightlight level was far less than other types. It is noteworthy that, according to the guidelines, Ia will be degraded or destroyed when subjected to all but very limited human impact [55]. In addition, the area of category Ia was generally small, and human-induced impacts contributed more in small PAs than to stochastic processes [58].
For type Ib PAs, the average NTL DN value was 0.04, which was the lowest among all types. Also, the NTL levels in and out of the PAs boundaries was lower than type Ia PAs. The 50-100 km buffer zone was the area with the most concentrated human population around Ib PAs. The NTL distribution for type II PAs were similar to type VI and both were lower among the types of PAs. The NTL levels of types II and VI were not significantly higher than those of Ia and Ib.
The main function of type IV PAs is to protect certain species and habitats. The NTL levels inside type IV PAs and buffer zones were significantly higher than all other types except V. NTL levels within and outside type III PAs were high, second only to type V. NTL levels within and outside of type V PAs were much higher than other types.
The result suggested that the NTL levels were significantly lower within the PAs than that of the surrounding 0-100 km buffer zones. In particular, there was a big difference between inside and outside the boundary of types Ia and III PAs. For most PAs, the surrounding areas with the highest DN value were in the 1-25 km buffer zone.

Skyglow as a Biodiversity Threat
Skyglow occurs when artificial light is scattered by water, atmospheric molecules, or aerosols and returned to Earth. In this study, we used the NTL data to show the impact of urbanization on protected areas. However, we did not apply any light propagation models to integrate the skyglow effect due to indirect lights. Skyglow is of increasing concern since it is able to multiply and extend the light pollution effect and affect those areas with no direct light pollution. Because of skyglow, the biological impacts of artificial light are not limited to the vicinity of the light source and may spread over much larger extents via several mechanisms. Individual lights may be visible kilometers away from their source, and the addition of artificial skyglow can extinguish such lunar light cycles and permanently remove dark nights from a landscape [59,60]. Therefore, artificial lighting may have an effect on natural ecosystems even when the source of lighting lies kilometers away. Most organisms have evolved molecular circadian clocks controlled by natural day-night cycles. These clocks play key roles in metabolism, growth, and behavior [61]. As the world grows ever-more illuminated, many light-sensitive species will be lost, especially in or near highly illuminated urban areas [62]. Light pollution threatens biodiversity through changed night habits, including organismal movements [63,64]; foraging [65][66][67][68]; interspecific interactions [69]; communication [70,71]; reproduction [72] of insects, amphibians, fish, birds, bats, and other animals; and it can disrupt plants resulting in phenological changes by distorting their natural day-night cycle [63].
Therefore, measures are necessary to prevent or reduce the ecological impact of night-time light pollution. Maintaining and increasing natural unlit areas is likely to be the most effective option for reducing the ecological effects of lighting [73]. From the lighting differences inside and outside the PAs, it can be seen that in the process of human development, PAs have greatly reduced human interference, so PAs are undoubtedly the best place to maintain darkness. More stringent control measures should be implemented within and around PAs, such as limiting the duration of lighting, reducing the trespass of lighting, changing the intensity of lighting, and changing the spectrum of lighting [73,74].
Setting an area surrounding PAs is an ideal option to provide a buffer against the light pollution released by human development. However, it cannot be solved by using only remote sensing data. Different types of PAs are in different natural and socio-economic conditions, and as such, the appropriate buffer radius should also be different. The success of planners in reducing the ecological impacts of light pollution will ultimately depend on an assessment of the critical mechanisms and thresholds that determine those impacts in a particular environment [73]. One possible way is to combine remote sensing data with biodiversity data to explore how biodiversity in protected areas with similar natural conditions respond to different lighting patterns. This may be of great significance for solving light pollution near the PAs and may also contribute to determining how much of a buffer distance should be set around the PAs.

Limitations of Night-Time Lighting Data
The 1992-2013 NTL data used in this study were derived from six different sensors. Therefore, it is difficult to discriminate whether the small difference between the light data of different years was caused by the sensor or by changes in the field light. Moreover, we are not sure how much data error was caused by the sensor. In this study, we compared the nightlight levels between 1992 and 2013. The time span and the difference in the NTL data were sufficiently large, such that the data error caused by the different sensor could be ignored.
In addition, the light data could not detect negative light growth (i.e., reductions in light) in each PA. To reduce the error caused by the sensor when processing the light data, we assumed that all regions were developing positively, and the DN value was increasing; thus, we assigned the larger DN value of the previous year to the darker pixels corresponding to the following year. If the actual nightlight of a certain area was dark each year, the DN value would remain unchanged. On the global scale, the nightlight within and around each type of PA was increasing each year, and there was no negative growth in the brightness of the PA. However, when assessing individual PAs, there may be negative growth due to increasingly strict regulations and improved awareness of protection, e.g., lighting tools could be replaced by those with a lower brightness and lighting time could be reduced. Therefore, the lighting data may not be suitable for the study of PAs on a small spatial scale, especially in developed countries where the management of PAs is more stringent.