Post-Earthquake Night-Time Light Piecewise (PNLP) Pattern Based on NPP / VIIRS Night-Time Light Data: A Case Study of the 2015 Nepal Earthquake

: Earthquakes are unpredictable and potentially destructive natural disasters that take a long time to recover from. Monitoring post-earthquake human activity (HA) is of great signiﬁcance to recovery and reconstruction work. There is a strong correlation between night-time light (NTL) and HA, which aid in the study of spatiotemporal changes in post-earthquake human activities. However, seasonal and noise impact from National Polar-Orbiting Partnership Satellite Visible Infrared Imaging Radiometer Suite (NPP / VIIRS) data greatly limits their application. To tackle these issues, random noise and seasonal ﬂuctuation of NPP / VIIRS from January 2014 to December 2018 is removed by adopting the seasonal-trend decomposition procedure based on loess (STL). Based on the theory of post-earthquake recovery model, a post-earthquake night-time light piecewise (PNLP) pattern is explored by employing the National Polar-Orbiting Partnership Satellite Visible Infrared Imaging Radiometer Suite (NPP / VIIRS) monthly data. PNLP indicators, including pre-earthquake development rate (k p ), recovery rate (k r1 ), reconstruction rate (k r2 ), development rate (k d ), relative reconstruction rate (k rp ) and loss (S), are deﬁned to describe the PNLP pattern. Furthermore, the 2015 Nepal earthquake is chosen as a case study and the spatiotemporal changes in di ﬀ erent areas are analyzed. The results reveal that: (1) STL is an e ﬀ ective algorithm for obtaining HA trend from the time series of denoising NTL; (2) the PNLP pattern, divided into four phases, namely the emergency phase (EP), recovery phase (RP-1), reconstruction phase (RP-2), and development phase (DP), aptly describes the variation in post-earthquake HA; (3) PNLP indicators are capable of evaluating the recovery di ﬀ erences across regions. The main socio-economic factors a ﬀ ecting the PNLP pattern and PNLP indicators are energy source for lighting, type of building, agricultural economy, and human poverty index. Based on the NPP / VIIRS data, the PNLP pattern can reﬂect the periodical changes of HA after earthquakes and provide an e ﬀ ective means for the analysis and evaluation of post-earthquake recovery and reconstruction.


Introduction
Earthquake is a natural disaster that could cause great damage to human society and its effects linger for a long time [1]. Due to its sudden onset and unpredictability, the damage caused by an earthquake to the environment and social development is imminent. Timely and accurate understanding of post-earthquake recovery and reconstruction is of great significance in disaster plan deployment, reconstruction policy formulation, and capital investment.
Post-earthquake recovery and reconstruction is a complicated and lengthy process involving natural recovery [2], reconstruction of buildings and structures [3], and psychological counselling [4].

NPP/VIIRS Monthly Composite Data
Currently, Defence Meteorological Satellite Program-Operational Linescan System (DMSP/OLS) and NPP/VIIRS are the two kinds of night-time light data that are widely utilized. DMSP/OLS images have been available since the early 1970s, but they have certain shortcomings, such as insufficient radiation accuracy, low spatial resolution, and oversaturation [41,42]; moreover, only OLS annual data are freely available. Therefore, OLS data have certain limitations in studying the time series changes in HA intensity. In March 2015, the National Geophysical Data Center (NGDC) of the National Oceanic and Atmospheric Administration (NOAA) released the VIIRS data product for several months (http://ngdc.noaa.gov/eog/vii-rs/download_monthly.html). Compared with

NPP/VIIRS Monthly Composite Data
Currently, Defence Meteorological Satellite Program-Operational Linescan System (DMSP/OLS) and NPP/VIIRS are the two kinds of night-time light data that are widely utilized. DMSP/OLS images have been available since the early 1970s, but they have certain shortcomings, such as insufficient radiation accuracy, low spatial resolution, and oversaturation [41,42]; moreover, only OLS annual data are freely available. Therefore, OLS data have certain limitations in studying the time series changes in HA intensity. In March 2015, the National Geophysical Data Center (NGDC) of the National Oceanic and Atmospheric Administration (NOAA) released the VIIRS data product for several months (http://ngdc.noaa.gov/eog/vii-rs/download_monthly.html). Compared with DMSP/OLS, NPP/VIIRS Remote Sens. 2020, 12,2009 4 of 21 data have more beneficial characteristics. The spatial resolution of NPP/VIIRS is about 500 m and the data are irradiated and not oversaturated [43,44]. Since April 2012, its monthly data has been freely available; thus, the VIIRS data have a higher potential than the OLS data in analyzing and evaluating studies related to long time series of HA.
In this paper, in order to analyze the variation of HA intensity, the NPP/VIIRS monthly data from January 2014 to December 2018 are utilized.

Auxiliary Data
Data on Nepal's administrative boundaries (GADM data) at the county level were obtained from the Global Administrative Areas (http://www.gadm.org/). Population data from the European Union Global Human Settlement (https://ghsl.jrc.ec.europa.eu/download.php?ds=pop) were used to remove the background noise. Monthly NDVI data were download from Modis/Terra Vegetation Indices Monthly L3 0.05Deg CMG (MOD13C2) collection (https://ladsweb.modaps.eosdis.nasa.gov/). Statistical data on house damage caused by the earthquake (indicating the degree of the damage in different areas) were acquired from the Ministry of Home Affairs (MHA), Government of Nepal. The socio-economic statistical data were obtained from Central Bureau of Statistics (CBS), 2011 (http://cbs.gov.np).

Methods
This paper includes the following steps: (1) data collection and (2) pre-processing; (3) construction and calculation of the PNLP pattern; (4) the analysis of feature importance of PNLP indicators based on random forest regressor (RFR). Figure 2 shows a detailed process of the method.
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 23 DMSP/OLS, NPP/VIIRS data have more beneficial characteristics. The spatial resolution of NPP/VIIRS is about 500 m and the data are irradiated and not oversaturated [43,44]. Since April 2012, its monthly data has been freely available; thus, the VIIRS data have a higher potential than the OLS data in analyzing and evaluating studies related to long time series of HA.
In this paper, in order to analyze the variation of HA intensity, the NPP/VIIRS monthly data from January 2014 to December 2018 are utilized.

Auxiliary Data
Data on Nepal's administrative boundaries (GADM data) at the county level were obtained from the Global Administrative Areas (http://www.gadm.org/). Population data from the European Union Global Human Settlement (https://ghsl.jrc.ec.europa.eu/download.php?ds=pop) were used to remove the background noise. Monthly NDVI data were download from Modis/Terra Vegetation Indices Monthly L3 0.05Deg CMG (MOD13C2) collection (https://ladsweb.modaps.eosdis.nasa.gov/). Statistical data on house damage caused by the earthquake (indicating the degree of the damage in different areas) were acquired from the Ministry of Home Affairs (MHA), Government of Nepal. The socio-economic statistical data were obtained from Central Bureau of Statistics (CBS), 2011 (http://cbs.gov.np).

Methods
This paper includes the following steps: (1) data collection and (2) pre-processing; (3) construction and calculation of the PNLP pattern; (4) the analysis of feature importance of PNLP indicators based on random forest regressor (RFR). Figure 2 shows a detailed process of the method.

Data Pre-Processing
Raw NPP/VIIRS monthly data cannot be directly applied. It contains negative light intensity values caused by sensor noise and calibration noise during data synthesis, as well as outliers caused by night fire flares and oil and gas flares [11,12]. To remove the background noise, Shi et al. assumed that DMSP/OLS and NPP/VIIRS have the same illumination region; they adopted DMSP/OLS data as a mask to correct the NPP/VIIRS dataset [44,45]. NPP/VIIRS corrected by this method can be better applied to population, urbanization, and other aspects, but two problems still persist after data correction. Firstly, their spatial resolution is different, leading to problems in adopting DMSP/OLS data correction. Secondly, there are still seasonal variations in the corrected data. To avoid the above problems, the Global Human Settlement Layer (GHSL) population data was utilized as mask data to remove the background noise. On the other hand, since 2017 the way NASA corrected the data for airglow has changed, which clearly increased the trend [46]. Hence, we first removed background noise and outliers. The main steps were as follows: 1.
The study area is clipped based on GADM data.
The normalized GHSL data are divided into several categories at 0.2 intervals. 4.
Samples of different population densities are generated based on random sampling and the mean value of DN is counted under different population densities. 5.
The monthly background noise threshold is chosen according to mean value of DN. NPP/VIIRS data below threshold are replaced by zero in order to remove the background noise. 6.
In order to eliminate the outliers, the maximum value of Kathmandu is taken as the maximum value of the study area and a value greater than the DN value in the study area is replaced by zero.

PNLP Conceptual Pattern
The process of disaster recovery is extremely complex and uncertain and it is difficult to construct the comprehensive framework of disaster recovery [48]. Haas et al. [49] provided a framework for disaster recovery. In this framework, a community experiences four post-disaster stages in a regular and predictable order. The first phase, called the emergency stage, is about several weeks after the disaster. The second phase, called the restoration stage, is about one month to six months after the disaster. The third and fourth stages are called the Reconstruction-1 and Reconstruction-2, respectively.
The intensity of HA in each phase has different characteristics. Based on the variation in intensity of HA and the theory of post-earthquake recovery, this paper constructs a PNLP pattern ( Figure 3). Most studies have shown a strong correlation between changes in night-time light and HA intensity factors, such as population and GDP. We assumed that the intensity of HA varied linearly over time. Hence, the PNLP pattern is ordered and predictable linear sequence is as defined by Haas.
At a certain period before the earthquake, according to the regional development situation, the intensity of HA increased or decreased linearly. The first post-earthquake phase is defined as the emergency phase (EP), which takes a relatively short period of time, with a sudden increase in the intensity of HA. The second phase after the earthquake is defined as the recovery phase (RP-1), in which emergency relief measures are gradually withdrawn, resettlement services gradually stabilize, and the intensity of HA gradually decreases. The third stage is defined as the reconstruction stage Figure 3. PNLP pattern. T1 represents the month of the earthquake; EP, RP-1, RP-2, and DP represent the emergency phase, recovery phase, reconstruction phase, and development phase, respectively; the defined PNLP indicators include k p , k r1 , k r2 , k d , and S. At a certain period before the earthquake, according to the regional development situation, the intensity of HA increased or decreased linearly. The first post-earthquake phase is defined as the emergency phase (EP), which takes a relatively short period of time, with a sudden increase in the intensity of HA. The second phase after the earthquake is defined as the recovery phase (RP-1), in which emergency relief measures are gradually withdrawn, resettlement services gradually stabilize, and the intensity of HA gradually decreases. The third stage is defined as the reconstruction stage (RP-2). In this stage, a series of work, such as life reconstruction, is carried out and the intensity of HA gradually increases. The fourth stage is defined as the development phase (DP), which may show different characteristics, that is, it may continue to develop on the basis of the reconstruction stage, or the intensity of HA enters a decline stage at a certain period after reconstruction. PNLP pattern ( Figure 3) describes various post-earthquake stages. The data we employed in this study are monthly data, however, as we know, the EP takes 24 h to 2 weeks [50]. Therefore, this pattern may not detect the EP changes, so EP and RP-1 based on monthly data are defined as phases (from the earthquake month to the lowest HA level).
The speed and the quality of recovery process must be monitored [51]. From Figure 3, we know that the PNLP pattern was divided into two parts: pre-earthquake and post-earthquake. We developed the K p (pre-earthquake development rate), which reflects the changing levels of HA intensity before the earthquake. K r1 (recovery rate) reflects the changing levels of HA intensity in RP-1. K r2 (reconstruction rate) reflects the changing levels of HA intensity in RP-2. K d (development rate) reflects the changing levels of HA intensity in DP. K rp (relative reconstruction rate) is a relative quantity, and it reflects the changing levels of HA intensity relative to the pre-earthquake HA level. S (loss) reflects the amount of HA loss during the EP and RP-1. In order to better analyze the changes in HA at various stages, we defined the parameters in Table 1. Table 1. PNLP pattern parameter definition.

Parameter Descriptions
K p Reflects the changing levels of human activity intensity before the earthquake K r1 Reflects the changing levels of human activity intensity in RP-1 K r2 Reflects the changing levels of human activity intensity in RP-2 K d Reflects the changing levels of human activity intensity in DP K rp Reflects the changing levels of HA intensity relative to the pre-earthquake level S Reflects the amount of HA loss during the EP and RP-1

HA Intensity Indicator
The total night-time light (TNL) value has a significant correlation with the economics of the region, reflecting the level of HA in the region [45], as shown below: where DN i represents the value of i-th pixel. TNL can intuitively reflect the socio-economic situation of a region. The average light index (ALI) is utilized to compare the differences in HA in different regions [52], as shown in Formula (2): where N is the number of pixels with the DN value greater than 0 in a region. In this study, ALI values from January 2014 to December 2018 are adopted to represent the short-term dynamic trend of HA. However, current studies prove that the monthly composite data of NPP/VIIRS have seasonal variations due to the influence of snow caver, NDVI, and albedo; hence, there are some problems in the continuous monitoring research.

ALI Trend Extraction Based on STL Algorithm
STL algorithm is widely used in remote sensing, especially in monitoring NDVI changes [53,54]. In NTL, Liu et al. [55] successfully monitored two religious events in Saudi Arabia and Iraq based on the STL algorithm. STL algorithm divides the original time series data into seasonal component, trend component, and remainder component [38], as shown in Formula (3): where Y v represents the origin time series data, and T v , S v , R v represent the trend component, seasonal component, and remainder component, respectively. The main influencing factors of continuous monitoring of NTL are inter-annual seasonal variation and noise, which make time series decomposition algorithm ideal in the study of NTL.
To remove the seasonal differences and noise in ALI sequences, the STL algorithm was employed. In this study, the STL algorithm is applied to decompose the ALI value and the trend component is taken as the changing trend of HA. The algorithm was implemented in R 3.6.3.

Solution of PNLP Parameters
We assumed that the stages of the PNLP pattern were linear. After the trend component was obtained by the STL algorithm, we determined the separation points of each stage through the extremum and then constructed the piecewise linear model through the least square method (LSM).
Based on the time interval between the separation points, the time parameters (T 1 T 2 and T 2 T 3 ) were solved. The slope of each stage was obtained through LSM, the parameters of k p , k r1 , k r2 , and k d , were obtained.
The calculation formula of k rp is as follows: The calculation formula of S is as follows: In this study, we assumed that the PNLP pattern was linear and the EP and RP-1 were merged, so the simplified formula was as follows:

Ineffective Factors of PNLP Indicators Based on Ramdom Forest Regression (RFR)
In order to determine which socio-economic indicators had a great influence on the PNLP indicators, we chose the random forest regression (RFR) algorithm to fit the selected socio-economic indicators and PNLP indicators. RF is an ensemble of decision tree based on classification and regression tree (CART) [56].
In this study, for the evaluation of feature importance, we utilized the mean decrease accuracy importance score. The main steps were as follows:

1.
For each decision tree, the prediction error (err1) of out-of-bag (OOB) is recorded. For a predictor variable, its importance is calculated as the mean of the difference between the transformed prediction error and the original. The formula is as follows: where importance represents the importance of predictor variable, N represents the number of trees.

ALI Sequence Test Based on Mann-Kendall Mutation Test
In order to test whether there was a mutation point in the original ALI sequence, the Mann-Kendall (MK) mutation test was adopted. For a time series (x 1 , x 2 , . . . , x n ), construct an order column: where under the assumption that the time series is random and independent, define the statistics: where UF 1 = 0; E(A k ) and Var (A k ) are the mean and variance of the A k , respectively. In the formula: UF k is a standard normal distribution. Arrange the time series in reverse order, and then calculate according to the above formula, while making: By analyzing the statistical sequences UF k and UB k , the time of mutation can be clarified. If the intersection points of the two curves UF k and UB k appear, and the intersection point is between the critical straight lines, then the moment corresponding to the intersection is the moment when the mutation starts [57].

The Level of Background Noise
Utilizing GHSL population data as auxiliary data, the mean DN value of each month is obtained at different population density (Figure 4).
From Figure 4, the mean of DN value can be observed in different population density. We know that as population density increased, so did the DN value. When we compare the DN values between different months, we can see significant differences. However, for one month, there was no obvious difference in the mean value between no population (=0) and low population (0-0.2). This may be due to the inability of NTL to monitor low HA areas [58]. As the population density increased, the mean value increased. In this study, low population density was not considered. Therefore, we utilized the mean value of 0.2-0.4 population density level as the threshold to distinguish background noise from HA. Background noise trend is shown in Figure 5.
Utilizing GHSL population data as auxiliary data, the mean DN value of each month is obtained at different population density (Figure 4).
From Figure 4, the mean of DN value can be observed in different population density. We know that as population density increased, so did the DN value. When we compare the DN values between different months, we can see significant differences. However, for one month, there was no obvious difference in the mean value between no population (=0) and low population (0-0.2). This may be due to the inability of NTL to monitor low HA areas [58]. As the population density increased, the mean value increased. In this study, low population density was not considered. Therefore, we utilized the mean value of 0.2-0.4 population density level as the threshold to distinguish background noise from HA. Background noise trend is shown in Figure 5.  From Figure 5, we know that the trend of background noise is dynamic. Moreover, it was clear that since 2017, the original noise level had increased significantly. We can see that after we subtracted the step factor, the corrected noise level was basically the same. In order to compare the denoising effect, the images of January 2015 and January 2018 were selected, as shown in Figure 6. From Figure 5, we know that the trend of background noise is dynamic. Moreover, it was clear that since 2017, the original noise level had increased significantly. We can see that after we subtracted the step factor, the corrected noise level was basically the same. In order to compare the denoising effect, the images of January 2015 and January 2018 were selected, as shown in Figure 6. From Figure 5, we know that the trend of background noise is dynamic. Moreover, it was clear that since 2017, the original noise level had increased significantly. We can see that after we subtracted the step factor, the corrected noise level was basically the same. In order to compare the denoising effect, the images of January 2015 and January 2018 were selected, as shown in Figure 6. As can be seen from Figure 6, compared with original images, the corrected images removed some noise area and we can see the lit area clearly. From the area in the black rectangular box, we can observe that in the original image, the human illumination and background noise were mixed, As can be seen from Figure 6, compared with original images, the corrected images removed some noise area and we can see the lit area clearly. From the area in the black rectangular box, we can observe that in the original image, the human illumination and background noise were mixed, but the human illumination can be clearly seen in the corrected image. The results showed our denoising method could remove some noises and eliminate the step change in the background noise.

Trend Analysis of STL Algorithm
The original ALI curve has obvious fluctuations and there are certain limitations in studying the dynamic trend of HA (Figure 6c,e,g and Figure 7a). The ALI trend acquired from the STL algorithm clearly showed the change trend of HA (Figure 6d From the ALI trend, it was observed that the intensity of HA showed an upward or downward trend during a certain period of time before the earthquake and then the intensity of HA first declined and then increased due to the impact of the earthquake. Furthermore, some regions showed a downward trend after an initial upward trend, but others showed an upward trend throughout the study period. In order to better describe the impact of human activity intensity caused by the Nepal earthquake on different regions, we developed the PNLP indicators to quantitatively analyze the similarities and differences in HA in various areas of Nepal.   NTL is highly correlated with HA and earthquakes can cause changes in HA intensity [34]. We analyzed the intensity of HA in different areas of Nepal by studying the trend in light changes. Figure  8 shows the PNLP indicators of the Nepal earthquake in different areas.  NTL is highly correlated with HA and earthquakes can cause changes in HA intensity [34]. We analyzed the intensity of HA in different areas of Nepal by studying the trend in light changes. Figure 8 shows the PNLP indicators of the Nepal earthquake in different areas. Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 23 Figure 8. PNLP indicators. The red dotted line is an auxiliary line to analyze regional differences of PNLP indicators. Kp was in a relatively stable range (from −0.01 to 0.01) in most areas. The fastest developing area before the earthquake was Bhaktapur (kp = 0.009), while Rasuwa, Parsa, and Dolakha were the areas with the slowest development (kp < -0.01).
In RP-1, recovery speed was slow in all regions, with kr1 ranging from 0 to 0.03. Kathmandu and Lalitpur had the highest recovery rate, with kr1 around 0.025, followed by regions Sarlahi and Lamjung, with kr1 greater than 0.01. In some remote areas kr1 was relatively small, which may have been due to the slow pace of Nepal's rescue of remote areas and the government services that cannot meet regional rescue requirements [59]. The HA losses in most areas were basically the same, ranging from 0 to 0.3. The largest loss was in Kathmandu and Lalitpur, with S greater than 0.2.
In RP-2, there was little difference in kr2 in most areas. Kathmandu and Lalitpur had the highest reconstruction rate, with the kr2 of 0.039 and 0.029, respectively. However, Nuwakot and Kavrepalanchok had the lowest reconstruction rate, with the kr2 of 0.005. For the relative reconstruction rate, the largest area was Kathmandu and the smallest area was Nuwakot.
In DP, Kathmandu and Lalitpur still maintain a high rate, with the kd of 0.039 and 0.029, respectively. Dhanusa and Makwanpur had the lowest development rate, with the kd of -0.02 and -0.021, respectively.
Throughout the recovery process, Kathmandu and Lalitpur showed high recovery capabilities compared to other areas. We know that Kathmandu and its surroundings are the core areas of Nepal. The PNLP indicators showed that the resilience of developed regions was significantly greater than that of less development regions

Influencing Factors of PNLP Pattern
From Figure 8, we can see that there were obvious differences in the recovery process. In order to explore the indicators caused by the differences in the recovery process, we collected Nepal's socioeconomic data from CBS and disaster data from MHA. ALI shows the level of development in a region. We represented the regional development according to the ALI mean in 2014. In general, we chose economic indicators, disaster data, population data, energy source for lighting, and type of house to analyze PNLP influencing factors (Figure 9). Indicator selection and meaning are shown in Table A1. K p was in a relatively stable range (from −0.01 to 0.01) in most areas. The fastest developing area before the earthquake was Bhaktapur (k p = 0.009), while Rasuwa, Parsa, and Dolakha were the areas with the slowest development (k p < −0.01).
In RP-1, recovery speed was slow in all regions, with k r1 ranging from 0 to 0.03. Kathmandu and Lalitpur had the highest recovery rate, with k r1 around 0.025, followed by regions Sarlahi and Lamjung, with k r1 greater than 0.01. In some remote areas k r1 was relatively small, which may have been due to the slow pace of Nepal's rescue of remote areas and the government services that cannot meet regional rescue requirements [59]. The HA losses in most areas were basically the same, ranging from 0 to 0.3. The largest loss was in Kathmandu and Lalitpur, with S greater than 0.2.
In RP-2, there was little difference in k r2 in most areas. Kathmandu and Lalitpur had the highest reconstruction rate, with the k r2 of 0.039 and 0.029, respectively. However, Nuwakot and Kavrepalanchok had the lowest reconstruction rate, with the k r2 of 0.005. For the relative reconstruction rate, the largest area was Kathmandu and the smallest area was Nuwakot.
In DP, Kathmandu and Lalitpur still maintain a high rate, with the k d of 0.039 and 0.029, respectively. Dhanusa and Makwanpur had the lowest development rate, with the k d of −0.02 and −0.021, respectively.
Throughout the recovery process, Kathmandu and Lalitpur showed high recovery capabilities compared to other areas. We know that Kathmandu and its surroundings are the core areas of Nepal. The PNLP indicators showed that the resilience of developed regions was significantly greater than that of less development regions.

Influencing Factors of PNLP Pattern
From Figure 8, we can see that there were obvious differences in the recovery process. In order to explore the indicators caused by the differences in the recovery process, we collected Nepal's socio-economic data from CBS and disaster data from MHA. ALI shows the level of development in a region. We represented the regional development according to the ALI mean in 2014. In general, we chose economic indicators, disaster data, population data, energy source for lighting, and type of house to analyze PNLP influencing factors (Figure 9). Indicator selection and meaning are shown in Table A1. From Figure 9, we can observe that the development before the earthquake was influenced by many factors. The factors affecting kp were mainly HPI (35.36%) and PG (10.17%). For kr1, the main influencing factor was CEM (13.81%), WOOD (12.76%), and ET (8.56%). The kr2 was also influenced by multiple factors, which were NHU (19.30%), BG (12.66%), ET (9.48%), and SD (8.33%). Compared with kp and kr1, disaster data (SD) became the main factor affecting kr2. In RP-1 and RP-2, energy source (ET) was always the main factor. Compared with the RP-1, although the influence of building type reduced in RP-2, the influence of the number of housing units (NHU) occupied the main position. In other words, architecture was the main factor affecting the recovery process [40]. The krp was influenced by multiple factors, including BG (15.38%), WOOD (9.95%), CEM (8.36%), PD (7.97%), and From Figure 9, we can observe that the development before the earthquake was influenced by many factors. The factors affecting k p were mainly HPI (35.36%) and PG (10.17%). For k r1 , the main influencing factor was CEM (13.81%), WOOD (12.76%), and ET (8.56%). The k r2 was also influenced by multiple factors, which were NHU (19.30%), BG (12.66%), ET (9.48%), and SD (8.33%). Compared with k p and k r1 , disaster data (SD) became the main factor affecting k r2 . In RP-1 and RP-2, energy source (ET) was always the main factor. Compared with the RP-1, although the influence of building type reduced in RP-2, the influence of the number of housing units (NHU) occupied the main position. In other words, architecture was the main factor affecting the recovery process [40]. The k rp was influenced by multiple factors, including BG (15.38%), WOOD (9.95%), CEM (8.36%), PD (7.97%), and TN (7.09%). The main factors affecting k d were KE (32.81%) and DAHA (26.12%). Energy source was always the main factor affecting the recovery process. However, compared with other phases, agriculture became the main factor affecting DP. This may have been because agriculture is the main economy mode of Nepal, so the level of agriculture affects the level of post-earthquake development. The main factors affecting S were HPI (20.68%), CEM (13.11%), ET (11.05%), and WOOD (10.45%).
In short, buildings (such as CEM and WOOD) and types of energy source (such as BG and ET) were the main factors affecting post-earthquake recovery. Agriculture had a certain influence in the recovery process, especially in DP. The gap between the rich and the poor was obviously a factor affecting the development rate before the earthquake, but had little effect on the post-earthquake process. Disaster data had an impact on RP-2.

Analysis of Noise Level
In this study, there are two different types of noise, one is background noise, and the other is the remainder and seasonal component of the STL algorithm. In order to analyze whether the two kinds of noise in the image affect each other, the correlation between the two noises is analyzed ( Figure 10).
Remote Sens. 2020, 12, x FOR PEER REVIEW 15 of 23 TN (7.09%). The main factors affecting kd were KE (32.81%) and DAHA (26.12%). Energy source was always the main factor affecting the recovery process. However, compared with other phases, agriculture became the main factor affecting DP. This may have been because agriculture is the main economy mode of Nepal, so the level of agriculture affects the level of post-earthquake development.
In short, buildings (such as CEM and WOOD) and types of energy source (such as BG and ET) were the main factors affecting post-earthquake recovery. Agriculture had a certain influence in the recovery process, especially in DP. The gap between the rich and the poor was obviously a factor affecting the development rate before the earthquake, but had little effect on the post-earthquake process. Disaster data had an impact on RP-2.

Analysis of Noise Level
In this study, there are two different types of noise, one is background noise, and the other is the remainder and seasonal component of the STL algorithm. In order to analyze whether the two kinds of noise in the image affect each other, the correlation between the two noises is analyzed ( Figure 10). From Figure 10, we can observe that there was a positive correlation (p <0.05, R = 0.32) between background noise and STL remainder and seasonal component. There is a certain relationship between the selection method of background noise and the STL algorithm. In this study, the application of the two methods is more reasonable.

Mutation Test
STL algorithm divides the original time series data into seasonal component, trend component, and remainder component based on a locally weighted regression smoother (LOESS). This method assumes that the trend component is smooth and slowly changing. Compared with the STL algorithm, breaks for additive seasonal and trend (BFAST) can monitor the trend break in the trend component, that is, multiple mutation detection of the trend [60]. From Figure 10, we can observe that there was a positive correlation (p <0.05, R = 0.32) between background noise and STL remainder and seasonal component. There is a certain relationship between the selection method of background noise and the STL algorithm. In this study, the application of the two methods is more reasonable.

Mutation Test
STL algorithm divides the original time series data into seasonal component, trend component, and remainder component based on a locally weighted regression smoother (LOESS). This method assumes that the trend component is smooth and slowly changing. Compared with the STL algorithm, breaks for additive seasonal and trend (BFAST) can monitor the trend break in the trend component, that is, multiple mutation detection of the trend [60].
An earthquake is a sudden disaster that may cause the structure of the ALI sequence. In this study, we utilized the STL algorithm instead of the BFAST algorithm for analysis and research. In order to determine whether the STL algorithm was applicable in this study, we utilized the MK mutation test to check whether there was a mutation point in the original ALI sequence after the earthquake (Figure 11). An earthquake is a sudden disaster that may cause the structure of the ALI sequence. In this study, we utilized the STL algorithm instead of the BFAST algorithm for analysis and research. In order to determine whether the STL algorithm was applicable in this study, we utilized the MK mutation test to check whether there was a mutation point in the original ALI sequence after the earthquake (Figure 11). Figure 11. Mk mutation test. The red curve represents the UF (UFk statistical series) and the blue curve represents the UB (UBk statistical series). The dotted lines represent 95% confidence intervals. In the confidence interval, the intersection of UF and UB is the breakpoint.
We performed the MK mutation test on the original ALI sequence in each region. From Figure  11, we can see that the UF and UB curves had no intersection near the earthquake month (April 2015). That is, the earthquake did not cause the mutation of the original ALI sequence. This may have been due to the spatial-time scale applied in this study. Li et al. [35] adopted NTL data to study the annual change of the pixel scale. They concluded that the restoration led to an increase in the intensity of HA after the earthquake, but the increase level was not high. In this study, the spatial scale was a regional scale that was coarser than the pixel scale. Although we utilized monthly data, the emergency time was only a few weeks.
Although there was no breakpoint around April 2015 in Figure 11, we can clearly observe that there were breakpoints after January 2017. It may have been due to the step factor. Through the step factor, the ALI values were reduced from January 2017 to December 2018, but the ALI trend was not be affected. We found that in most regions the breakpoints appeared around June 2017 and were Figure 11. Mk mutation test. The red curve represents the UF (UF k statistical series) and the blue curve represents the UB (UB k statistical series). The dotted lines represent 95% confidence intervals. In the confidence interval, the intersection of UF and UB is the breakpoint.
We performed the MK mutation test on the original ALI sequence in each region. From Figure 11, we can see that the UF and UB curves had no intersection near the earthquake month (April 2015). That is, the earthquake did not cause the mutation of the original ALI sequence. This may have been due to the spatial-time scale applied in this study. Li et al. [35] adopted NTL data to study the annual change of the pixel scale. They concluded that the restoration led to an increase in the intensity of HA after the earthquake, but the increase level was not high. In this study, the spatial scale was a regional scale that was coarser than the pixel scale. Although we utilized monthly data, the emergency time was only a few weeks.
Although there was no breakpoint around April 2015 in Figure 11, we can clearly observe that there were breakpoints after January 2017. It may have been due to the step factor. Through the step factor, the ALI values were reduced from January 2017 to December 2018, but the ALI trend was not be affected. We found that in most regions the breakpoints appeared around June 2017 and were close to 2018. If the breakpoint was due to the step factor, it should have been near January 2017. So, we do not think it is the data breakpoint caused by the step factor. As we can see from Figure 7, near January 2018 it was the dividing point between RP-2 and DP. Here, the ALI trend changed. We think that the data breakpoint was caused by the phase change. In the PNLP pattern, the differences between the phases were not considered. We ignored the changes between phases caused by this breakpoint. So, we can infer that the STL algorithm can be better applied to this study.

Seasonal Impact Analysis
NTL time series data is seasonally affected by snow cover, albedo and NDVI. The NDVI data was selected to check whether the trend is affected by seasonal changes (Figure 12).
Remote Sens. 2020, 12, x FOR PEER REVIEW 17 of 23 close to 2018. If the breakpoint was due to the step factor, it should have been near January 2017. So, we do not think it is the data breakpoint caused by the step factor. As we can see from Figure 7, near January 2018 it was the dividing point between RP-2 and DP. Here, the ALI trend changed. We think that the data breakpoint was caused by the phase change. In the PNLP pattern, the differences between the phases were not considered. We ignored the changes between phases caused by this breakpoint. So, we can infer that the STL algorithm can be better applied to this study.

Seasonal Impact Analysis
NTL time series data is seasonally affected by snow cover, albedo and NDVI. The NDVI data was selected to check whether the trend is affected by seasonal changes (Figure 12). In the research of Levin [36], he showed that there was a weak negative correlation between VIIRS monthly data and mean NDVI values. In our study, from Figure 12a, we observed that there was no correlation between mean NDVI values and trend component. However, from Figure 12d, we found that there was a weak negative correlation between mean NDVI values and sum of seasonal and random component (P <0.01, R = −143). This result is basically consistent with the result of Levin. From this perspective, we believe that through the STL algorithm, the seasonal effects in the VIIRS data were decomposed into the seasonal component and random component of STL, and the trend component was less affected by the seasons. However, it must be mentioned that Levin obtained a higher correlation between VIIRS data and snow cover, albedo, and NDVI through principal component analysis. This was not our focus, so we only utilized NDVI data to observe seasonality.

Trend Analysis
As we mentioned above, the step change was removed. In order to observe if the trend was caused by an earthquake, we selected three regions and analyzed their trends. Figure 13 is shown below. From Figure 13a-c, HA intensity showed an upward trend. Their manifestations were In the research of Levin [36], he showed that there was a weak negative correlation between VIIRS monthly data and mean NDVI values. In our study, from Figure 12a, we observed that there was no correlation between mean NDVI values and trend component. However, from Figure 12d, we found that there was a weak negative correlation between mean NDVI values and sum of seasonal and random component (P <0.01, R = −143). This result is basically consistent with the result of Levin. From this perspective, we believe that through the STL algorithm, the seasonal effects in the VIIRS data were decomposed into the seasonal component and random component of STL, and the trend component was less affected by the seasons. However, it must be mentioned that Levin obtained a higher correlation between VIIRS data and snow cover, albedo, and NDVI through principal component analysis. This was not our focus, so we only utilized NDVI data to observe seasonality.

Trend Analysis
As we mentioned above, the step change was removed. In order to observe if the trend was caused by an earthquake, we selected three regions and analyzed their trends. Figure 13 is shown below. From Figure 13a-c, HA intensity showed an upward trend. Their manifestations were different from those in Nepal. Moreover, we learned that if no major disaster occurred, HA showed a steady upward trend. In Nepal, the trend changes were monitored around April 2015. That is, we believe that the earthquake caused the trend change in Nepal.
Remote Sens. 2020, 12, x FOR PEER REVIEW 18 of 23 different from those in Nepal. Moreover, we learned that if no major disaster occurred, HA showed a steady upward trend. In Nepal, the trend changes were monitored around April 2015. That is, we believe that the earthquake caused the trend change in Nepal. However, a problem must be mentioned here. We did not use ground truth data to verify the accuracy of the trend. However, we know that the changes of NTL are closely related to HA, and here we directly use the trend changes to explain the HA changes. In the future, what factors NTL can reflect should be considered.
It must be said that not only seasonal changes and the step affect our trend. There are other factors that affect the trend, such as image angle, errors in calibration of the imaging sensor, and "blooming". For example, Xiong et al. [61] performed on-orbit calibration on VIIRS DNB band. In our research, we did not consider the influence of these factors on the trend. Through the STL algorithm, we believe that the factors that affect the HA trend were eliminated. In future studies, understanding which factors affect the VIIRS data will help to remove noise.

Diversity of PNLP Conceptual Pattern
Linear and orderly PNLP pattern reflects the change trend of HA intensity after the earthquake. It can compare the recovery process in different areas after the earthquake. In this study, we developed the NPP/VIIRS monthly data to construct the PNLP pattern. However, the time taken in the EP is relatively short and its change characteristics cannot be described in the PNLP pattern. Therefore, the application of daily data can make the pattern more refined and complete. VIIRS daily Black Mable product has great potential in constructing the PNLP pattern [24].
The construction of PNLP pattern is based on Haas's disaster framework. However, with regard to Haas's framework, subsequent case studies questioned the idea of orderly and inevitable progress in the recovery [62,63]. In other words, if there are new external shocks during the recovery process, PNLP pattern will change, or the entire PNLP may lack a certain phase due to disaster intensity and social development.
PNLP is a linear pattern, but disaster recovery is often a non-linear and multidimensional process [64]. Therefore, each phase of the pattern may no longer be linear, but a non-linear process.
PNLP is a linear and orderly pattern based on the NTL data. It can provide an effective means for monitoring the recovery differences between different regions. In order to better monitor the progress of earthquake recovery, the PNLP pattern can be developed into many different forms based on more in-depth study of NTL data and disaster recovery theory. Therefore, considering multiple However, a problem must be mentioned here. We did not use ground truth data to verify the accuracy of the trend. However, we know that the changes of NTL are closely related to HA, and here we directly use the trend changes to explain the HA changes. In the future, what factors NTL can reflect should be considered.
It must be said that not only seasonal changes and the step affect our trend. There are other factors that affect the trend, such as image angle, errors in calibration of the imaging sensor, and "blooming". For example, Xiong et al. [61] performed on-orbit calibration on VIIRS DNB band. In our research, we did not consider the influence of these factors on the trend. Through the STL algorithm, we believe that the factors that affect the HA trend were eliminated. In future studies, understanding which factors affect the VIIRS data will help to remove noise.

Diversity of PNLP Conceptual Pattern
Linear and orderly PNLP pattern reflects the change trend of HA intensity after the earthquake. It can compare the recovery process in different areas after the earthquake. In this study, we developed the NPP/VIIRS monthly data to construct the PNLP pattern. However, the time taken in the EP is relatively short and its change characteristics cannot be described in the PNLP pattern. Therefore, the application of daily data can make the pattern more refined and complete. VIIRS daily Black Mable product has great potential in constructing the PNLP pattern [24].
The construction of PNLP pattern is based on Haas's disaster framework. However, with regard to Haas's framework, subsequent case studies questioned the idea of orderly and inevitable progress in the recovery [62,63]. In other words, if there are new external shocks during the recovery process, PNLP pattern will change, or the entire PNLP may lack a certain phase due to disaster intensity and social development.
PNLP is a linear pattern, but disaster recovery is often a non-linear and multidimensional process [64]. Therefore, each phase of the pattern may no longer be linear, but a non-linear process.
PNLP is a linear and orderly pattern based on the NTL data. It can provide an effective means for monitoring the recovery differences between different regions. In order to better monitor the progress of earthquake recovery, the PNLP pattern can be developed into many different forms based on more in-depth study of NTL data and disaster recovery theory. Therefore, considering multiple factors, such as time scale, stage diversity, and discontinuity, a more refined and comprehensive model should be established in the future.

Conclusions
It is important to study the characteristics of phase changes during the process of earthquake recovery and reconstruction for post-earthquake work deployment arrangements. NTL has a strong correlation with HA. Based on the NPP/VIIRS monthly composite data and the disaster recovery theory, this paper constructs a PNLP pattern.
This paper evaluates the application of the PNLP pattern to the Nepal earthquake. The results show that the PNLP model can reflect the changes in HA intensity after the earthquake in Nepal. The PNLP pattern divides the intensity of HA after the Nepal earthquake into EP, RP-1, RP-2, and DP. Based on the PNLP indicators, we can conclude that the recovery process in developed regions in significantly higher than in less developed regions. There are certain differences in the factors that affect different PNLP indicators. But in general, it is mainly related to the type and number of the building, the type of energy source, economic (especially the agricultural economy) and disaster damage.
Based on the NPP/VIIRS data, this paper proposes a PNLP pattern to analyze the spatiotemporal characteristics of short-term dynamic sequences after an earthquake. In order to improve the quality of NTL data, this paper applies the STL algorithm to remove the seasonal and noise effects of ALI continuous data from January 2014 to December 2018. Nevertheless, research deficiencies must be recognized. Due to the coarse space-time scale, the earthquake recovery mutation cannot be detected. The PNLP pattern may produce diversity changes due to regional economic differences. Therefore, the NPP/VIIRS daily data and Lujia-1 data have greater potential for application research in this context because of their spatiotemporal characteristics.

Conflicts of Interest:
The authors declare no conflicts of interest.