**1. Introduction**

A protected area (PA) is defined as a geographical space that is recognized, dedicated, and managed, through legal or other e ffective 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 e fforts 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 grea<sup>t</sup> 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 e fforts 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 e ffects 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 bu ffers against the changing environment. However, they are not ye<sup>t</sup> 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–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–21]. The effectiveness of PAs has been a ffected 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 e ffective data source for monitoring human development and mapping the dynamics of land cover on a regional or global scale [25–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–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 e ffective indicator to evaluate the conservation of protected areas [36–38]. Reported studies have investigated the di fference 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 km2, 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. 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].


**Table 1.** Seven types of terrestrial protected areas defined by IUCN.


**Table 1.** *Cont.*

#### *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]:

$$DN\_{(\tiny{u,j})} = \begin{cases} 0, & DN\_{(\tiny{u+1,j})} = 0 \\ DN\_{(\tiny{u-1,j})'} & DN\_{(\tiny{u+1,j})} > 0 \text{ nd } DN\_{(\tiny{u-1,j})} > D \, DN\_{(\tiny{u,j})} \text{ (n = 1992, 1993, \dots, 2013), (1)} \\ & DN\_{(\tiny{u,j})'} & \text{otherwise} \end{cases}$$

where *DN*(*<sup>n</sup>*−1, *i*), *DN*(*<sup>n</sup>*, *i*), and *DN*(*n*<sup>+</sup>1, *i*) are the DN values of the *i*th lit pixel of the NTL data in the (*n* − 1)th, the *n*th, 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:

$$\text{F12} = 0.870 \times \text{F10} - 0.078 \left(\text{R}^2 \, \, \, = \, 0.925\right) \tag{2}$$

Correct F14 based on F12:

$$\text{F14} = 0.927 \times \text{F12} - 1.709 \left(\text{R}^2 \, \, \, = \, 0.965\right) \tag{3}$$

Correct F15 based on F14:

$$\text{F15} = 0.961 \times \text{F14} + 1.408 \left(\text{R}^2 = 0.951\right) \tag{4}$$

Correct F16 based on F15:

$$\text{F16} = 1.062 \times \text{F15} + 1.672 \left( \text{R}^2 \, \text{ } = \, 0.960 \right) \tag{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].

$$DN\_{(\tiny\rm (n,j)} = \begin{cases} 0, & DN\_{(\tiny\rm (n,j)}^a = 0 \middle| DN\_{(\tiny\rm (n,j)}^b = 0 \\\ \cdot DN\_{(\tiny\rm (n,j)}^a + \cdot DN\_{(\tiny\rm (n,j)}^b) & \text{otherwise} \end{cases} \quad \text{(\$n = 1992, 1993, ..., 2013)} \text{,} \tag{6}$$

where *DNa* (*<sup>n</sup>*,*<sup>i</sup>*) and *DN<sup>b</sup>* (*<sup>n</sup>*,*<sup>i</sup>*) are the DN values of the *i*th lit pixel from two NTL data in the *n*th year, and *DN*(*<sup>n</sup>*,*i*) is the DN value of the *i*th lit pixel of the intra-annual composite in the *n*th year.

#### 2.3.4. Correction for Di fferent 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–48] and have used gross domestic product (GDP) and electricity consumption (EC) to evaluate the performance of intercalibration techniques [28,49–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.

**Figure 1.** Gross domestic product (GDP; unit: trillions of dollars) and global the sum of NTL (SNTL) (1992–2013) from (**a**) the raw NTL and (**b**) the calibrated NTL. Electricity consumption (EC; unit: billion kWh) and global SNTL (1992–2013) from (**c**) the raw NTL and (**d**) the calibrated NTL. DMSP/OLS—Defense Meteorological Satellite Program/ Operational Linescan System.

**Figure 2.** The correlation of GDP (**a**) and EC (**b**) with the raw NTL and the calibrated NTL. The black line is the 1:1 line.

#### *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 = \operatorname{median} \frac{X\_j - X\_i}{j - i} \qquad \qquad 0 < i < j \le n \tag{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{cases} \begin{array}{l} Q\_{\left[\left(N-1\right)/2\right]} \\ Q\_{\left(N/2\right)} + Q\_{\left[\left(N+2\right)/2\right]} \end{array} & N \text{ is odd} \\ \end{cases} \tag{8}$$

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