Identifying and Classifying Shrinking Cities Using Long-Term Continuous Night-Time Light Time Series

: Shrinking cities—cities suffering from population and economic decline—has become a pressing societal issue of worldwide concern. While night-time light (NTL) data have been applied as an important tool for the identiﬁcation of shrinking cities, the current methods are constrained and biased by the lack of using long-term continuous NTL time series and the use of unidimensional indices. In this study, we proposed a novel method to identify and classify shrinking cities by long-term continuous NTL time series and population data, and applied the method in northeastern China (NEC) from 1996 to 2020. First, we established a long-term consistent NTL time series by applying a geographically weighted regression model to two distinct NTL datasets. Then, we generated NTL index (NI) and population index (PI) by random forest model and the slope of population data, respectively. Finally, we developed a shrinking city index (SCI), based on NI and PI to identify and classify city shrinkage. The results showed that the shrinkage pattern of NEC in 1996–2009 (stage 1) and 2010–2020 (stage 2) was quite different. From stage 1 to stage 2, the shrinkage situation worsened as the number of shrinking cities increased from 102 to 162, and the proportion of severe shrinkage increased from 9.2% to 30.3%. In stage 2, 85.4% of the cities exhibited population decline, and 15.7% of the cities displayed an NTL decrease, suggesting that the changes in NTL and population were not synchronized. Our proposed method provides a robust and long-term characterization of city shrinkage and is beneﬁcial to provide valuable information for sustainable urban planning and decision-making.


Introduction
Shrinking city, a term proposed in the early 1990s [1], refers to cities suffering from population and economic decline. City shrinkage has been identified in more than a quarter of the world's cities, and this percentage is still rising [2], particularly in developed countries [3][4][5][6]. Wolff et al. [7] confirmed that 54% of the urban areas in Europe were experiencing population loss. Pallagst et al. [8] noted that 13% of the urban regions in the US shrank, among which many historical industrial cities, such as Detroit, were in the most serious condition [9]. Recently, shrinking cities have also been found in developing countries [10][11][12]. Among them, cities in northeastern China where old industrial bases are located were typical cases [13,14]. City shrinkage has a profound impact on economic development, social differentiation, urban renewal, government decision-making, infrastructure construction and urban planning development; thus, this phenomenon has received increasing attention worldwide in recent decades [15][16][17][18][19].
describe the characteristics of shrinking cities, and variations among different factors cannot be identified, thus limiting detailed analysis. In addition, the construction of an NTL index in most previous studies was merely based on the average digital number (DN) value [27,36]. There are many NTL features at the regional scale, but little attention has been given to selecting representative features to generate a comprehensive indicator. This limitation may cause some information related to classifying shrinking cities with NTL data to not be thoroughly extracted.
To fill these research gaps, this study proposed a novel method to identify and classify shrinking cities by combining long-term time series of NTL data from 1996 to 2020 with auxiliary data. To achieve these goals, we (1) cross-calibrated different NTL sensors and generated long-term consistent NTL time series; (2) developed robust standards to identify and classify shrinking cities based on multiple variables; and (3) analyzed the spatiotemporal patterns of shrinking cities. This study aims to reveal the long-term characteristics of city shrinkage by using continuous sets of NTL data and auxiliary data and provides valuable insight for government decision makers regarding shrinking city development.

Study Area
Heilongjiang, Jilin and Liaoning, three provinces in northeastern China (NEC), were selected as the study area ( Figure 1). Benefitted by vast plains and rich resources, this region started its urbanization at an early age. NEC includes 34 prefecture-level cities today and covers an area of 806,100 km 2 . Considering the access to statistical data, we chose county (district) level as the research scale. We merged all the districts belonging to a specific prefecture level city and used the prefecture city's name to name its corresponding merged districts. Thus, the basic unit of our study consists of three parts: county, county-level city and merged districts. We made some adjustments to assure that statistical data in the whole period matched the administrative boundary. To be specific, the population and GDP of withdrawn counties (districts) were allocated into the surrounding cities by using the county (district)-level administrative boundary of the NEC in 2018 as standard. We finally obtained 185 basic study unites. Among them, 9 cities are big cities (urban resident population is over 1 million), 17 cities are medium cities (urban resident population is between 0.5 million and 1 million), and others are small cities (urban resident population is less than 0.5 million).
NEC was known as the old industrial basement of China and used to be an important pillar of China's economy. Many of the cities in NEC are important resource-based industrial cities such as Fushun, Anshan, Jixi and Hegang, and the economy of these cities are severely dependent on heavy industry. Exemplified by Fushun, the value added of the second industry accounts for 43% of its GDP in 2018. In recent decades, due to resource exhaustion and global economic downturns, many of the cities in these provinces have suffered from population decline and economic recession. According to the China National Statistical Yearbooks, the total population in NEC decreased from 109.7 million in 2012 to 105.8 million in 2018 [39,40], while the national population increased rapidly during the same period. In terms of the economy, the total GDP of NEC in 2018 was USD 8812.1 billion with a growth rate of 5.5%; both were lower than the national average level. Therefore, NEC are ideal sites to study city shrinkage and test our method.

Datasets
The datasets used in the study are summarized in Table 1, which include NTL data, social census data and human settlement data. NTL data include the global radiance-calibrated DMSP-OLS (DMSPgrc) and NPP/VIIRS products. DMSPgrc are yearly composites with a spatial resolution of 30 arc-seconds. Compared to the frequently used stable light product of DMSP-OLS, DMSPgrc products have better data quality. By combining sparse data obtained at low-gain settings with the operational data obtained at high-gain settings, the dynamic range of DMSPgrc was improved, and the saturation effects were eliminated [38]. DMSPgrc were only available for eight years, 1996, 1999, 2000, 2003, 2004, 2006, 2010 and 2011, and all the data in these years were included in this study. Version 2 (V2) yearly composite VIIRS data were collected from 2012 to 2020. VIIRS provide NTL measurements in radiance and have a higher spatial resolution of 15 arc-seconds than DMSP. The effects of lunar illumination, stray light, lightning and cloud cover were eliminated [33]. Compared to Version 1 VIIRS monthly data, the V2 annually VIIRS data applied outlier removal to exclude background noise [41]. Social census data, including population data and GDP data at the county scale, were collected from the China Statistical Yearbook. To further analyze the land expansion condition in shrinking cities, we used human settlement data produced by Gong et al. [42]. The data have a resolution of 30 m and an overall accuracy of 90%.

Datasets
The datasets used in the study are summarized in Table 1, which include NTL data, social census data and human settlement data. NTL data include the global radiancecalibrated DMSP-OLS (DMSP grc ) and NPP/VIIRS products. DMSP grc are yearly composites with a spatial resolution of 30 arc-seconds. Compared to the frequently used stable light product of DMSP-OLS, DMSP grc products have better data quality. By combining sparse data obtained at low-gain settings with the operational data obtained at high-gain settings, the dynamic range of DMSP grc was improved, and the saturation effects were eliminated [38]. DMSP grc were only available for eight years, 1996, 1999, 2000, 2003, 2004, 2006, 2010 and 2011, and all the data in these years were included in this study. Version 2 (V2) yearly composite VIIRS data were collected from 2012 to 2020. VIIRS provide NTL measurements in radiance and have a higher spatial resolution of 15 arc-seconds than DMSP. The effects of lunar illumination, stray light, lightning and cloud cover were eliminated [33]. Compared to Version 1 VIIRS monthly data, the V2 annually VIIRS data applied outlier removal to exclude background noise [41]. Social census data, including population data and GDP data at the county scale, were collected from the China Statistical Yearbook. To further analyze the land expansion condition in shrinking cities, we used human settlement data produced by Gong et al. [42]. The data have a resolution of 30 m and an overall accuracy of 90%.

Methodology
In this paper, "city shrinkage" was defined as a phenomenon in which a city experiences either an NTL decline or a population decrease in a certain stable period. The rationale of using population data from censuses instead of population grids based on remote sensing (e.g., Landscan) was that population grids based on remote sensing are more unstable at the regional scale than population data based on censuses. The method comprised three parts: (1) building consistent long-term time series (from 1996 to 2020) of NTL data based on DMSP grc data and VIIRS data, (2) identifying shrinking cities and (3) classifying shrinking cities. Figure 2 shows the flowchart of this study.

Methodology
In this paper, "city shrinkage" was defined as a phenomenon in which a city experiences either an NTL decline or a population decrease in a certain stable period. The rationale of using population data from censuses instead of population grids based on remote sensing (e.g., Landscan) was that population grids based on remote sensing are more unstable at the regional scale than population data based on censuses. The method comprised three parts: (1) building consistent long-term time series (from 1996 to 2020) of NTL data based on DMSPgrc data and VIIRS data, (2) identifying shrinking cities and (3) classifying shrinking cities. Figure 2 shows the flowchart of this study.

Figure 2.
Flowchart of this study (NI refers to the NTL index, PI refers to the population index and SCI refers to the shrinking city index).

Data Pre-Processing
There were two tasks in the pre-processing stage: (1) DMSPgrc and VIIRS data crosssensor calibration and (2) trend assessment for both datasets followed by stage division, if necessary. Geographically weighted regression (GWR) was chosen as the cross-sensor calibration method; GWR is a commonly used regression model that considers spatial non-stationarity [43] and has been proven to be effective in the calibration of NTL data [38]. A schematic diagram of the whole process is shown in Figure 3.
In the pre-processing of DMSPgrc data, first, saturation effects and geometric errors were excluded, and the invariant region method was then applied for inter-calibration. DMSPgrc only include data for eight years, and thus, a logistical model was used to estimate data for the remaining years from 1996 to 2013. The expression is as follows:

Data Pre-Processing
There were two tasks in the pre-processing stage: (1) DMSP grc and VIIRS data crosssensor calibration and (2) trend assessment for both datasets followed by stage division, if necessary. Geographically weighted regression (GWR) was chosen as the cross-sensor calibration method; GWR is a commonly used regression model that considers spatial non-stationarity [43] and has been proven to be effective in the calibration of NTL data [38]. A schematic diagram of the whole process is shown in Figure 3.
In the pre-processing of DMSP grc data, first, saturation effects and geometric errors were excluded, and the invariant region method was then applied for inter-calibration. DMSP grc only include data for eight years, and thus, a logistical model was used to estimate data for the remaining years from 1996 to 2013. The expression is as follows: where DMSP grc is the DN value of a specific pixel, e is the base of the natural logarithm, t is the year, and a, b, c, d are logistical regression coefficients, and they were determined by iterating the parameters continuously until the optimal group was reached; the fitted Remote Sens. 2021, 13, 3142 6 of 19 data were labeled as DMSP grc '. Considering that the fitting effect of some pixels may have estimation errors, pixels within fitting curves that have a coefficient of determination (R 2 ) lower than 0.85 were replaced by the average of the surrounding eight pixels according to Zheng et al. [38]. Although most of the background noises of the V2 VIIRS data have been removed [41], we made some pre-processing to further improve the data. Since Shenyang, Changchun and Harbin are the three biggest and most developed cities in NEC, the pixel values of the other areas in NEC should not exceed these cities. The highest DN value of those three cities in every year was used as a threshold to correct the outliers. In addition, values less than zero or other extreme values were removed. Then, the GWR model was applied to DMSP grc ' (resampled to 500 m) and VIIRS in the overlapping years 2012 and 2013; the GWR expression is as follows: where u i and v i denote the spatial location of the ith pixel, β 0 and β 1 stand for the intercept and slope, respectively, DN V I IRS and DN DMSP grc indicate the DN value of VIIRS and DMSP grc ' and ε i is the residual of the model. By applying the regression coefficients for all other years, NTL data for a complete time series from 1996 to 2020 with a resolution of 500 m were finally obtained. where DMSPgrc is the DN value of a specific pixel, e is the base of the natural logarithm, t is the year, and a, b, c, d are logistical regression coefficients, and they were determined by iterating the parameters continuously until the optimal group was reached; the fitted data were labeled as DMSPgrc '. Considering that the fitting effect of some pixels may have estimation errors, pixels within fitting curves that have a coefficient of determination (R 2 ) lower than 0.85 were replaced by the average of the surrounding eight pixels according to Zheng et al. [38]. Although most of the background noises of the V2 VIIRS data have been removed [41], we made some pre-processing to further improve the data. Since Shenyang, Changchun and Harbin are the three biggest and most developed cities in NEC, the pixel values of the other areas in NEC should not exceed these cities. The highest DN value of those three cities in every year was used as a threshold to correct the outliers. In addition, values less than zero or other extreme values were removed. Then, the GWR model was applied to DMSPgrc' (resampled to 500 m) and VIIRS in the overlapping years 2012 and 2013; the GWR expression is as follows: where and denote the spatial location of the ith pixel, 0 and 1 stand for the intercept and slope, respectively, and ′ indicate the DN value of VIIRS and DMSPgrc' and is the residual of the model. By applying the regression coefficients for all other years, NTL data for a complete time series from 1996 to 2020 with a resolution of 500 m were finally obtained. The trends of both datasets were checked, and we found that although the light curves did not display an obvious trend, most of the city population trajectories exhibited a noticeable two-stage pattern in which the population increased initially and declined afterwards ( Figure 4). It is worth noting that the population peaks in 125 cities, accounting for 67.8% of all the cities studied, appeared between 2009 and 2011. Considering that a serious world financial crisis occurred in 2008, it is reasonable to assume that the city shrinkage in NEC was affected by this crisis [13]. Previous studies also noted that cities in NEC entered a transition period in 2010 [13]. Therefore, we divided the time series into two stages: S1 (1996-2009) and S2 (2010-2020). The trends of both datasets were checked, and we found that although the light curves did not display an obvious trend, most of the city population trajectories exhibited a noticeable two-stage pattern in which the population increased initially and declined afterwards ( Figure 4). It is worth noting that the population peaks in 125 cities, accounting for 67.8% of all the cities studied, appeared between 2009 and 2011. Considering that a serious world financial crisis occurred in 2008, it is reasonable to assume that the city shrinkage in NEC was affected by this crisis [13]. Previous studies also noted that cities in NEC entered a transition period in 2010 [13]. Therefore, we divided the time series into two stages: S1 (1996-2009) and S2 (2010-2020).

Construction of the NTL Index (NI) and Population Index (PI)
Reasonable quantitative indices must be established to identify shrinking cities via long-term time series of population data and NTL data. We generated the following indices in this study.
(1) NTL index (NI) NTL data have many feature values, but most previous studies only chose the average DN as the indicator for identifying shrinking cities. In this study, a random forest (RF) method was applied to search the NTL features with the highest discrimination abilities for shrinking cities. The RF approach is a classic machine learning method that uses multiple trees for training and sample prediction and has been widely used in classification [44]. We calculated eight candidate NTL features, including the slope of the total DN, the intercept of the total DN, the slope of the average DN, the intercept of the average DN, the slope of the top 10% of DN values, the intercept of the top 10% of DN values, the slope of light pixel areas and the intercept of light pixel areas. Since there is no recognized shrinking city dataset, ten typical shrinking cities and thirty non-shrinking cities identified by Long et al. [12] were chosen as training samples. Notably, Long et al. [12] used census data to identify shrinking cities in China, and the results of their study have been supported by many related studies. The population census data and the GDP data for the ten typical shrinking cities were further checked to ensure their representativeness as shrinking cities. The R 2 was used to judge the accuracy of the model, and the out-of-bag (OOB) error was used to assess the relative importance of all features. According to experience, NTL features with a normalized OOB higher than 0.5 were chosen as important features. Then, the NI was constructed by using the selected features to calculate modular distance: where denotes the number of selected features, 1 2 , …, and are the top-n features of highest importance (all normalized) and 1 , 2 , and 3 are the weights. Here, we used the reciprocals of the OOB errors as the weights for different features. For a specific city, if either of the value of the selected features is negative, we define the NI of this city is negative, otherwise it is positive.
(2) Population index (PI) We defined the slope changes for the registered residence population in each city during a period by the PI (population index):

Construction of the NTL Index (NI) and Population Index (PI)
Reasonable quantitative indices must be established to identify shrinking cities via long-term time series of population data and NTL data. We generated the following indices in this study.
(1) NTL index (NI) NTL data have many feature values, but most previous studies only chose the average DN as the indicator for identifying shrinking cities. In this study, a random forest (RF) method was applied to search the NTL features with the highest discrimination abilities for shrinking cities. The RF approach is a classic machine learning method that uses multiple trees for training and sample prediction and has been widely used in classification [44]. We calculated eight candidate NTL features, including the slope of the total DN, the intercept of the total DN, the slope of the average DN, the intercept of the average DN, the slope of the top 10% of DN values, the intercept of the top 10% of DN values, the slope of light pixel areas and the intercept of light pixel areas. Since there is no recognized shrinking city dataset, ten typical shrinking cities and thirty non-shrinking cities identified by Long et al. [12] were chosen as training samples. Notably, Long et al. [12] used census data to identify shrinking cities in China, and the results of their study have been supported by many related studies. The population census data and the GDP data for the ten typical shrinking cities were further checked to ensure their representativeness as shrinking cities. The R 2 was used to judge the accuracy of the model, and the out-of-bag (OOB) error was used to assess the relative importance of all features. According to experience, NTL features with a normalized OOB higher than 0.5 were chosen as important features. Then, the NI was constructed by using the selected features to calculate modular distance: where n denotes the number of selected features, Fea 1 Fea 2 , . . . , and Fea n are the top-n features of highest importance (all normalized) and a 1 , a 2 , and a 3 are the weights. Here, we used the reciprocals of the OOB errors as the weights for different features. For a specific city, if either of the value of the selected features is negative, we define the NI of this city is negative, otherwise it is positive.
(2) Population index (PI) We defined the slope changes for the registered residence population in each city during a period by the PI (population index): where x i denotes the ith year, x denotes the average over all years, y i denotes the population in the ith year and y denotes the average population in the whole period.

Shrinking City Classification
Based on the characteristics of the NI and PI in a given period, all the cities were divided into four categories marked as non-shrinkage (the value of NI is positive, and the value of PI is positive), P-shrinkage (the value of NI is positive, and the value of PI is negative), N-shrinkage (the value of NI is negative, and the value of PI is positive) and NP-shrinkage (the value of NI is negative, and the value of PI is negative). To further evaluate the severity of city shrinkage, the shrinking city index (SCI) was established to combine the PI and NI and reflect the degree of shrinkage according to the modular distance between the two indices: where NI is the NTL index and expresses NTL shrinkage and PI is the population index and reflects population shrinkage (all normalized). After the SCI values of all the cities were calculated, the natural breaks method was used to categorize the cities into three levels: slight shrinkage, medium shrinkage and severe shrinkage. Whether an indicator can appropriately reflect the characteristics of shrinking cities is mainly determined by its sensitivity to population decreases and declines in economic activities. Considering that the SCI was generated from the PI, the SCI itself can reflect population declines. Thus, we need to verify the consistency of this index with changes of economic activities to further test its robustness. GDP data were applied to express the economic activities of the cities, and the slope of changes in each city's GDP during a period were calculated as the GI (GDP index): where x i denotes the ith year, x denotes the average over all years, y i denotes the GDP in the ith year and y denotes the average GDP in the whole period. Then, linear regression was used to verify the consistency of the two indices.

Performance Assessment of the Proposed Method
The effect of the cross-sensor calibration was checked by two aspects: the curve of the time series and the quality of the calibrated image. Figure 5, exemplified by Changchun city (the capital city of Jilin Province), illustrates the NTL time series after calibration. The NTL between the different satellite sensors is well calibrated (Figure 5a). Additionally, the calibrated DMSP grc in Figure 5d has a better data quality than its uncalibrated counterpart (Figure 5b) in terms of a higher spatial resolution and more detailed spatial information.
For the RS results, the R 2 of the model was 0.97, which met the accuracy requirement. The OOB errors of all the candidate NTL features are shown in Figure 6. Three NTL features have a normalized OOB error higher than 0.5, including the slope of the total DN, the intercept of the total DN and the slope of light pixels. Thus, they were chosen to build the NI.
The R 2 value of the linear regression model is 0.71 (Figure 7), which demonstrates that there was a certain consistency between the GI and SCI. The main deviation points are associated with cities with high GI values; this can be explained by the fact that areas where NTL was not detected may still generate economic activity during the day. Remote Sens. 2021, 13, x FOR PEER REVIEW 9 of 20 For the RS results, the R 2 of the model was 0.97, which met the accuracy requirement. The OOB errors of all the candidate NTL features are shown in Figure 6. Three NTL features have a normalized OOB error higher than 0.5, including the slope of the total DN, the intercept of the total DN and the slope of light pixels. Thus, they were chosen to build the NI. The R 2 value of the linear regression model is 0.71 (Figure 7), which demonstrates that there was a certain consistency between the GI and SCI. The main deviation points are associated with cities with high GI values; this can be explained by the fact that areas where NTL was not detected may still generate economic activity during the day.  For the RS results, the R 2 of the model was 0.97, which met the accuracy requirement. The OOB errors of all the candidate NTL features are shown in Figure 6. Three NTL features have a normalized OOB error higher than 0.5, including the slope of the total DN, the intercept of the total DN and the slope of light pixels. Thus, they were chosen to build the NI. The R 2 value of the linear regression model is 0.71 (Figure 7), which demonstrates that there was a certain consistency between the GI and SCI. The main deviation points are associated with cities with high GI values; this can be explained by the fact that areas where NTL was not detected may still generate economic activity during the day.

Spatial Patterns of Shrinking Cities
The numbers and distributions of the four types of cities were successfully identified

Spatial Patterns of Shrinking Cities
The numbers and distributions of the four types of cities were successfully identified by the proposed framework (Figure 8a,d). Overall, most of the cities have been shrinking continuously over the past 25 years. In S1, 102 shrinking cities were identified, accounting for 55.1% of the total number of cities studied. This finding indicates that city shrinkage in NEC has occurred for a long time and not just in recent years. Most of the shrinking cities were located in the northeastern and central part of NEC. The three types of shrinkage, namely, NP-shrinkage, P-shrinkage and N-shrinkage, were evenly distributed, and the numbers of the three types were 26, 37 and 39, respectively. A total of 34.1% of the cities in NEC displayed a population decline, and 35.1% of them were found to shrink according to NTL. In S2, there was a sharp increase in the number of shrinking cities, and the shrinkage pattern displayed considerable heterogeneity. Among the 185 cities, 162 cities were found to be shrinking, and only 12.4% were non-shrinking cities. A total of 133 cities were Pshrinking cities, accounting for 72% of all cities in this scenario, and only 4 cities were N-shrinking cities. The proportions of cities with decreases in population and NTL were 85.4% and 15.7%, respectively. The study results highlight a significant difference in the shrinking patterns between the two stages, suggesting that the stage division was rational to some extent. ing cities also displayed a significant transformation; 71 cities changed to other types, accounting for 85.5% of the original number. The main shift was toward P-shrinkage, with 102 cities displaying this pattern, much more than any other change type observed.  As for the shrinkage situation at the provincial level (Figure 8e,f), in S1, all three provinces had shrinking city proportions higher than 50%. Liaoning had 31 shrinking cities out of 58 total, or 53.4% of all corresponding cities, which was the lowest percentage among the three provinces. There were 7, 11 and 8 NP-shrinking cities in the three provinces, respectively. The NP-shrinking cities in Jilin accounted for 37.9% of all the shrinking cities, and this value was the highest among those for the three provinces. Eighteen N-shrinking cities and seventeen P-shrinking cities were found in Heilongjiang, more than in other provinces. In S2, NP-shrinking cities were nearly equally distributed in the three provinces. The numbers of P-shrinkage were 58, 38 and 37 in Heilongjiang, Jilin, and Liaoning, respectively, all of which accounted for large proportions of cities in these provinces; among them, Heilongjiang displayed the highest percentage. There were few N-shrinking cities in the three provinces, and there was a sharp increase in the P-shrinkage type. These results indicate that the leading cause of shrinkage in most cities in S2 was a population decrease.
Considering that the patterns of shrinkage in the two stages were quite different, we further studied the transformation situation of the four types of shrinkage from stage 1 to stage 2 ( Figure 9). The largest number of type changes was found among non-shrinking cities shifting to P-shrinkage, and the number was 58, accounting for 70% of all the nonshrinkage cities in S1. The N-shrinkage category had the largest number of cities that shifted to other shrinkage types, with all the 39 cities changing to other types. Non-shrinking cities also displayed a significant transformation; 71 cities changed to other types, accounting for 85.5% of the original number. The main shift was toward P-shrinkage, with 102 cities displaying this pattern, much more than any other change type observed.

Degree of Shrinking Cities
The degree of city shrinkage in the NEC in the two periods is presented in Figure  10a,b. It is obvious that the extent of shrinkage became more severe from S1 to S2. During S1, most (77.5%) of the shrinking cities in S2 displayed medium shrinkage, and slightly and severe shrinking cities accounted for 6% and 16.7% of all cities, respectively. The severe shrinking cities were mainly located in Jilin Province. Figure 10c displays the shrinkage degree at the provincial level: the proportion of medium shrinkage was the highest in all three provinces, and that of severe shrinkage was the lowest. The intensity of city shrinkage notably increased in S2: 56 cities exhibited severe shrinkage, 71 cities exhibited medium shrinkage and only 35 cities exhibited slight shrinkage. Severely shrinking cities were mainly located in the northern part of NEC. Among the different provinces ( Figure  10c,d), the shrinking degree were similar in three provinces in S1. In S2, the Liaoning Province displayed the highest proportion of non-shrinking cities, and the least proportion of severe shrinking cities among the three provinces, accounting for 13.2% and 22.4%

Degree of Shrinking Cities
The degree of city shrinkage in the NEC in the two periods is presented in Figure 10a,b. It is obvious that the extent of shrinkage became more severe from S1 to S2. During S1, most (77.5%) of the shrinking cities in S2 displayed medium shrinkage, and slightly and severe shrinking cities accounted for 6% and 16.7% of all cities, respectively. The severe shrinking cities were mainly located in Jilin Province. Figure 10c displays the shrinkage degree at the provincial level: the proportion of medium shrinkage was the highest in all three provinces, and that of severe shrinkage was the lowest. The intensity of city shrinkage notably increased in S2: 56 cities exhibited severe shrinkage, 71 cities exhibited medium shrinkage and only 35 cities exhibited slight shrinkage. Severely shrinking cities were mainly located in the northern part of NEC. Among the different provinces (Figure 10c,d), the shrinking degree were similar in three provinces in S1. In S2, the Liaoning Province displayed the highest proportion of non-shrinking cities, and the least proportion of severe shrinking cities among the three provinces, accounting for 13.2% and 22.4% of all the cities. The shrinking conditions in Heilongjiang Province were opposite those in Liaoning Province. Over 90% of cities were shrinking, and 38.5% were severely shrinking, the most in the three province

The Strengths and Limitations of the Proposed Method
Compared with existing methods, the merits of our proposed framework are threefold. First, consistent long-term time series of NTL data from 1996 to 2018 were constructed. The lack of long-term NTL time series is a pivotal bottleneck restricting studies of shrinking cities based on NTL data. Some recent studies have made some attempts to apply long-term time series NTL data to shrinking cities [45]. Continuous time series of data can avoid the instabilities associated with single temporal observations caused by multiple factors, and thus promote the robustness of the results. In addition, studies have shown that city shrinkage is a continuous process [23], and long-term time series provide a basis for studies of such complex and progressive issues. Most existing studies were not able to reveal the long-term trends and continuity changes related to city shrinkage due to the use of short time series [34,35]. By applying our method, the overall trend was successfully tracked, and the characteristics of city shrinkage in the early stage were detected. Moreover, the stages of city shrinkage can be divided based on consistent long-term time series. Thus, the dynamic transformations between different stages were analyzed, and the driving forces behind these shifts could be revealed. Second, to reveal the complex phenomenon of city shrinkage, NTL data and population data were integrated to establish a comprehensive indicator different from the single or basic indicators used in previous studies [34,35]. Third, our NTL index is more robust, as they were built based on representative NTL features that were extracted by the RF model.
Although this method has many advantages, there are still some aspects that can be improved. For instance, some studies pointed out that NTL data may not be a good proxy

The Strengths and Limitations of the Proposed Method
Compared with existing methods, the merits of our proposed framework are threefold. First, consistent long-term time series of NTL data from 1996 to 2018 were constructed. The lack of long-term NTL time series is a pivotal bottleneck restricting studies of shrinking cities based on NTL data. Some recent studies have made some attempts to apply long-term time series NTL data to shrinking cities [45]. Continuous time series of data can avoid the instabilities associated with single temporal observations caused by multiple factors, and thus promote the robustness of the results. In addition, studies have shown that city shrinkage is a continuous process [23], and long-term time series provide a basis for studies of such complex and progressive issues. Most existing studies were not able to reveal the long-term trends and continuity changes related to city shrinkage due to the use of short time series [34,35]. By applying our method, the overall trend was successfully tracked, and the characteristics of city shrinkage in the early stage were detected. Moreover, the stages of city shrinkage can be divided based on consistent long-term time series. Thus, the dynamic transformations between different stages were analyzed, and the driving forces behind these shifts could be revealed. Second, to reveal the complex phenomenon of city shrinkage, NTL data and population data were integrated to establish a comprehensive indicator different from the single or basic indicators used in previous studies [34,35].
Third, our NTL index is more robust, as they were built based on representative NTL features that were extracted by the RF model.
Although this method has many advantages, there are still some aspects that can be improved. For instance, some studies pointed out that NTL data may not be a good proxy for time series measures [46,47]. We acknowledge that there are some defects in using NTL for change detection. However, we did not assert that NTL only represents economic activities, since it was regarded more as a comprehensive indicator of city shrinkage. In addition, the strength and the scope of the NTL are both important dimensions in identifying shrinking cities, and NTL is sensitive to spatial changes [48][49][50][51]. In this study, we considered both the average DN value and the light pixel area when calculating NI so that the change of spatial range and the strength variation have all been taken into consideration. We hope we can further study this topic in the future. It is reported that census data from the China City Statistical Yearbook have many problems [52]. We admit that these data have multiple defects, but they have been widely used in shrinking city studies including studies in NEC [12][13][14][53][54][55][56], and have obtained many valuable research results. Moreover, since this study is more from a long time series, large scale and macro point of view to evaluate the city shrinkage, obtaining large numbers of fine-scale confidential population data is not necessary and not feasible. Thus, we think these data are capable to be used as the auxiliary population data in this study. If more accurate or detail population grid appear in the future, we will further improve this method. S1 and S2 were divided based on the intuitive judgement of data and the support of auxiliary information; thus, the accuracy of stage division could potentially be improved. Moreover, the division of stages was not detailed enough. We may apply more scientific methods or complex models such as time series breaking point detection method in the future study to solve these problems. When choosing the training set before applying the RF model, because a uniform standard for shrinking cities does not currently exist, we selected ten typical shrinking cities from Long's study [12]. This caused the credibility of our results to be influenced by Long's study. Similarly, there is currently no recognized validation dataset, which makes it difficult to verify the identification results. Previous results are not suitable for verification in this study because of the different study areas and study periods. Therefore, we merely assessed the consistency of the SCI with economic changes to verify the results to some extent. In the future, if a recognized shrinking city dataset is established, we will further verify our method by the validation dataset. In addition, we will try to find a more scientific method to verify the results by other perspectives, such as applying social media big data as validation data.

The Asynchrony of NTL and Population Changes
The most significant finding in this study is the asynchrony of NTL and population changes. This phenomenon was first observed in S1. A total of 41% of the cities in S1 were P-shrinkage or N-shrinkage cities, indicating that the changes in NI and PI were not consistent. In S2, the asynchrony of the two factors greatly increased and exhibited a significant pattern. A total of 133 cities in S2 were identified as P-shrinkage cities, accounting for 72% of all cities. These cities all experienced a decrease in population and an increase in NTL at the same time. The findings above, especially the shrinkage pattern detected in S2, are interesting. More than half of the cities in the region experienced population loss, but their NTL were still increasing. This finding contradicts intuition.
To explain these findings, the exact representation of the NTL data used in this study needs to be clarified. NTL data are normally used to reflect the intensity of human activities. Except for showing correspondence with demographic data, NTL data are highly correlated with economic activities. Several studies have used NTL data to indicate the economic activities of cities [27,37]. Since the changes in NTL data were not consistent with the changes in population data in our study, we inferred that the NI in our study possibly reflects more about the changes in city economic activities. To verify this hypothesis, we adopted linear regression for the NI with the GI and PI. The results are presented in Figure 11. The R 2 value for the GI and NI is 0.67, and the R 2 value for the PI and NI is 0.57. Thus, compared to population, the NI in our study has a higher correlation with economic activities. This finding suggests that focus should be placed on the scope of application when using NTL data. Some previous studies believed that population data can be directly replaced by NTL data. However, as the results of this study show, using a city shrinkage study as an example, the results based on NTL data do not always exhibit consistency with results based on population data. According to different research contents, huge differences may exist between the results derived from NTL data and population data.

Current Situation and Policy Implications
As the results in S2 show, approximately 85% of the cities in NEC are suffering from population decline, and the annual slop change shows that this decline is accelerating. There is a high probability that the population will continue to decline in the coming years. Unlike the pessimistic population forecast, in S2, 84% of the cities in NEC showed increasing trends associated with economic activities. Thus, we can assume that the period of simultaneous declines in the economy and population in NEC is over. The population is expected to continuously decrease, and the economy is projected to continuously develop; these conditions will be the new norm for most cities in NEC.
Although the economy in some cities has recovered, the city shrinkage situation in NEC is still grim. On the one hand, the continuous loss of population is still accelerating, which may cause potential problems. On the other hand, some unreasonable policies still exist, which may worsen the shrinkage situation in NEC. The main problems are currently as follows: (1) High-quality talent may play an important role in current urban development. However, some cities in NEC only formulated policies to increase the number of residents and have ignored the proportion of high-quality talent. (2) Due to the old age of the cities in NEC, urban infrastructures are commonly poor. This issue has seriously affected the convenience of living and further reduced people's willingness to stay. (3) Although the population has been declining for many years, some local governments consistently focused on rapid land development [59]. This unreasonable urban development strategy has further exacerbated the problem of city shrinkage. Table 2 shows the human settlement expansion rates in shrinking cities in NEC. High-speed expansion means that the land expansion rate of cities is faster than the average rate in China, while low-speed expansion means the opposite. Many NP-shrinkage cities in the two stages were also highspeed expansion cities, and this phenomenon is worthy of further assessment.  The uncoupling of economic activities and population growth in shrinking cities has also been noted in past studies. Bartholomae et al. [22] found that some German cities were experiencing population declines and increases in economic activities at the same time. Du et al. [10] showed that one-third of Chinese cities displayed the paradox of economic vitality growth and population decline. Shan also observed this phenomenon in the middle reaches of the Yangtze River. Our study results confirmed the existence of this phenomenon based on long-term continuous time series for the first time. Many hypotheses have been proposed by researchers to explain the cause of this phenomenon. Long et al. [12] theorized that economic decline in a shrinking city always lags behind population loss and suggested that population decline and economic growth uncoupling is a temporary phenomenon that happens in the early stage of city shrinkage. Jin et al. [57] proposed that economic development in many cities does not depend on an intensive labor force, and thus, the impact of population loss on the city economy is not very significant. Bartholomae et al. [22] proposed the concept of smart increasing; by improving the industrial structure and increasing the production capacity, city economic activities can be maintained during periods of population loss. All these theories have distinct rationality and are applicable in different situations. Additionally, these theories can be combined to explain the phenomenon observed in NEC.
According to the different causes, the P-shrinkage cities in S2 were divided into three categories, namely, economic revival cities, early-stage shrinking cities and populationinsensitive cities. The detailed explanation is as follows: (1) Economic revival cities were N-shrinkage cities or NP-shrinkage cities in S1, accounting for 33.1% of all cities. The economies of these cities underwent a recovery process. Since 2005, the Chinese government began its plan to revitalize Northeast China [58], and local governments in NEC cities actively responded to the plan. The transition period from S1 to S2 coincided with the implementation period of the "Revitalization of Northeast China" policy. In this period, local governments implemented many practical measures to improve economic activities in cities. Through industrial structure transformation, these cities gradually became smartly growing cities, and their labor-intensive industries and low-efficiency industries were replaced by high-tech industries and high-efficiency industries; therefore, a reduced population could generate a sufficient production capacity and economic activities.
(2) Early-stage shrinking cities were non-shrinkage cities in S1, accounting for 43.6% of all cities. Population loss began in these cities in S2; thus, these cities were in the early stages of city shrinkage. Population loss had not yet led to a further economic recession.
(3) Population-insensitive cities were P-shrinkage cities in S1, accounting for 23.3% of all cities. These cities steadily experienced population loss and economic growth throughout the study period and functioned differently from other cities in the urban system. The economies in these cities do not depend on labor; thus, the effect of population loss was insignificant.

Current Situation and Policy Implications
As the results in S2 show, approximately 85% of the cities in NEC are suffering from population decline, and the annual slop change shows that this decline is accelerating. There is a high probability that the population will continue to decline in the coming years. Unlike the pessimistic population forecast, in S2, 84% of the cities in NEC showed increasing trends associated with economic activities. Thus, we can assume that the period of simultaneous declines in the economy and population in NEC is over. The population is expected to continuously decrease, and the economy is projected to continuously develop; these conditions will be the new norm for most cities in NEC.
Although the economy in some cities has recovered, the city shrinkage situation in NEC is still grim. On the one hand, the continuous loss of population is still accelerating, which may cause potential problems. On the other hand, some unreasonable policies still exist, which may worsen the shrinkage situation in NEC. The main problems are currently as follows: (1) High-quality talent may play an important role in current urban development. However, some cities in NEC only formulated policies to increase the number of residents and have ignored the proportion of high-quality talent. (2) Due to the old age of the cities in NEC, urban infrastructures are commonly poor. This issue has seriously affected the convenience of living and further reduced people's willingness to stay. (3) Although the population has been declining for many years, some local governments consistently focused on rapid land development [59]. This unreasonable urban development strategy has further exacerbated the problem of city shrinkage. Table 2 shows the human settlement expansion rates in shrinking cities in NEC. High-speed expansion means that the land expansion rate of cities is faster than the average rate in China, while low-speed expansion means the opposite. Many NP-shrinkage cities in the two stages were also high-speed expansion cities, and this phenomenon is worthy of further assessment. City shrinkage is not necessarily a bad thing for urban development and management; instead, it depends on whether reasonable strategies can be adopted. In view of the city shrinkage problem in NEC, we put forward the following suggestions: (1) Change the traditional concept of urban development that focus more on quantity rather than quality is a priority for local governments, and the concept of efficient shrinkage needs to be introduced. Efficient shrinkage provides a path for solving city shrinkage problems, as noted by Popper et al. [60], who emphasized the importance of improvements in inventory, efficiency and development quality for shrinking cities. (2) Improve the management of human resources, and introduce high-tech talent. For example, Hegang has achieved great results on attracting high-tech talent through high-welfare policies, and economic activities in the city have therefore continually increased. (3) Adjust the industrial structure of traditional industrial cities, and use new technology to replace labor-intensive industries. For instance, by prohibiting new expansion of traditional industries and providing reasonable assistance to withdraw from low-efficiency, low-tech and high-energy-consumption industries, Anshan successfully reversed its shrinking trend. (4) Redefine the function of cities, and develop new plans for cites. The experience of western countries is worthy of attention in this context. In the face of continuous city shrinkage, Buffalo, NY, redefined its functional position and provided a large number of jobs. Youngstown, OH, reformulated urban planning based on small-scale urban areas and demolished built-up areas to improve the quality of urban life [60]. Because different shrinking cities have different features, practical policies should be chosen to suit their local conditions.

Conclusions
In this study, a novel method that combines long-term time series of NTL data and population data was applied to identify and classify shrinking cities in NEC from 1996 to 2020. In detail, we applied the GWR regression model to achieve cross-sensor calibration, and obtained long-term continuous NTL time series. Moreover, we generated NI by RF model and PI by the slope of population to indicate NTL and population change, respectively. We also combined the two indices into SCI, and identified and classified shrinking cities in NEC. The results show that most of the cities in NEC have been shrinking for the past 25 years. Moreover, the shrinkage patterns in the two stages were quite different, and economic and population changes were not synchronous. From S1 to S2, the shrinkage situation worsened, and the heterogeneity of shrinkage types increased. In S2, most of the cities experienced accelerated population decline and increases in economic activities, suggesting that governments need to maintain implementing effective measures and take the route of efficient shrinkage.
In future works, multiple sources datasets, such as industrial data, traffic data, ecological data and other types of data that are possibility related to city shrinkage, can be used to further improve the robustness of the identification of city shrinkage in this method. With multiple fine-scale auxiliary data, this method can also be improved from region scale into grid scale. In addition, this method should be further examined with other regions and cities, and establish a credible shrinking city dataset. Moreover, it is meaningful to further study the stage division of shrinking cities, more robust method should be applied to have a more detailed and accurate stage division result, and thus the driving forces of different stages can be revealed.