A Scheme for the Long-Term Monitoring of Impervious−Relevant Land Disturbances Using High Frequency Landsat Archives and the Google Earth Engine

Abstract

Impervious surfaces are commonly acknowledged as major components of human settlements. The expansion of impervious surfaces could lead to a series of human−dominated environmental and ecological issues. Tracing impervious surface dynamics at a finer temporal−spatial scale is a critical way to better understand the increasingly human-dominated system of Earth. In this study, we put forward a new scheme to conduct long-term monitoring of impervious−relevant land disturbances using high frequency Landsat archives and the Google Earth Engine (GEE). First, the developed region was identified using a classification-based approach. Then, the GEE-version LandTrendr (Landsat-based detection of Trends in Disturbance and Recovery) was used to detect land disturbances, characterizing the conversion from vegetation to impervious surfaces. Finally, the actual disturbance areas within the developed regions were derived and quantitatively evaluated. A case study was conducted to detect impervious surface dynamics in Nanjing, China, from 1988 to 2018. Results show that our scheme can efficiently monitor impervious surface dynamics at yearly intervals with good accuracy. The overall accuracy (OA) of the classification results for 1988 and 2018 are 95.86% and 94.14%. Based on temporal−spatial accuracy assessments of the final detection result, the temporal accuracy is 90.75%, and the average detection time deviation is −1.28 a. The OA, precision, and recall of the sampling inspection, respectively, are 84.34%, 85.43%, and 96.37%. This scheme provides new insights into capturing the expansion of impervious−relevant land disturbances with high frequency Landsat archives in an efficient way.

1. Introduction

Reflecting a substantial demographic transformation from rural to urban, urbanization has become a worldwide phenomenon in recent decades [1]. As one of the predominant consequences of global urbanization, the Earth’s terrestrial surface has been experiencing an accelerated conversion to urban areas to meet the increasing demands of human activities [2,3,4,5,6]. China, which is no exception, has been undergoing rapid urbanization and associated land cover changes since the 1980s [6].

Impervious surfaces are commonly acknowledged to be major components of human settlements [6]. Impervious surfaces are defined as artificial materials that prevent the infiltration of water into the soil, such as roads, rooftops, sidewalks, and patios [6,7,8]. Impervious surfaces are significant factors impacting the ecosystem [6,7,8,9]. The mounting expansion of impervious surfaces could lead to a series of human-dominated environmental and ecological issues [9,10,11]. For example, the increasing threatens to biodiversity and ecosystem productivity [3], the aggravation of the urban heat island effect [12], the changes in the surface’s microclimate, and hydrology [13]. Therefore, monitoring impervious surface expansion is significant to understand the impacts of urbanization and associated land cover changes on the Earth’s ecosystem [3,9,10,14]. Moreover, the datasets documenting impervious surface dynamics both statistically and temporal−spatially are highly essential for policy making, environmental change research, demographic research, resource management, and urban sustainability [1,6,15,16,17].

The advancement of satellite-based remote sensing has been a remarkable revolution in the Earth observations. Almost 50 years of archived data made the Landsat program earn the record of the longest running terrestrial satellite [18,19]. Importantly, the free and open access of all Landsat archives in 2008 provided a valuable opportunity to capture the geospatial information at a fine spatial resolution [18,19,20,21]. With progress in automated atmospheric correction algorithms and cloud detection algorithms [19,22], tracing target dynamics using Landsat time series has proven to outperform methods that use a few scenes and discrete dates [16,19,23].

Generally, analytical approaches towards detecting changes with Landsat time series can be roughly divided into three aspects: classification, thresholding-based extraction, and spectral trajectory-based detection.

• Classification performance is dependent on advanced classification algorithms. Over the past decades, supervised classifiers have especially grown in utility to accommodate complex feature [24,25], such as Random Forests (RF) [26]. These classifiers have been widely used for annual classification to generate time series maps of the informed land covers, in order to determine both the states and trends of the changes [6,11,14,24,27,28,29]. However, sensitivity to training data is a common issue for most supervised classifiers [24,25]. On the other hand, post-classification approaches are indispensable for improving spatial-temporal consistency of the classification results [11,28,29].
• Thresholding-based extraction has been widely employed in identifying specific land covers in the time series [23]. Significant deviations between the target and other land covers are the key to detecting the target. Thus, images need to be transformed into the dimension that is sensitive to the target [21,23], such as the NUACI (Normalized Urban Areas Composite Index) [21], which was proposed for the multi-temporal mapping of global urban lands. Thresholding-based extraction is indeed convenient but highly dependent on the predefined threshold [23].
• Spectral trajectory-based detection is a milestone in monitoring changes with Landsat time series [15,19,23]. Detection approaches using high temporal frequency satellite data have become a frontier of research, enabling a more nuanced understanding of changes on the Earth [15,19], for instance, the vegetation change tracker (VCT) [30], Landsat-based detection of Trends in Disturbance and Recovery (LandTrendr) [31], an algorithm based on statistical quality control charts [32], Continuous Change Detection and Classification (CCDC) [33], Continuous Subpixel Monitoring (CSM) [10], and Continuous monitoring of Land Disturbance (COLD) [34]. These spectral trajectory-based methods are well adopted in diverse thematic domains. However, implementing these methods often requires rigorous data pre-processing, affordable computation, and time-consuming operations, which limits their more extensive applications.

Considerable tools and platforms have been applied to process diverse remote sensing data (for instance, ENVI and ERDAS). At present, the platform offered by the Google Earth Engine (GEE) (https://earthengine.google.com/) has been developed to cope with complex data and algorithms. GEE is a cloud-based platform with the ability to execute various algorithms and access petabyte-scale geospatial datasets [35]. An increasing number of studies have adopted GEE to research issues of the Earth’s terrestrial surface [6,21,36,37,38,39].

Land disturbances characterizing the conversion from vegetation to impervious surfaces have frequently occurred in most areas of southern China. Such land disturbances are a breakthrough that reflects the expansion of impervious surfaces. Because of its capability to target both discrete events and gradual trends [40,41], LandTrendr is known to be an effective spectral trajectory-based algorithm in detecting the changes associated with vegetation loss or gain. After LandTrendr was ported to the GEE platform, the GEE-based LandTrendr (LT-GEE) [37] could shorten the computing time from days to minutes. LT-GEE also simplifies onerous data management and image-preprocessing by directly accessing the geospatial datasets in GEE. These advantages of the cloud-based platform presented an opportunity to fast and effectively monitor impervious–relevant land disturbances over a long time span. However, employing only detection methods would generate a certain amount of noises that do not belong to actual changes. Thus, performing a detection method in the identified regions associated with relevant land conversions could offer new insights to improve detection performance and efficiency.

In this study, we suggest a new scheme to rapidly implement the long-term monitoring of impervious–relevant land disturbances at a regional scale and annual bias. This scheme aims to both improve the detection performance in heterogeneous backgrounds and reduce noises by combining spatial constraints with a temporally-dense detection method. Therefore, we sought to (1) establish the developed region using a classification-based approach; (2) employ LT-GEE and high temporal frequency Landsat archives, to capture the impervious–relevant land disturbances at a fine spatial resolution and annual bias; (3) derive the actual disturbances within the developed region and quantitatively evaluate the temporal–spatial accuracy of the results. A case study was conducted in Nanjing, China, to detect the impervious surface dynamics from 1988 to 2018.

2. Study Area and Data

2.1. Study Area

The capital city of Jiangsu Province, Nanjing was selected as our study area. Covering an area of 6587.02 km² [42], Nanjing (31°14′–32°37′N, 118°22′–119°14′E) is located in eastern China and the middle area of the lower reaches of the Yangtze River (Figure 1a). The administrative boundary of Nanjing contains 11 districts, which are Xuanwu, Qinhuai, Jianye, Gulou, Pukou, Qixia, Yuhuatai, Jiangning, Liuhe, Lishui and Gaochun (Figure 1b). Nanjing belongs to a subtropical and humid climate. The main region of Nanjing is covered by one scene (path/row: 120/038) in the World References System-2 (WRS-2), and the remaining areas are covered by three scenes (path/row: 120/037, 121/038, 121/037) (Figure 1a).

Nanjing is an important gateway to the regions surrounding the Yangtze River Delta. Nanjing is a typical city in the Yangtze River Delta Urban Agglomeration that has experienced rapid growths in its population and economy over the past decades. The resident population of Nanjing in 2017 was 8.335 million, reaching almost 1.7 times the population in 1988. Economically, the regional gross domestic product (GDP) increased to RMB 1171.51 billion in 2017, accounting for 1.4% of the national GDP [42]. During this period, the urbanization of Nanjing has been stimulated by population growth, economic development, and industrialization. To meet the increasing demands of human activities, a great number of lands have been converted from agricultural or natural lands to built-up areas covered by artificially impervious materials. The drastic expansion of impervious surfaces has become a noteworthy issue due to their impact on the environment and ecosystem.

Conducting our scheme in Nanjing will help understand the characteristics and mechanisms of impervious surface expansion and land conversions in the Yangtze River Delta urban agglomeration. Moreover, our local knowledge is relatively sufficient to support data collection and result analysis.

2.2. Data

2.2.1. Landsat Data

We accessed the Landsat Collection 1 Tier 1 products from the United States Geological Survey (USGS) (https://www.usgs.gov/land-resources/nli/landsat) via the GEE platform, which includes the USGS Landsat Collection 1 Tier 1 Raw Scenes and the USGS Landsat Surface Reflectance Tier 1. The Landsat scenes with the highest data quality are placed into Tier 1 and are spatially referenced using the World Reference System-2 (WRS-2) in the UTM projection [43]. In GEE, images from a single sensor were regrouped into an individual “ImageCollection” to feasibly filter specific images from the multi-petabyte geospatial datasets [35].

1988 was chosen as the initial stage because of the poor georegistration of images from 1985 and 1986, and insufficient observations with lower cloud coverages in 1984 and 1987 (Figure 2). All available Raw Scenes in 1988 and 2018 were applied to establish the developed region using classification-based approaches. These data are the products processed by georegistration, and calibrated with at-sensor radiance. The ImageCollection IDs are LANDSAT/LT05/C01/T1 and LANDSAT/LC08/C01/T1.

Surface Reflectance products from 1988 to 2018 were adopted to detect impervious–relevant land disturbances using LT-GEE. All data are the atmospherically corrected surface reflectance from the Landsat 5 TM (Thematic Mapper), Landsat 7 ETM+ (Enhanced Thematic Mapper Plus), and Landsat 8 OLI (Operational Land Imager) sensors. The atmospheric correction approaches include LEDAPS (Landsat Ecosystem Disturbance Adaptive Processing System) for TM/ETM+ [44] and LaSRC (Landsat 8 Surface Reflectance Code) for OLI [45].

2.2.2. Auxiliary Data

The auxiliary data include the ground truth samples for classification and the reference data for the temporal–spatial accuracy assessment.

(1) Land cover category and samples for classification:

Based on the Level 1 Type in the classification system proposed by Gong et al. [27] and the geographical environment of Nanjing, the main land covers in the study area were classified as bareland, cropland, forest, impervious, and water.

The ground truth samples were collected via visual interpretation, aided by very high resolution (VHR) images from Google Earth, historical Landsat images, and prior knowledge. A total sum of 2571, 2733, and 2832 samples were yielded for 1988, 2010, and 2018 using a stable collection approach [9,27,28,29,46] (

Table 1. Details of samples for classification.

Year	Bareland	Cropland	Forest	Impervious	Water	Sum
2018	123	1494	469	598	148	2832
2010	105	1487	468	528	145	2733
1988	157	1556	465	247	146	2571

). For each year, 70% of the samples were used to train the classifier, and the remaining 30% were applied for verification.

(2) Reference data for the temporal–spatial accuracy assessment:

Temporal samples and a spatial reference were applied to evaluate the temporal and spatial accuracy of our change detection results. In detail, based on the stratified random sampling strategy, a total of 227 temporal samples were collected throughout the historical VHR images from Google Earth. All samples were collected from the developed areas that were identified by the approach in Section 3.1 and were labeled with the change time converted from other land covers to impervious surfaces (mainly from 2005 to 2017). Only the first change was used if multiple changes occurred within one sample spot. Based on the classification result in 2010 and the VHR images from Google Earth, the impervious surface in a rectangle measuring 175.2 km² (left upper: 118°43′55.67′′E, 32°1′54.04′′N; right bottom: 118°53′8.1′′E, 31°55′8.52′′N) was manually rechecked as a spatial reference.

2.2.3. Other Open-Access Products

Two kinds of open-access products were collected to compare against our detection results. These products include diverse global land cover maps and three statistical yearbooks. However, most of the global land cover maps were only mapped for one epoch. We collected these products as much as possible, including the Finer Resolution Observation and Monitoring of Global Land Cover (FROM-GLC) series products [27,36,47,48,49], the high-resolution multi-temporal maps of global urban land [21], and the Global Human Settlement Layers (GHSL) [50]. Moreover, the statistical data of built-up areas were extracted from three yearbooks [42,51,52]. Product details are listed in

Table A1. Details of the open-accessing products used in this study.

Data	Year	Data Structures	Method	Coverage	Reference
FROM-GLC10	2017	Raster (10 m)	Pixel-based Classification	Global	[36]
FROM-GLC Version2	2015	Raster (30 m)	Pixel-based Classification	Global	[49]
FROM-GLC	2010	Raster (30 m)	Pixel-based Classification	Global	[27]
FROM-GLC-agg	2010	Raster (30 m)	Multisource data aggregation	Global	[48]
FROM-GLC-seg	2010	Raster (30 m)	A segmentation-based approach	Global	[47]
NUACI-based mapping	1990, 1995, 2000, 2005, 2010, 2015	Raster (30 m)	NUACI-based classification	Global	[21]
GHSL	1975, 1990, 2000, 2015	Raster (30 m)	A classification based on symbolic machine learning	Global	[50]
China city statistical yearbook	1985–2010	Records	Statistical methods	China	[51]
Jiangsu statistical yearbook	2011–2018	Records	Statistical methods	Jiangsu	[52]
Nanjing statistical yearbook	1978–2018	Records	Statistical methods	Nanjing	[42]

.

3. Methodology

As shown in Figure 3, the scheme we developed has two main parts, which include (1) the establishment of the developed region using a classification-based approach, and (2) long-term monitoring of the impervious–relevant land disturbance. Specifically, the developed region was first extracted to build spatial constraints. Then, a trajectory-based detection algorithm was applied to detect when and where impervious–relevant land disturbances occurred in the developed region. More explanations of the scheme are detailed in the following sections.

3.1. Establishment of the Developed Region Using a Classification-Based Approach

Figure 4 shows the flow chart of the establishment of the developed region using a classification-based approach. The developed region was extracted based on the layers of impervious surfaces in the initial and terminal stages. In this study, we designed a classification framework (shown in Figure 4) to obtain the land cover maps with accuracy in each epoch. First, a pixel-based compositing approach was adopted to reconstruct seasonal cloud-free composites. Second, image features were extracted from seasonal composites. Then, the land cover classifications were conducted using the Random Forests (RF) classifier. Finally, the land cover maps for 1988, 2010, and 2018 were generated. The classification results for 1988 and 2018 were adopted to extract the developed region, and the classification result for 2010 was used to obtain the spatial reference mentioned in Section 2.2.2.

3.1.1. Pixel-Based Compositing

A best-available-pixel (BAP) compositing method provided by GEE’s algorithm library was employed to generate the least cloud-free composites for each season in 1988, 2010, and 2018. Pixel-based compositing methods can provide a new analysis paradigm instead of depending on a single scene [53]. We used a function named simpleComposite to compute and return a top-of-atmosphere (TOA) reflectance [54] composite from the given collection of raw scenes. The returned TOA is the composite assigned with the least-cloud-covered pixels.

3.1.2. Feature Extraction and Classification Using RF

Six spectral bands and three spectral indices were extracted from each seasonal composite. These bands and indices include green, blue, red, NIR (Near Infrared), SWIR 1 (Short Wavelength Infrared 1), SWIR 2 (Short Wavelength Infrared 2), the Normalized Difference Vegetation Index (NDVI) [55], the normalized difference moisture index (NDMI) [56], and the modified Normalized Difference Water Index (mNDWI) [57]. The image features extracted from the four seasonal composites constructed annual image feature stacks, so a total of 36 image features were used for classification in each year. Calculations of the indices are specified below:

$$\mathrm{NDVI}=\frac{{\rho }_{\mathrm{NIR}}-{\rho }_{\mathrm{Red}}}{{\rho }_{\mathrm{NIR}}+{\rho }_{\mathrm{Red}}}$$

(1)

$$\mathrm{NDMI}=\frac{{\rho }_{\mathrm{NIR}}-{\rho }_{\mathrm{SWIR}1}}{{\rho }_{\mathrm{NIR}}+{\rho }_{\mathrm{SWIR}1}}$$

(2)

$$\mathrm{mNDWI}=\frac{{\rho }_{\mathrm{Green}}-{\rho }_{\mathrm{SWIR}1}}{{\rho }_{\mathrm{Green}}+{\rho }_{\mathrm{SWIR}1}}$$

(3)

where${\rho }_{\mathrm{NIR}}$,${\rho }_{\mathrm{Red}}$,${\rho }_{\mathrm{Green}}$, and ${\rho }_{\mathrm{SWIR}1}$ represent the near infrared, red, green, and short wave infrared (B6: 1.57–1.65 $\mathsf{\mu }\mathrm{m}$ for OLI, B5: 1.55–1.75 $\mathsf{\mu }\mathrm{m}$ for TM/ETM+) bands.

The Random Forests (RF) [26] classifier was adopted for classification. RF has been widely applied in land cover mapping by using on-the-fly images with reliable accuracy [24,25,58]. RF is also a robust classifier that can use high dimensional data involved in processing [27,59]. In this study, the number of decision trees to be generated ($Ntree$) was set as 1000, the number of variables to be selected and tested for the best split of each node when growing the trees ($Mtry$) was set as the square root of the number of the input variables and other parameters were set as defaults.

3.1.3. Establishment of the Developed Region

The developed region layer was derived using the classification results from 1988 and 2018. Based on the general assumption that the expansion of an impervious surface is irreversible [11,17,28,60], we inherited the hypothesis that conversions from impervious surfaces to other land covers are impossible. We defined developed regions as regions that have experienced conversions from natural land to impervious surfaces. Therefore, the developed region equals the impervious surfaces in 2018 subtracted the intersected parts in both 1988 and 2018 (Figure 4). No additional post-processes were applied to the results.

3.2. Monitoring of the Impervious–Relevant Land Disturbance

3.2.1. LT-GEE Algorithm

LT-GEE is a new version of the LandTrendr (hereafter LT) algorithm ported to the GEE platform [31,37]. LT is a spectral trajectory-based change detection algorithm. It is capable of capturing change events and smoothing entire trends from the annual Landsat time series stacks. Basically, LT needs data pre-processing and trajectory segmentation [31,41]. In detail, LT first extracts the spectral trajectories at the pixel scale from a pixel’s spectral history, like a band or an index. Then, the temporal segment algorithm captures the broad features of the trajectory by modeling every pixel’s spectral time series as a sequence of straight line segments [31]. Then, it identifies the breakpoints separating periods of durable change or stability in the spectral trajectory. Finally, LT records the year when changes occurred. Figure 5 illustrates how LT finds a breakpoint (vertex) using the trajectory segment. LT was initially implemented in IDL (Interactive Data Language, hereafter LT-IDL) (https://github.com/KennedyResearch/LandTrendr-2012). Onerous data management and image pre-processing are almost eliminated by LT-GEE with the help of GEE (https://github.com/eMapR/LT-GEE).

Implementing LT-GEE requires two parts:

• The construction of an annual image collection:

LT-GEE applies USGS Landsat Surface Reflectance Tier 1 to generate an annual image collection using the medoid compositing approach [31]. For a specific pixel in the image stacks during a given period, the medoid is the value for a given band that is numerically closest to the median of all corresponding pixels [37]. Meanwhile, to solve the spectral property differences between OLI (bands 2–6, 7) and TM/ETM+ (bands 1–5, 7), the OLI image is transformed into the properties of the TM/ETM+ using the method proposed by Roy et al. [61] (hereafter the Thematic Mapper-equivalent Band 1–5 and 7). Moreover, clouds or cloud shadows in each required Surface Reflectance image are masked out using CFMASK [22,62].

• Parameter setting and trajectory segmentation:

After data pre-processing, LT-GEE will first calculate the input of the spectral band or index to extract the spectral trajectory. In this annual image collection, the calculated input becomes the first band to be segmented, and the annual fitted-to-vertex (FTV) data for subsequent bands are also generated [31,37,41]. Controlling parameter is the key to trajectory segmentation. LT-GEE’s parameters are detailed in

Table A2. Details of LT-GEE parameters.

Parameter	Type	Definition	Default
startYear	Integer	The initial year of detection duration	1984
endYear	Integer	The terminal year of detection duration	2017
startDay	Integer	The time of initial year	0620
endDay	Integer	The time of terminal year	0910
index	String	The first band in the image collection and generate annual fitted-to-vertex (FTV) data for each subsequent band	NBR
maxSegments	Integer	Maximum number of segments to be fitted on the time series	6
spikeThreshold	Float	Threshold for dampening the spikes (1.0 means no dampening)	0.9
vertexCountOvershoot	Integer	The initial model can overshoot the maxSegments + 1 vertices by this amount. Later, it will be pruned down to maxSegments + 1	3
preventOneYearRecovery	Boolean	Prevent segments that represent one year recoveries	true
recoveryThreshold	Float	If a segment has a recovery rate faster than 1/recoveryThreshold (in years), then the segment is disallowed	0.25
pvalThreshold	Float	If the p-value of the fitted model exceeds this threshold, then the current model is discarded and another one is fitted using the Levenberg-Marquardt optimizer	0.05
bestModelProportion	Float	Takes the model with most vertices that has a p-value that is at most this proportion away from the model with lowest p-value	0.75
minObservationsNeeded	Integer	Minimum observations needed to perform output fitting	6
timeSeries	ImageCollection	Collection from which to extract trends (it is assumed that each image in the collection represents one year). The first band is used to find breakpoints, and all subsequent bands are fitted using those breakpoints
delta	String	Vegetation change type, ensure that the band or index to be segmented is oriented with vegetation “loss” and “gain”	lose
sort	String	Vegetation change sort, which contains “Least”, “Newest”, “Greatest”, “Oldest”, “Fastest”, and “Slowest”	greatest
year	Boolean, Integer	Filter detection result by year of detection	checked: true, start: 1986, end: 2017
mag	Boolean, Integer, String	Filter detection result by magnitude	checked: true, value: 200, operator: ‘>’
dur	Boolean, Integer, String	Filter detection result by change event duration	checked: true, value: 4, operator: ‘<’
preval	Boolean, Integer, String	Filter detection result by pre-change spectral value	checked: true, value: 300, operator: ‘>’
mmu	Boolean, Integer	Filter detection result by minimum mapping patch unit	checked: true, value: 11

.

3.2.2. Change Detection Framework for Impervious–relevant Land Disturbances

The change detection framework of the impervious–relevant land disturbance shown in Figure 6 includes the input of the spectral band or index for detection, the parameter configuration, and integration of the detection results. First, an annual image collection was constructed through the preprocessing of LT-GEE. The selected spectral band or index was calculated and input to LT-GEE for trajectory segmentation. Then, we derived the optimal parameter configurations of each input band or index for detection via parameter adjustments. Finally, the majority vote algorithm was adopted to integrate all detection outputs into a final result.

(1) Inputs for spectral bands or indices for detection:

LT-GEE provides multiple choices to construct a historical spectral trajectory as an input [37]. We adopted five individual inputs as the first band in the annual image collection to extract and segment the trajectory. These inputs included Thematic Mapper-equivalent Band 3 and 7 (hereafter B3 and B7), the Normalized Burn Ratio (NBR) [63], NDVI [55], and NDMI [56].

$$\mathrm{NBR}=\frac{{\rho }_{\mathrm{NIR}}-{\rho }_{\mathrm{SWIR}2}}{{\rho }_{\mathrm{NIR}}+{\rho }_{\mathrm{SWIR}2}}$$

(4)

where${\rho }_{\mathrm{NIR}}$ and ${\rho }_{\mathrm{SWIR}2}$ present the Thematic Mapper-equivalent Band 4 and 7 after medoid compositing. Other indices’ formulas are mentioned in Section 3.1.2.

(2) Parameter configuration:

A parametric tuning approach [64] was used to obtain the optimal parameter configurations for LT-GEE. We literally divided these parameters (Table A2) into five parameter sets based on their functions: (1) controls for the detection time range (startYear, endYear, startDay, and endDay), (2) inputs for the spectral band or index (index), (3) controls for the trajectory segment performance (maxSegments, spikeThreshold, vertexCountOvershoot, preventOneYearRecovery, recoveryThreshold, pvalThreshold, bestModelProportion, and minObservationsNeeded), (4) orientation of the vegetation change tendency (delta and sort), and (5) options for pixel filtering (year, mag, dur, preval, and mmu).

Under the same detection time range and vegetation change tendency, we mainly configured the parameters in the third set with a different index. Without filtering any pixels, the configurations were performed as follows: (1) set the variation range ($\left[Pa{r}_{n.min},Pa{r}_{n.max}\right]$) and step length for each parameter. The step length was set as 0.05 for the float, or 1 for the integer, (2) create the sequences of the parameter values to be tested at equal step length intervals ($Pa{r}_n=\left\{Pa{r}_{n.min}, Pa{r}_{n.2}, Pa{r}_{n.3},\dots , Pa{r}_{n.max}\right\}$), (3) derive the corresponding detection outputs using the parameter sequence of maxSegments ($Pa{r}_1$). Other parameters were set as the default values during this processes, (4) calculate the temporal accuracy for all outputs of maxSegments ($P_1=\left\{P_{1.1}, P_{1.2}, P_{1.3},\dots , P_{1.x}\right\}$), (5) the value in $Pa{r}_1$ that yields the highest $P_1$ will be the final configuration of maxSegments. (6) According to (3)–(5), use the values of the configured parameters to participate in the configuration of the next parameter, in the order of the parameters listed in Table A2.

(3) Output and integration of the detection results:

Multiple outputs can be generated using the optimal parameter configurations for each input (index). The majority vote algorithm [62] was used to integrate all outputs into a final result, assuming that $\left\{R_1,R_2 , \dots , R_{n-1}, R_n\right\}$ is an output sequence of a single point-of-view in row $i$ and column$j$. This method allows one to check whether a majority element exists in the sequence. If so, the majority element’s value,$R_m$, is assigned to the pixel in$\left(i,j\right)$ in the final result when the number of majority elements is greater than$n/2$. Conversely, when no majority element exists, or the number of majority elements is less than$n/2$, the value of $\left(i,j\right)$ in the final result is assigned with the value in the sequence with the highest temporal accuracy.

3.3. Accuracy Assessment

3.3.1. Accuracy Assessment for Classification

Using independent validation samples, three confusion matrixes were constructed to assess the accuracy of the classification results for the initial stage in 1988, the terminal stage in 2018, and the comparable stage in 2010. The overall accuracy (OA), producer accuracy (PA) and user accuracy (UA) were calculated. The OA, PA and UA of impervious (${\mathrm{PA}}_{\mathrm{imp}}$ and${\mathrm{UA}}_{\mathrm{imp}}$) were required to exceed 90% for convincing results.

3.3.2. Accuracy Assessment for Change Detection

Assessment metrics were designed to quantify the differences between the detection results and actual changes in temporal–spatial domains.

(1) Temporal accuracy assessment

Temporal accuracy was assessed using temporal samples, assuming that $Y_T$ and$Y_D$, respectively, are the actual change time and the detection time in a spot. The detection time deviation ${\Delta }_Y$ for each spot was calculated as${\Delta }_Y=Y_D-Y_T$. The time delay between the land cover change and the actual emergence of a built-up area commonly exists in our study area. Therefore, we consider the detected change to be correct if$\left|{\Delta }_Y\right|\le 3$, based on prior observations. The temporal accuracy ($P$) refers to the proportion of correct spots (C) in all 227 temporal samples, which can be represented as$P=\frac{card C}{227}$. A higher temporal accuracy, lower average and lower standard deviation (SD) of ${\Delta }_Y$ are expected.

(2) Spatial accuracy assessment

Employing rechecked spatial reference in 2010, we preferred to apply a sampling inspection to evaluate the outputs’ performance in spatial distribution. The Precision, recall, and OA were calculated by the confusion matrix (

Table 2. The confusion matrix for the change detection results in 2010.

Pixels		Change Detection Result
		Impervious	Non-Impervious
Spatial Reference	Impervious	AA	AB
	Non-Impervious	BA	BB

).

$$\mathrm{Precision}=\frac{AA}{AA+BA}$$

(5)

$$\mathrm{Recall}=\frac{AA}{AA+AB}$$

(6)

$$\mathrm{Overall}\mathrm{accuracy}=\frac{AA+BB}{AA+BB+AB+BA}$$

(7)

4. Results

4.1. Classification Results and the Developed Region

The result of establishing the developed region is displayed in Figure 7. The accuracy assessments for the classification results are listed in

Table 3. Accuracy assessments for classification results in 1988, 2010, and 2018 (%). UA, user accuracy; PA, producer accuracy; OA, overall accuracy.

Class	1988			2018			2010
	UA	PA	OA	UA	PA	OA	UA	PA	OA
Bareland	89.47	79.07	95.86	73.08	61.29	94.14	79.17	65.52	93.36
Cropland	96.71	98.54		95.19	96.81		94.08	95.60
Forest	98.51	93.62		96.48	93.20		95.00	92.36
Impervious	90.91	93.75		92.57	93.10		93.10	92.05
Water	91.67	93.62		94.00	95.92		90.20	97.87

.

As shown in Table 3, all classification results achieved a high accuracy, with OAs over 90%. The focused impervious surface shows a relatively high accuracy, with UAs over 90% and PAs over 92%. It was noted that cropland is the most recognizable type, with UAs over 94% and PAs of over 95%, followed by forest (UAs of over 95% and PAs of over 92%) and water (UAs of over 90% and PAs of over 93%). In addition, the results illustrated there was a relatively poor identification for bareland, with a variation of UAs from 73.08% to 89.47% and PAs from 61.29% to 79.07%.

Figure 7a,b shows the land cover maps in 1988 and 2018. Dark blue polygons in Figure 7c represent the developed regions for change detection. Increasing by approximately 805.29 km², the impervious surfaces in 2018 occupied about 14.22% of the total study area. During this period, the conversions from cropland to impervious surfaces show the largest proportion, with almost 11.69% of the total study area. Forest is the second highest contributor to conversions to impervious surfaces, with about 0.25% of the study area.

4.2. Change Detection Results

4.2.1. Optimal Parameter Configurations

The optimal parameter configurations for each input (index) are listed in

Table 4. Optimal parameter configurations for the different inputs (index).

Parameter	Input Spectral Band or Index
	NBR	NDVI	NDMI	B3	B7
startYear	1988	1988	1988	1988	1988
endYear	2018	2018	2018	2018	2018
startDay	0601	0601	0601	0601	0601
endDay	0930	0930	0930	0930	0930
index	NBR	NDVI	NDMI	B3	B7
maxSegments	10	10	10	10	10
spikeThreshold	0.2	0.4	0.2	0.2	0.1
vertexCountOvershoot	3	5	15	15	3
preventOneYearRecovery	true	true	true	true	true
recoveryThreshold	0.35	0.35	0.25	0.75	0.50
pvalThreshold	0.05	0.05	0.05	0.05	0.05
bestModelProportion	0.75	0.50	0.50	0.75	0.50
minObservationsNeeded	6	6	6	6	6
timeSeries	-	-	-	-	-
delta	lose	lose	lose	lose	lose
sort	greatest	greatest	greatest	greatest	greatest
year	false	false	false	false	false
mag	false	false	false	false	false
dur	false	false	false	false	false
preval	false	false	false	false	false
mmu	false	false	false	false	false

. Local summer was adopted as the time window to detect changes. More subtle changes associated with vegetation loss could be captured due to the abundant information about vegetation from June to September. The year, mag, dur, preval, and mmu were set as “false” to maintain the integrity of the detection outputs without filtering any patch.

4.2.2. Change Detection Results and Quantitative Evaluations

Running LT-GEE with the optimal parameter configurations of different inputs, five optimal outputs were derived. These outputs were applied to generate an integrated result using the majority vote. The integrated map of the impervious–relevant land disturbances is displayed in Figure 8. Bright yellow and dark blue, respectively, represent the impervious surfaces for the initial stage in 1988 and the terminal stage in 2018. The expansions of impervious surfaces dynamics are accumulated using gradient colors from bright yellow to dark blue.

(1) Temporal accuracy

Temporal accuracy assessments are listed in

Table 5. Temporal accuracy assessments of all change detection results.

Input	Sample Number	${\mathbf{\Delta }}_{\mathit{Y}}\le 3$		${\mathbf{\Delta }}_{\mathit{Y}}\le 2$		SD	Mean
		Number of errors	Accuracy (P) (%)	Number of errors	Accuracy (P) (%)
NBR	227	21	90.75	39	82.82	3.11	−1.34
NDVI		23	89.87	41	81.94	3.15	−1.36
NDMI		30	86.78	55	75.77	3.86	−1.83
B3		38	83.26	61	73.13	3.95	−1.21
B7		26	88.55	34	85.02	3.66	−0.63
Majority Vote		21	90.75	40	82.38	3.08	−1.28

. Before all optimal outputs were integrated into one, NBR’s output had the highest temporal accuracy (P = 90.75%) and a relatively lower SD (3.11) when limited with${\Delta }_Y\le 3$. One the other hand, the B7′s output possessed the highest temporal accuracy (P = 85.02%) and the lowest average (−0.63) when limited with${\Delta }_Y\le 2$. The B3′s output had relatively poor temporal accuracy. After the majority vote, the integrated map had the same temporal accuracy as the NBR’s output (P = 90.75%), when checked by${\Delta }_Y\le 3$. However, the SD of the integrated map decreased to 3.08, and its mean reached a balance with all means (−1.28 a). All the means in Table 5 are negative, which demonstrated the leading tendency in detecting change time. The temporal errors are mostly due to the fact that our detection mainly focused on the change time of the impervious–relevant land disturbances caused by vegetation loss instead of the actual appearance of a built-up areas.

The evaluations are also displayed in a violin plot [65] (Figure 9) to illustrate the distribution of detection time deviation. The medians of the outputs are all −1, except for B7′s output, which is 0. The time deviations of all outputs are symmetrically distributed along both sides of the median. Valid value ranges for the outputs detected by NBR, NDVI, NDMI, B3, B7, and the integrated results are [−5, 3], [−5, 3], [−5, 3], [−5, 3], [−4, 4], and [−5, 3]. The separate sums of the corresponding valid values are214, 213, 206, 201, 206, and 214. The inputs’ differences in detection are shown based on the frequency of the data near the median although the number of valid values is close.

(2) Spatial accuracy

Spatial accuracy assessments are shown in Figure 10. NBR’s output exhibited the best spatial accuracy with an OA of 84.34%, a precision of 85.43%, and a recall of 96.37%. B3′s output had the lowest precision (83.50%). B7′s output possessed the lowest recall and OA (91.47% and 80.68%, respectively). Compared to NBR’s accuracy, the accuracy of the integrated result increased slightly. All spatial accuracy assessments are visualized in Figure 11. Matching the spatial reference in most regions (the shallow green and purple pixels), all outputs ideally detected land disturbances, especially for NBR and the integrated result. The integrated result effectively removed some isolated patches that obviously differed from the spatial reference (the shallow blue and dark blue pixels).

4.3. Expansion of Impervious Surface Dynamics in Nanjing

The spatial-temporal expansion of impervious surface dynamics in Nanjing was analyzed using the integrated map (Figure 8). Impervious surfaces in Nanjing has increased from 189.23 Km² in 1988 to 994.52 km² in 2018 (a more than five-fold increase). Figure 12 indicates the expansion differences in eleven districts. The expansion intensities of impervious surfaces in old towns are lower because of their limited space. These districts include Gulou, Qinhuai, and Xuanwu, which only increased approximately 1.37, 1.72 and 1.89 times compared to 1988. On the other hand, the expansion in Jiangning, Pukou, and Liuhe was much greater. Especially, the impervious surfaces in Jiangning increased from 23.39 Km² in 1988 to almost 235.37 km² in 2018, followed by Pukou and Liuhe (which increased by 124.70 Km² and 123.99 Km², respectively). The expansion patterns of impervious surfaces are also shown in Figure 8. City cores were also expanded from old towns to suburban areas. As a result, more suburban areas have served as the main functional zones of the city. Old towns, such as Gulou, Qinhuai, and Xuanwu, are more concentrated along the Yangtze River and Zhongshan Mountain, which are still the cultural and economic cores of the city. The parts of Jianye near Gulou, the regions of Pukou along the Yangtze River, the northern Jiangning, most areas of Yuhuatai, and the regions of Qixia next to Zhongshan Mountain became areas of further expansion. Since Nanjing is divided by the Yangtze River but connected by each district core, the final distribution pattern of impervious surfaces in Nanjing was shaped into a cross, from southwest to northeast and from north to south.

5. Discussions

5.1. Comparisons with Other Open-Access Products

Product comparison has become a common practice due to the lack of sufficient time-series references [9,14,21]. We compared our final results with those of other open-access products (Table A1) to discuss the differences in area estimates and mapping details (shown in Figure 13, Figure 14, Figure 15 and Figure 16). It is worth noting that none of these products were considered as the ground truth. This comparison aims to provide insightful improvements to our future work.

Figure 13 indicates the area estimates of the impervious surfaces from all products. Unlike NUACI-based mapping, the variation trend of the statistical yearbook records is similar to our study. GHSL shows the minimum gap (compared to this study) for 2000, with areas of 480.07 Km² and 450.66 Km², respectively. Compared to our area estimate in 2015 (916.44 Km²), FROM-GLC version2 and NUACI-based mapping in 2015 were found to have similar estimates to this study, with 1045.72 Km² and 1055.46 Km², respectively. Further, NUACI-based mapping has the highest estimates in 1990, 2000, 2005, and 2010. On the other hand, the FROM-GLC series in 2010 have relatively lower estimates.

Mapping details were amplified in four typical regions (Figure 14, Figure 15 and Figure 16). All regions are rectangles (with 14 × 13 Km), which represent the city center, suburb, countryside, and large infrastructure. Generally, all maps show similar spatial patterns. Figure 14 shows that FROM-GLC10 is more delicate than our study for 2017. FROM-GLC10 keeps more entire shapes of roads and built-up areas because of its spatial resolution with 10 m. As illustrated in Figure 15, all products in 2015 show a similar pattern in city cores. FROM-GLC version2 and GHSL exhibit a similar imperviousness distribution to this study for the city core and suburb but show small differences for rural and large infrastructure. NUACI-based mapping seems to lose patches in the suburb, countryside, and large infrastructure. In Figure 16, compared with the results of this study for 2010, FROM-GLC series and NUACI-based mapping show obvious differences in the distribution of impervious surfaces.

We admit that gaps caused by many factors indeed exist among these products, for instance, gaps in spatial scales, mapping methods, and data processing. More convincing reference products are still expected to evaluate mapping performance. The application of our scheme to Nanjing generated a dataset featuring the impervious–relevant land disturbance dynamics from 1988 to 2018, at yearly intervals. This dataset provides detailed information about the expansion of impervious surfaces in Nanjing both temporally and spatially. Thus, this dataset provides more change details at yearly intervals to fill the gaps that other products with long temporal intervals have missed. Over 30 years of detected data made our dataset helpful in understanding the characteristics and mechanisms of impervious surface expansion in Nanjing, and even the characteristics of the Yangtze River Delta urban agglomeration.

5.2. Parameter Sensitivity to the Output

Taking NBR as the example input, the sensitivity of the main parameters was discussed further. Giving the variation range for each chosen parameter to be tested, the value of each parameter was increased by one step length each time, until the value equaled the maximum. Other parameters were set as defaults when being evaluated for the chosen parameter, which is similar to the parameter configuration in Section 3.2.2 (2). All outputs derived with the increasing values were evaluated using temporal samples.

Figure 17a–h displays the temporal accuracy assessments of maxSegments, spikeThreshold, vertexCountOvershoot, preventOneYearRecovery, recoveryThreshold, pvalThreshold, bestModelProportion, and minObservationsNeeded, respectively. Figure 17d,h indicates that variations in preventOneYearRecovery and minObservationsNeeded barely affected the outputs with unchanged detection accuracies ($P=$ 85.46% and 85.46%, namely). Figure 17f shows that the alterations in pvalThreshold barely led to changes in detection accuracies and the median of ${\Delta }_Y$ ($P=$ 85.46% and median = −1) but slightly decreased the average (from −2.31 to −2.26) and the SD of ${\Delta }_Y$ (from 4.60 to 4.52). For bestModelProportion (Figure 17 g), $P$ abruptly decreased from 85.46% to 85.02% after an increased value greater than 0.90. Its SD experienced both a sharp decrease and a slight increase, but its average experienced an opposite variation. Obviously, variations in the other four parameters’ values sensitively affected the detection accuracies of the NBR runs. As shown in Figure 17b,c, both evaluations of spikeThreshold and vertexCountOvershoot show a peak-like variation, which may be because saturations exist in their increasing parameter values. Figure 17a,e illustrates that the values of maxSegments and recoveryThreshold are crucial to yield outputs with higher accuracies (the $P$ values are 88.11% and 88.99%).

Even if the parameters in LT-GEE were modified, our results still share some points with the LT-IDL’s parameter evaluation discussed by Kennedy et al. [31]. Success at capturing the most disturbances can be achieved by increasing maxSegments, and prohibition on over fast recovery could be enforced by setting recoveryThreshold to more stringent values [31]. Moreover, reasonable values of spikeThreshold and vertexCountOvershoot should be set to avoid saturation leading to lower detection performance. More prudent choices for the values of maxSegments and recoveryThreshold could significantly enhance the outputs.

5.3. The Proposed Scheme and Future Work

In this study, the proposed scheme extracted a developed region to implement long-term change detection in spatial constraints. To create the developed region, the unlimited ability of Landsat was used to employ the best-available-pixel (BAP) compositing method to reconstruct seasonal composites. In this way, any pixel with the best quality can be selected even if it belongs to the scene covered by heavy clouds. Moreover, extracting image features for classification using internal seasonal cloud-free composites can significantly improve the classification accuracy [9,66,67]. Our scheme also demonstrates that LT can be innovatively employed in the long-term monitoring of the impervious–relevant land disturbances associated with vegetation loss, although LT was adopted to detect changes and trends in forest disturbances and recovery [31,41,68]. On the other hand, our scheme is capable of capturing more subtle details because the dense Landsat archives allow us to find out where and when changes occur. Further, both abrupt and gradual changes can be derived by using straight line segments from historical spectral trajectories at the pixel level [23]. Moreover, the heterogeneity of the pixels in the same position caused by yearly classification results [11,28,29] can be avoided.

Our scheme facilitates the long-term monitoring of impervious–relevant land disturbances associated with vegetation loss at the pixel level. On the one hand, by combing the spatial constraints with a temporally-dense detection method, our scheme is capable of improving the temporal–spatial consistency of detection results. Thus, our scheme can efficiently trace other land cover dynamics with a tendency for continuous change (such as forest/grassland degradation). On the other hand, we focused on impervious–relevant land disturbances to reflect the actual expansion of built-up areas, aided by the understanding of the land conversion patterns in our study area. Therefore, our scheme also has a positive impact on strategies to capture the expansion of built-up areas in other tasks.

However, some issues of application still need to be considered pertinently in our future works. First, classification errors and images quality may (more or less) affect the stability of the developed region [9]. Better algorithms or integrated classifiers could be feasible to improve these flaws. Second, considering that all Landsat-7 ETM+ images were adopted in LT-GEE, we still need to explore how the stripe flaws caused by the failure of the Scan Line Corrector affect the output [69,70]. Third, the reference data for evaluation were collected using visual interpretation from Google Earth. Geometric correction differences between the VHR images and the Landsat images inevitably exist and may affect the evaluation results. Thus, an optimization of evaluation approaches is desirable. Finally, reversible and repeated changes between impervious surfaces and other land covers also wait to be identified.

6. Conclusions

To obtain a better understanding of the regional expansion of impervious surfaces at a finer temporal–spatial scale, we developed a new scheme to rapidly implement long-term monitoring of impervious–relevant land disturbances at yearly intervals. We first extracted the developed region using a classification-based approach. Then, we employed LT-GEE and high temporal frequency Landsat archives to detect impervious–relevant land disturbances with an annual bias. Finally, we derived the actual disturbances within the developed region and quantitatively evaluated the temporal–spatial accuracy of the results. Our study was conducted in Nanjing, which is a typical region that has experienced rapid urbanization in the Yangtze River Delta Urban Agglomeration in China.

The main conclusions are summarized as follows. First, based on the temporal–spatial accuracy assessments, our scheme has proven to be an effective paradigm with ideal accuracy. Second, monitoring of the impervious surface-relevant land disturbances provided a positive sight to capture the expansion of the built-up areas. Third, the detection performance and efficiency in heterogeneous backgrounds can be improved by combining spatial constraints with temporally-dense detection methods. Fourth, the scheme based on the GEE platform could save a significant amount of time in data preprocessing, algorithm implementation, and result optimization. Finally, our scheme performed well in Nanjing, which also has the potential to trace other land cover dynamics in an efficient way.

Our scheme can generate a geospatial dataset that includes impervious–relevant land disturbance dynamics during a given period, at yearly intervals. These datasets could help policymakers engage in macroscopic urban planning. For instance, these datasets are able to reasonably allocate lands for urban expansion, plan the main functional zones of a city with relevant infrastructures, and seek a balance between fast urban area expansion and natural resources. These datasets also make contributions to exploring urban ecosystem issues for sustainable urban development, such as understanding the characteristics of regional impervious surface expansions.