Identification and temporally-spatial quantification of geomorphic relevant changes by construction projects in loess landscapes: case study Lanzhou City, NW China

ABSTRACT

Sustainable development of urban areas demands, among others, a holistic approach for identification and monitoring of environmental changes caused by human activities. In this study, the human-made impact on the geo-environment is studied through the application of the change detection study based on Principal Component Analysis on Landsat imageries and comparison of the digital elevation models taken at different times. The analysis is set up to identify and quantify the anthropogenic geomorphological changes in the loess landscape in Lanzhou city, Northwest China. Since 2002 Lanzhou has been undergoing a rapid economic development associated with construction boom. Due to limited flat building ground in the narrowed Yellow River Valley and subsequent expansion into the surrounding loess mountains, massive earthworks are conducted for reclamation of the suitable building ground. The results of the change detection analysis show that approximately 10% of the semi-natural study area corresponding to 35 km² has been reshaped by leveling and terracing since 1994. In particular, the geomorphology was significantly changed in these and adjacent areas. For the single developing area Taipingyang, a moving volume of up to 57 million m³ was roughly estimated.

KEYWORDS:

• Change detection
• principal component analysis
• Chinese loess plateau
• economic development zone

1. Introduction

Human history records many cases of rising and fall of civilizations (e.g. Butzler, 2012; Ehrlich & Ehrlich, 2013). Among the complex socio-economic causes for the collapse, many of them also apparently suffered from unwise exploitation of environmental resources (e.g. McNeil, Burney, & Burney, 2010).

In modern times, since the industrial revolution, human activity has become an increasingly important global factor in ecological changes. In the environmental sciences, the impact of human activity is under discussion since March’s “Man and Nature” (Marsh, 1865). The term of “man as a geological agent” was introduced by Sherlock and Woodward (1922). Recent studies, e.g. Wilkinson (2005), show that compared to many natural processes acting on the earth surface, the human’s activity also becomes the most critical geomorphic agent. Price, Ford, Cooper, and Neal (2016) emphasized that the annual shift of sediments by human activities exceeding that transported by rivers almost by the factor of three. Therefore, in terms of sustainable development, it seems to be a logical step to identify, monitor, and evaluate the consequences of the environmental changes caused by human activity.

Today, with the development of the satellite remote sensing, imageries from various satellite missions with repetition rates ranging between several days up to a couple of weeks are available. Utilizing these advanced data, monitoring of the environment has become an increasingly standardized process. With a focus on multispectral satellite data and open access policy to the products, the most suitable data are provided by missions such as Moderate-resolution Imaging Spectroradiometer (MODIS), Landsat and Sentinel-2, with repetition rates of 16 days for MODIS and Landsat and 10 days for Sentinel-2, respectively (Heute et al., 2002; Li & Roy, 2017). Table 1 introduces some technical aspects of the mentioned remote sensing satellites.

Table 1. Technical data of most prominent remote sensing missions.

	MODIS	Landsat	Sentinel-2
Continuity since	1999	1972	2015
Swath wide [km]	2330	185	290
Spatial resolution	250 m (bands 1–2) 500 m (bands 3–7) 1000 m (bands 8–36)	30 m (pan 15 m)	10 m (bands 2–4, 8) 20 m (bands 5–7, 11,12) 60 m (bands1, 9, 10)
Band numbers	36	7 TM, 8 ETM+, 11 LDCM	13
Repeat cycle [d]	16	16	10 (5)*

TM – Thematic mapper, ETM+ – Enhanced thematic mapper, LDCM – Landsat Data Continuity Mission; *Bisection of the repeat cycle is facilitated by acting of two Sentinel-2.

Also, different methods have been developed in the past to identify and map environmental changes based on remote sensing imagery. Comprehensive overviews of the existing methodologies are given, among others, by Lu, Mausel, Brondízio, and Moran (2004) and Ban and Yousif (2016). Methodological approaches range from the simple algebraic methods such as image differencing, rationing or index calculation, to more or less complex transformation algorithms such as Change Vector Analysis (e.g. Dewi, Bijker, & Stein, 2017; Xiaolu & Bo, 2011) or Principal Component Analysis (PCA) (e.g. Deng, Wang, Deng, & Qi, 2008; Gong, 1993; Munyati, 2004). Also, unsupervised and supervised classification algorithms provide valuable options for change detection.

The current study focuses on Lanzhou, the capital of the Province Gansu. Due to its geographically and strategically important location (“New Silk Road – Node”), the surroundings of Lanzhou are subject to rapid social and economic development enforced by the establishment of a so-called Economic-Technological Development Zone (ETDZ). ETDZs are selected areas dedicated to push the development of specific regions and attract foreign investment (Liu & Wu, 2011). Such ETDZs were established countrywide by the Chinese national-level programs started in 1978. Since 1984 the number of ETDZs has been continuously growing. According to the report of the Ministry of Commerce, P.R.C. (Ministry of Commerce, 2016), recently, there are about 54 national-level ETDZs. The development activities in the Lanzhou area comprehensively started in the year 2002. Since then, new urban and industrial facilities have been constructing on the rural outskirts of the city. For the preparation of the suitable building ground, large scale bulk earthworks are conducted before infrastructure construction. The extent of the activities at which the natural thick loess deposits shaping a hilly landscape in the northern Lanzhou area has been leveled, and terraced gave them the name of the largest “mountain-moving project” in Chinese history. In total, about 700 mountainous peaks have been leveled in the entire Lanzhou area to create more than 250 square kilometers of flat land (Li, Qian, & Wu, 2014). Such tremendous changes are supposed to have a significant impact on the natural environmental conditions controlling groundwater supply or slope stability.

As Lanzhou area suffers from recurring hazards such as earthquakes, mass movements, and debris flows, sustainable development of new urban areas is only possible by implementing appropriate assessment of the hazards and consequent risks in the early stage of the spatial planning process. At regional scales, the landslide susceptibility representing the spatial distribution of the landslide occurrence is usually assessed to approximate the landslide hazard. Recently popular are data-driven (statistical) approaches using observations of the landslides from the past to establish a predictive model for the future. The application of statistically-based approaches requires several assumptions to be made concerning the existing environmental conditions. One of the critical points is the assumption of the stationarity of the environmental factors. Given the rapid development of the ETDZ of Lanzhou, stationary assumptions appear missconductive (Torizin et al., 2018). Instead, there is a substantial demand for detection and quantification of the environmental changes to be considered in the landslide susceptibility assessment (Tian et al., 2017; Torizin et al., 2018).

The presented paper introduces a change detection study based on free available satellite remote sensing imagery and a pair of digital elevation models (DEM) aimed to detect and spatially and temporally quantify the large scale bulk earthworks in Lanzhou area.

2. Study area

The study area is located on the western margin of the Loess Plateau, one of the largest loess deposits of the world, covering an area of about 600,000 km² (Wang, 2015). The thickness of the loess deposits in the surroundings of Lanzhou exceeds over vast areas 10 m and reaches up to 300 m in the northern part of Lanzhou (Derbyshire, Dijkstra, Unwin, & Wang, 2000; Dijkstra, 2000). Approximately 60% of the study area is covered by loess of different Pleistocene age (Figure 1).

Figure 1. Project area depicting the generalized lithostratigraphic settings based on Geological Map (Gansu Bureau of Geology and Mineral Resources (GBGMR), 1987) with underlying shaded relief produced in 2012.

The natural topography of the study area is characterized by undulating furrowed loess deposits wrapping a paleo-landscape influenced by the river incision of the Yellow River and tectonic activity induced by the ongoing Himalayan orogenesis. The elevation of the mountainous area is oscillating between 1500 and 2200 m above sea level (m a. s. l.).

As introduced above, due to its narrowed position in the Yellow River Valley, extensive building activity and expansion since 2002, Lanzhou suffers from limited suitable building ground. Therefore, the ambitious plan of developing the Lanzhou New Area (LNZ) can only be realized by expansion into the semi-natural hilly landscape. Being comparably soft sediment loess can be excavated, moved, and densified to leveled building ground by heavy machinery such as bulldozers.

3. Methodology and analysis

3.1. Principal component analysis in change detection

Specific materials at the surface are exhibiting measurable spectral differences, which might be of interest for the application of change detection. The distinguishing of freshly exposed loess material as a result of earthwork activities from other spectral classes such as vegetation of varying vitality or urban area is the primary scope.

In the current study, PCA was chosen among the methods mentioned in the introduction part as the potentially appropriate method given the local conditions of the study area. Change detection PCA became a standard methodology in remote sensing over the past decades (e.g. Byrne, Crapper, & Mayo, 1980; Gong, 1993). Nevertheless, it is worth it to be briefly introduced.

PCA is a technique applied to reduce dimensions in multidimensional data and to identify patterns highlighting their similarities and differences. PCA is an orthogonal transformation method to convert a set of possibly correlated variables into a linear combination of a set of pair-wise uncorrelated variables called principal components (PCs). Thereby, the covariance of the variables subdivides into two parts: the scale part and direction part. The so-called eigenvectors characterize the direction part forming an orthonormal basis, and thus, constitute the independence (no correlation among eigenvectors) of the PCs. The scalar values of these eigenvectors the so-called eigenvalues are considered as a magnitude describing the variance in a particular direction. The PCs are defined in such a way that the maximum variance of the variables is held in the first PCs and the amount of variance is decreasing with subsequent PCs. This means, for the set of correlated variables, the sufficient amount of information is contained in the first few PCs allowing a reduction of the dimensions of the dataset. The normalized eigenvalues of the eigenvectors provide the proportion of the overall variance explained by the PC, in some packages dealing with PCA also denoted as importance. The decision of how many PCs should be considered in the analysis is frequently made applying the Kaiser Criterion. This criterion compares the normalized eigenvalues (importance) and provides a cutoff threshold for the PCs with importance ≤ 1.

Applied for change detection on a multi-temporal multispectral image stack, the first PC is bearing the information that is common to all images in the stack (unchanged areas). The differences marking changed areas are contained in the following PCs. Thus, changes are identifiable by the interpretation of the loadings of the eigenvectors and the binning of the factor values of the PCs (Drury, 1993; Gong, 1993).

To support the decision on how much difference is contained in which PC, a contrast of the eigenvector loadings C can be introduced. The contrast C can be calculated as:

$C=abs\left({\sum }_{i=1}^n{c}_{{ }_A}^i-c_{{ }_B}^i\right),$

where cⁱ_A is the ith band of scene A and cⁱ_B is the ith band of scene B.

The interpretation of the contrast is straightforward. The larger the contrast, the higher are the differences between the two scenes contained in the specific PC. Moreover, one can focus on the contrast of the eigenvector components (bands) most sensitive to particular surface changes and use them for the selection of the most appropriate PC to identify the changed areas. For geological mapping, short wave infrared (SWIR) bands are known to provide high contrast values. Focusing on the SWIR contrasts in the PCs supports the decision on the selection of the PC that is best for identification of areas modified by bulk earthworks. Once the PC with the most considerable SWIR contrasts among the two scenes is found, the factor values of this PC – rescaled between 0 and 255 – can be used to determine the threshold to separate the changed from unchanged areas. Because the factor values are usually nearly normally distributed, a threshold T for the delineation of such areas can be found, e.g. considering standard deviation σ with a value smaller or greater than T (Gong, 1993).

In this study, the 2nd and the 98th percentile were proved as applicable thresholds (Figure 2). These percentiles are used as thresholds to separate approximately the extreme values in a normal distribution. The changed areas can be now interpreted in both directions, exposed areas due to earthwork activities or newly sealed areas, e.g. by greening or construction work. The decision which side of the distribution is to take to mark the exposure is obtainable from the loadings of the second time scene for the SWIR components. Negative loadings indicate exposed areas; positive loadings indicate sealing of the areas by construction or greening. Because sealing activities are post-earthwork activities, they are not in the focus of the current study.

Figure 2. Classification steps of changed areas, illustrated on PC2 from image stack 2013–10-11 to 2015–10-10: a) Selected greyscale coded PC containing the changes; b) Binning statistics of the factor values with 2nd and 98th percentiles representing changes; c) Mapping of the changed areas according to the binning statistics.

3.2. Tools and data

The change detection analysis was performed on Landsat imaginary in GRASS GIS (GRASS Development Team, 2017), providing the PCA as a standard processing tool. Postprocessing of the PCA results was made in QGIS (QGIS Development Team, 2017).

Landsat imagery from comparable annual time without cloud cover was used to minimize the influence caused by phenological or climate conditions. In total, six scenes from the Landsat 5 mission called as Thematic Mapper (TM), and two scenes from the Landsat 8 mission with the Operational Land Imager (OLI) were selected covering a period between 1994 and 2015 (Table 2). The Landsat imagery acquired from USGS Earth Explorer is atmospherically corrected and transformed to the surface reflectance (NASA, 2016). Each of the selected scenes contains six solar reflectance bands: blue, green, red, near-infrared (NIR), and the two SWIR. The obtained data stack enables a subdivision of the 21 years in seven intervals ranging from two to seven years.

Table 2. Landsat scenes included in the study.

Mission	Instrument	Scene id	Acquisition date
Landsat 5	TM	LT51310351994296BJC00	1994–10-23
Landsat 5	TM	LT51310352001299BJC00	2001–10-26
Landsat 5	TM	LT51310352004276BJC00	2004–10-18
Landsat 5	TM	LT51310352006281IKR00	2006–10-08
Landsat 5	TM	LT51310352009289IKR02	2009–10-16
Landsat 5	TM	LT51310352011279IKR02	2011–10-06
Landsat 8	OLI	LC81310352013284LGN00	2013–10-11
Landsat 8	OLI	LC81300352015283LGN00	2015–10-10

The reflectance data of the blue, green, red, NIR, and both SWIRs were included in the PCA for each image and calculated pairwise. Thus, in each PCA, twelve bands were scrutinized starting with the scenes from the years 1994 and 2001, followed by scene combination for the years 2001 and 2004 and so on.

3.3. Analysis workflow

The general workflow for change detection in this study can be described in the following steps (Figure 3):

1. Generation of the GRASS GIS project (import data, set region extent and resolution);

2. Perform PCA subsequently for Landsat scene pairs;

3. Dimension reduction to the significant PCs based on Kaiser-Criterion (e.g. Backhaus, Erichson, & Plinke et al., 2015);

4. Screening the contrast to reveal the PC containing the most significant differences among the loadings of the eigenvectors;

5. Estimation of the factor values thresholds indicating changes in selected PC layer based on the loadings of the SWIR bands;

6. Exclude the city area located in the Yellow River valley;

7. Visual adjustment of the chosen factor values with the respective images in QGIS;

8. Validation of the mapped changed areas

9. Perform area statistics per time interval.

Figure 3. Analysis flow chart.

For the area, descriptive statistics were calculated. The loadings of the input bands in the new PC were utilized for that purpose. The loadings indicate the contribution of the single input bands in the process of forming the PC. Based on the algebraic sign of contributions of the SWIR bands, it can be decided whether the low or the high factor values in the newly formed PC layer is representing earthwork areas.

4. Results

4.1. PCA results

Table 3 lists the PCA results for PCs with eigenvalue importance greater than one. In particular, it displays the list of PCs for the seven scene pairs with corresponding eigenvector loadings describing the contribution of each component to the eigenvector.

Table 3. Important PCs, eigenvalues, loadings of the eigenvector for corresponding six bands, importance, and contrast C; the utilized PCs are in bold.

PC	Eigenvalue	Vector components scene 1						Vector components scene 2						Importance [%]	C
		Blue	Green	Red	NIR	SWIR1	SWIR2	Blue	Green	Red	NIR	SWIR1	SWIR2
		1994–10-23						2001–10-26
PC1	3,455,444	0.05	0.11	0.16	0.25	0.41	0.39	0.07	0.13	0.19	0.29	0.47	0.46	89.0	0.24
PC2	160,658	0.35	0.38	0.38	0.32	0.17	0.25	0.15	0.11	0.05	−0.06	−0.42	−0.41	4.1	2.43
PC3	144,181	0.19	0.15	0.08	−0.08	−0.45	−0.37	0.36	0.42	0.45	0.23	0.02	0.15	3.7	2.11
PC4	59,702	−0.08	−0.09	−0.14	0.46	0.01	−0.29	−0.11	−0.10	−0.13	0.71	0.05	−0.34	1.5	0.21
		2001–10-26						2004–10-18
PC1	3,622,996	−0.06	−0.13	−0.18	−0.27	−0.46	−0.45	−0.09	−0.15	−0.20	−0.26	−0.39	−0.41	88.33	0.05
PC2	204,673	−0.22	−0.28	−0.31	−0.40	−0.24	−0.24	0.07	0.12	0.21	−0.03	0.40	0.53	4.99	2.98
PC3	127,301	0.34	0.37	0.37	−0.15	−0.29	−0.07	0.31	0.36	0.42	−0.26	−0.16	0.09	3.1	0.19
PC4	88,934	−0.06	−0.04	0.02	−0.37	0.27	0.45	−0.09	−0.14	−0.12	−0.73	−0.06	0.11	2.17	1.29
		2004–10-18						2006–10-08
PC1	3,515,889	−0.09	−0.16	−0.21	−0.26	−0.40	−0.42	−0.10	−0.17	−0.22	−0.26	−0.41	−0.44	89.58	0.05
PC2	145,954	0.10	0.10	0.17	−0.66	−0.07	0.20	0.10	0.12	0.20	−0.58	0.01	0.28	3.72	0.30
PC3	126,064	−0.23	−0.28	−0.29	−0.04	0.35	0.27	−0.32	−0.39	−0.44	−0.28	0.19	0.17	3.21	0.85
PC4	93,876	0.19	0.27	0.35	0.17	0.32	0.41	−0.05	−0.09	−0.12	−0.22	−0.42	−0.47	2.39	3.07
		2006–10-08						2009–10-16
PC1	2,859,574	0.11	0.19	0.24	0.29	0.45	0.48	0.06	0.11	0.14	0.26	0.38	0.35	84.6	0.46
PC2	232,268	0.01	−0.04	−0.09	0.00	−0.37	−0.46	0.28	0.33	0.41	0.26	0.25	0.40	6.9	2.88
PC3	132,511	0.28	0.33	0.42	−0.30	−0.09	0.15	0.21	0.22	0.28	−0.53	−0.25	0.07	3.9	0.79
PC4	102,582	0.19	0.26	0.27	0.63	−0.07	−0.22	0.07	0.09	0.07	0.27	−0.33	−0.41	3.0	1.30
		2009–10-16						2011–10-06
PC1	2,498,583	0.08	0.13	0.17	0.29	0.41	0.39	0.10	0.15	0.20	0.27	0.44	0.45	84.4	0.14
PC2	169,068	0.30	0.32	0.40	−0.33	−0.24	0.10	0.30	0.34	0.41	−0.23	−0.23	0.06	5.7	0.10
PC3	138,161	0.18	0.24	0.29	0.42	0.27	0.28	−0.05	−0.08	−0.14	0.14	−0.41	−0.54	4.7	2.76
PC4	99,603	0.01	0.03	0.00	0.32	−0.34	−0.48	0.12	0.22	0.26	0.64	0.01	−0.11	3.4	1.60
		2011–10-06						2013–10-11
PC1	2,767,680	0.09	0.13	0.18	0.25	0.41	0.42	0.12	0.16	0.22	0.29	0.43	0.41	82.7	0.15
PC2	255,650	0.12	0.19	0.25	0.21	0.40	0.46	−0.15	−0.22	−0.30	−0.26	−0.32	−0.37	7.6	3.25
PC3	156,001	0.30	0.36	0.44	−0.17	−0.25	0.02	0.29	0.36	0.40	−0.20	−0.30	−0.02	4.7	0.17
PC4	114,213	−0.04	−0.11	−0.10	−0.68	0.12	0.31	−0.02	−0.07	−0.03	−0.52	0.17	0.31	3.4	0.34
		2013–10-11						2015–10-10
PC1	3,846,825	0.10	0.14	0.19	0.25	0.37	0.35	0.15	0.21	0.27	0.31	0.44	0.44	87.1	0.42
PC2	224,850	0.23	0.31	0.37	0.32	0.30	0.39	−0.04	−0.09	−0.16	−0.11	−0.39	−0.42	5.1	3.13
PC3	165,268	0.21	0.28	0.31	−0.23	−0.38	−0.16	0.30	0.42	0.49	0.01	−0.20	−0.10	3.7	0.89
PC4	126,999	−0.04	−0.01	−0.07	0.49	−0.18	−0.32	−0.04	0.02	−0.03	0.73	−0.04	−0.28	2.9	0.49

The interpretation of the results is straightforward. The respective first PC holds almost 90% of the variance of the original data and represents unchanged areas in the scenes of different dates emphasized by the loadings of the eigenvector, which are similar among the components (bands) of both scenes. Small variations in atmospheric and light and shadow conditions can explain small differences of the loadings. Also, the contrast value shows that the loading differences of the first PC are neglectable. For the majority of scene pairs, the second PC (PC2) possessed the highest contrast and was found appropriate to be used for change identification. For one scene pair, the third PC (PC3), for another, the fourth PC (PC4) showed the highest contrast and were consequently applied for further analysis. As expected, the contribution of both SWIR bands to the contrast value is comparably high in all PCs, which has been found suitable for identification of changes in loess, highlighted in Table 3.

Figure 4 depicts the changed areas for corresponding scene pairs based on the 2nd, and the 98th percentile classification. In total, the estimated study area is about 416 km². However, the central area of Lanzhou City was excluded from the analysis as it was already entirely urbanized. Therefore, the total area considered for the change detection and used for calculation of the statistics is about 353 km².

Figure 4. Changed areas in time intervals and validation spots; depicted on the Sentinel-2 scene [11/01/2016] (Natural colors: RGB bands 4, 3, 2). The city area is excluded.

The results illuminate that in the last two decades, about 10% of the semi-natural loess area has been affected by bulk earthworks and construction activities (Figure 5). These activities exhibit surficial removal of material at different scales. For large scale projects, these activities are associated with the massive reshaping of the landscape morphology, including leveling of hills and filling up of preexisting valleys.

Figure 5. The proportion of the leveled and bulldozed area in the respective periods.

4.2. Validation of PCA classification

Checking the accuracy of the results is based on the following validation sources:

• Examination in the field;

• Visual interpretation of the high-resolution imagery in Google Earth.

The spatial pattern of changed areas in Figure 4 shows that the majority of the leveling and reshaping activities have been conducted in the north of Lanzhou city.

The indicated geomorphological changes were validated in a field campaign on three spots, the Taipingyang-, Biguiyuan- and Jiuzhou development zone (Figure 4), which are all located to the north of Lanzhou City. All three locations, which were a grassland or bare ground with a semi-natural loess morphology before, are now under preparation for construction or already occupied by housing communities and commercial business districts. The new land created by the leveling and reshaping activities was compacted, and the artificial slopes were engineered.

One zone, the Taipingyang development zone in the north-eastern part of Lanzhou City, is depicted in Figure 6(a,b). In this area, the loess mountains were leveled in the period from 2012 to 2013. With an elevation model of 10 m resolution produced before the year 2013, the minimum and maximum altitude values in that area were estimated (NASG, 2012). At the beginning of 2012, the altitude ranged between 1564 m and 1840 m. In an area of about 4.4 km², the morphology was changed from a hilly loess landscape with some fluvial channels to a peneplain by leveling of more than 200 m in height. Massive construction activities have not been started yet in this area.

Figure 6. Taipingyang development zone, July 2017, a) leveling and b) terracing work to generate construction land.

Comparable changes are observable in the Biguiyuan development zone. However, here the flat area is already consumed for erecting of high-rise buildings. The slopes on the boundary of the project area have been engineered by step-slope excavation (Figure 7(a)).

Figure 7. a) Housing community in the Biguiyuan development zone, 27–04-2016, b) Commercial building district and housing community in the Jiuzhou development zone, 26–04-2016.

The Jiuzhou development zone, former an unused land located northeast of Lanzhou city, has now become a popular and busy region filled with housing communities and commercial buildings (Figure 7(b)).

For visual interpretation, we used the high-resolution imageries provided by Google Earth. The number of high-resolution images covering specific periods varies in Google Earth depending on the area of interest. For Lanzhou more than 60 images from different dates covering the period from 2001 until today are available allowing spot-wise validation of regions of interest for particular time intervals. Figure 8 shows an example of visual interpretation. In the eastern part of the today’s completely developed Biguiyuan project area changes were revealed between Landsat scenes from Oct 2001 and Oct 2004 by the PC 2 (Figure 8(a)). The high-resolution images available from November 2002 and January 2004 indicate that the leveling activities started in 2002 and were conducted mainly in 2003–2004 (Figure 8(b–c)). The shapes of automatically detected areas and visually recognizable changes are matching well, considering a different resolution of the compared imagery.

Figure 8. Example of a visual check with Google Earth, showing leveling activities in Biguiyuan development zone in 2002–2004: a) Factor values of the change detection between the scenes Oct 2001 – Oct 2004. b) The beginning phase of earthworks in 2002 [11/25/2002]; c) Shape of the development area in November 2004 [11/04/2004]. Differences in the brownish hues indicate the leveling area in the center of the picture.

4.3. Estimation of moved volumes

At first glance, a comparison of two DEMs of different dates seems to be a straightforward method to estimate geomorphological changes directly from the DEM comparison. So, high-resolution DEMs based on airborne LiDAR are recently used to visualize and to quantify geomorphological changes in dynamic regions such as coastal erosion processes, e.g. after storm events (e.g. Burvingt, Masselink, Russel, & Scott, 2016). Unfortunately, this proves to be not that easy for arbitrary areas of the world in which no consequent monitoring program and subsequently no repetitive airborne laser scanning (LiDAR) measurements were conducted during the last years. Although different satellite missions provide DEMs in a ground resolution of 10 m covering a large part of the world, some issues remain in a direct comparison of these data. For the project area, we acquired two DEMs from different dates with a ground resolution of 10 m. The first DEM is a product from the ZY-3 satellite mission (optical sensor) generated from a stack of scenes of 2012. The second DEM is a custom-tailored space born radar-based product from the TerraSAR-X satellite mission generated from the data acquired in the period from July to August 2016 (Astrium GEO-Information Services, 2011).

The first comparison of the two DEMs revealed that the differences between both DEMs are not interpretable without additional information due to different sensors and not comparable data processing approaches. Uncertainties in both DEMs cause general differences also in areas not affected by bulk earthwork activities. Especially in the southern part of the study area, high relief regions exhibit significant deviations between both DEMs that cannot be lead back to a specific change of the morphology. Because, the radar-based DEMs are generally known for issues in areas of high relief, the general interpretation of all computed differences between both DEMs as anthropogenic changes would be misleading. However, for large construction sites detected in the PCA as coherent areas and verified by field checks and visual interpretation of high-resolution imageries, the differences can be interpreted. Thus, based on the DEM data, the volumes of the displaced material during the development of the Taipingyang zone was roughly estimated (Figures 6 and 9). The results of the PCA change detection were used to generate a mask layer in which the differences between the two DEMs could be interpreted meaningfully. The obtained difference pattern inside the masked area reveals the former ridge and valley structure leveled during the reclamation earthwork (Figure 9(c)). Based on the obtained difference pattern, a rough balance for cut and fill activities was established. The results show that 1.7 km² of the 4.6 km² area was leveled by cutting and 1.9 km² by filling. Considering elevation uncertainties of both DEMs of about 10 m, the morphology in the Taipingyang development zone was changed by cutting hills up to 110 m and filling valleys up to 100 m. In this area, the Precambrian granite and the Cretaceous sediments former outcropped in the valleys were covered by remolded loess of several tens of meters. A rough volume of the displaced material was estimated at 57 million m³.

Figure 9. Elevation difference classes on the validation spot Taipingyang a) elevation in 2012, b) elevation in 2016, c) cutting and filling difference classes depicted on a Sentinel-2 scene [11/01/2016] (Natural colors: RGB bands 4, 3, 2).

5. Discussion and conclusions

The presented change detection applied with a one-sided interpretation of the PCA-results is suitable for identification areas in which the sediments were exposed and moved. In particular, in the loess dominated landscape, the bright reflectance of the freshly exposed loess in the SWIR bands contributes to the PCs indicating the changes between two dates. The contrasts of these bands contribute more than 50% of the total contrast in the two scenes.

Based on the change detection from multispectral satellite imageries, no information can be derived about the magnitude or type of earthwork activity that caused the detected change. However, the more coherent the changed areas are, the more likely is it that the earthworks were conducted as cut and fill activities. Additionally, change detection must be combined with other methods such as ground truth and a visual interpretation of high-resolution imageries.

Great care is needed comparing DEMs from different origins. The difference between the elevations of optical-based and radar-based DEMs cannot be interpreted successfully without additional information due to uncertainties of the DEMs caused by issues of data acquisition and processing. The introduced application of change detection can support the estimation of material volumes displaced by earthwork activity providing verified spots for which the differences of the DEMs can be interpreted more meaningfully. In this context, it should be noticed that, of course, elevation models based on LIDAR campaigns in fixed time intervals would provide more reliable and accurate results for volume estimation at those spots. However, such activities are costly and must have been done before the changes occurred to have a precise reference surface. Compared to LIDAR, utilization of increasingly popular photogrammetric solutions in the field of Unmanned Aerial Vehicles (UAV) could reduce the costs, but currently, the covered areas remain to be comparably small and still cannot serve for assessments at regional scales.

Economic zones, corridors, or spots in the P. R. of China are subject to continuous development. In particular, the northwestern part of the country, including the central loess plateau is affected by the increasing urbanization and establishment of new economic spots. The proposed method applied to identify changes caused by earthworks can be used for other loess dominated landscapes.

Characterization of the spatiotemporal environmental changes can provide constraints for many further analyses such as statistical landslide susceptibility assessment. Because statistical methods are not designed to reveal causal relations, but the association of specific variables, the knowledge of the locations and the speed of temporal changes have to be incorporated. The findings of this study help in the context of statistical landslide susceptibility analysis to create a framework describing the human-induced evolution of the landscape.

The continuation of the change detection time series based on Sentinel-2 data, or based on a combination of Landsat and Sentinel-2 would increase the probability to have access to cloudless scenes also in shorter time intervals.

Data availability statement

Derived data supporting the findings of change detection study are available from the corresponding author on request. Restrictions are given for the compared elevation models due to contract.

Disclosure statement

No potential conflict of interest was reported by the authors.