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

Abstract

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.

1. Introduction

A protected area (PA) is defined as a geographical space that is recognized, dedicated, and managed, through legal or other effective means, to achieve the long-term conservation of nature with associated ecosystem services and cultural values in mind [1]. In general, protected areas include national parks, national forests, national seashores, all levels of natural reserves, wildlife refuges and sanctuaries, and designated areas for the conservation of native biological diversity and natural and cultural heritage and significance. Protected areas also include some of the last frontiers that have unique landscape characteristics and ecosystem functions in wilderness conditions [2].

PAs reflect the efforts to protect the world’s threatened species and their habitats. PAs are increasingly recognized as essential providers of ecosystem services and biological resources, key components in climate change mitigation strategies, and in some cases, are vehicles for protecting threatened human communities or sites of great cultural and spiritual value.

Humans have created protected areas over the past millennia for a multitude of reasons [3,4]. The establishment of Yellowstone National Park in 1872 by the US Congress ushered in the modern era of governmental protection of natural areas, which catalyzed a global movement [5,6]. However, even with the implementation of a tremendous variety of monitoring programs, as well as conservation planning efforts and achievements, species’ population decline, biodiversity loss, extinction, system degradation, pathogen spread, and state change events are occurring at unprecedented rates [7,8]. The effects are augmented by continued changes in land use, urbanization, and invasive spread alongside the direct, indirect, and interactive effects of climate change and disruption [3,4]. Protected areas become more important in serving as indicators of ecosystem conditions and functions, either by their status and/or by comparison with unprotected adjacent areas. Protected areas are highly prized by society with diversified representative characteristics. Earth’s remaining wilderness areas are becoming increasingly important buffers against the changing environment. However, they are not yet an explicit target in international policy frameworks [9].

Over the past four decades, more PAs have been established around the world. PAs play a vital role in conserving biodiversity; specifically, PAs provide a paradise for endangered and species from declining habitats and a sanctuary for over-harvested and poached species [10,11,12]. PAs have a positive impact on the local environment, such as maintaining water resources, regulating climate, and preventing forest damage [13], among other benefits [14,15].

On the other hand, ecosystems of PAs have been disturbed by human development, in particular, urbanization, such that biodiversity and habitats in those PAs have been reduced [16,17,18,19,20,21]. The effectiveness of PAs has been affected by the land use of the surrounding areas. Nightlight may impose problems to PAs [22,23]. The threat from human development to PAs is a concern. Currently, assessment studies on the impact of human development on PAs are mainly based on data, such as human population [14] and housing [24] at regional or local scales. However, accurate, reliable, and comprehensive population or housing statistics are often not available in global scale monitoring.

The Defense Meteorological Satellite Program/Operational Linescan System (DMSP/OLS) night-time stable light (NTL) dataset is a timely and effective data source for monitoring human development and mapping the dynamics of land cover on a regional or global scale [25,26,27,28]. The NTL data have a wide time span and are suitable for dynamic analysis on long-term sequences. The NTL can be measured and used as a proxy of human development. NTL imposes impacts on several taxa in terrestrial and aquatic ecosystems, including mammals, birds, reptiles, amphibians, fish, invertebrates, and plants. It is necessary to study the lighting conditions of reserves [29,30,31]. DMSP/OLS NTL is an excellent data source for analyzing the impact of human development on various ecosystems [32], vegetation [33], and biodiversity [34] over long-term sequences and over large geographical areas and in remote locations [35].

DMSP/OLS NTL data have been used as an effective indicator to evaluate the conservation of protected areas [36,37,38]. Reported studies have investigated the difference in NTL between the interior and the surrounding areas and have revealed that the NTL level of the PAs is lower than that of the 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², 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.

2. Materials and Methods

2.1. 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. Seven types of terrestrial protected areas defined by IUCN.

IUCN PAs Type	IUCN PAs Code	Number of PAs	Area (km2)	Description
Strict Nature Reserve	Ia	11,921	3,874,328	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	Ⅱ	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	Ⅲ	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.
Habitat/Species Management Area	Ⅳ	56,836	7,250,215	PAs aiming to protect particular species or habitats and management reflects this priority. Many category IV PAs will need regular, active interventions to address the requirements of particular species or to maintain habitats, but this is not a requirement of the category.
Protected Landscape/Seascape	V	44,915	4,351,409	PAs where the interaction of people and nature over time have produced the area of distinct character with significant ecological, biological, cultural, and scenic value.
Protected Area with Sustainable use of Natural Resources	Ⅵ	7177	12,757,697	PAs that conserve ecosystems and habitats, together with associated cultural values and traditional natural resource management systems. They are generally large, with most of the area in a natural condition, where a proportion is under sustainable natural resource management and where low-level, non-industrial use of natural resources compatible with nature conservation is seen as one of the main aims of the area.

. 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].

2.2. 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].

2.3. 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.

2.3.1. 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 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]:

$$D{N}_{\left(n,i\right)}=\{\begin{array}{c}0, D{N}_{\left(n+1 ,i\right)}=0 \\ D{N}_{\left(n-1 ,i\right)}, D{N}_{\left(n+1 ,i\right)}>0 nd D{N}_{\left(n-1 ,i\right)}>D D{N}_{\left(n,i\right)}\\ D{N}_{\left(n,i\right)}, otherwise \end{array}(n=1992,1993,\dots ,2013),$$

(1)

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.

2.3.2. 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:

$$\mathrm{F}12=0.870\times \mathrm{F}10-0.078({\mathrm{R}}^2=0.925)$$

(2)

Correct F14 based on F12:

$$\mathrm{F}14=0.927\times \mathrm{F}12-1.709({\mathrm{R}}^2=0.965)$$

(3)

Correct F15 based on F14:

$$\mathrm{F}15=0.961\times \mathrm{F}14+1.408({\mathrm{R}}^2=0.951)$$

(4)

Correct F16 based on F15:

$$\mathrm{F}16=1.062\times \mathrm{F}15+1.672({\mathrm{R}}^2=0.960)$$

(5)

2.3.3. 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].

$$D{N}_{\left(n,i\right)}=\{\begin{array}{c}0, D{N}_{\left(n,i\right)}^a=0|D{N}_{\left(n,i\right)}^b=0\\ D{N}_{\left(n,i\right)}^a+D{N}_{\left(n,i\right)}^b, \mathrm{otherwise} \end{array}(n=1992,1993,\dots ,2013),$$

(6)

where $D{N}_{\left(n,i\right)}^a$ and $D{N}_{\left(n,i\right)}^b$ 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.

2.3.4. 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.

2.3.5. 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.

2.4. 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:

$$Q=median\frac{X_j-X_i}{j-i} 0<i<j\le n$$

(7)

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.

$$Q=\{\begin{array}{c}{Q}_{\left[\left(N-1\right)/2\right] } N is odd\\ \left(Q_{\left(N/2\right)}+Q_{\left[\left(N+2\right)/2\right]}\right)/2 N is even\end{array}$$

(8)

where Q > 0, Q < 0 and Q = 0 indicate that there is a rising trend, a downward trend, and no obvious trend, respectively.

3. Results

3.1. 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).

3.2. Average NTL Level in Different PAs and Buffers

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.

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).

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.

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. NTL level of interior and buffer zones for the seven types of PAs.

PA Type	Interior	0–1 km	1–5 km	5–10 km	10–25 km	25–50 km	50–100 km
Ia	0.098	1.572	2.372	2.598	2.478	2.344	1.916↑
Ib	0.039	0.496	0.920	1.380	1.897	2.317	2.452↑
II	0.147	1.129	1.602	1.920	1.932	1.702	1.427
III	0.594	4.868↑	6.289↑	6.038	4.734↑	3.272↑	2.013↑
IV	0.652	1.317	3.593↑	3.301	2.631	2.078	1.456
V	3.184↑	3.270↑	7.151↑	6.840↑	5.418↑	3.610↑	2.357↑
VI	0.160	0.169	1.217	1.483	1.739	1.811	1.692
Ave.	0.696	1.832	3.307	3.366	2.975	2.448	1.902

). 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).

3.3. 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.

3.4. 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 > 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.

4. Discussion

4.1. 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.

4.2. 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.

4.3. 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.