A New Method to Predict Gully Head Erosion in the Loess Plateau of China Based on SBAS-InSAR

Abstract

Gully head erosion causes serious land degradation in semiarid regions. The existing studies on gully head erosion are mainly based on measuring the gully volume in small-scale catchments, which is a labor-intensive and time-consuming approach. Therefore, it is necessary to explore an accurate method quantitatively over large areas and long periods. The objective of this study was to develop a model to assess gully head erosion in the Loess Plateau of China using a method based on the SBAS-InSAR technique. The gully heads were extracted from the digital elevation model and validated by field investigation and aerial images. The surface deformation was estimated with SBAS-InSAR and 22 descending ALOS PALSAR datasets from 2007 to 2011. A gully head erosion model was developed; this model can incorporate terrain factors and soil types, as well as provides erosion rate predictions consistent with the SBAS-InSAR measurements (R² = 0.889). The results show that gully head erosion significantly depends on the slope angle above the gully head, slope length, topographic wetness index, and catchment area. The relationship between these factors and the gully head erosion rate is a power function, and the average rate of gully head erosion is 7.5 m³/m²/year, indicating the high erosional vulnerability of the area. The accuracy of the model can be further improved by considering other factors, such as the stream power factor, curvature, and slope aspect. This study indicates that the erosion rate of gully heads is almost unaffected by soil type in the research area. An advantage of this model is that the gully head area and surface deformation can be easily extracted and measured from satellite images, which is effective for assessing gully head erosion at a large scale in combination with SBAS-InSAR results and terrain attributes.

1. Introduction

Gully erosion is one of the most hazardous forms of soil erosion [1,2]. It is a severe land degradation problem in almost all climatic and geographical areas, especially in the Loess Plateau [3,4]. The DongZhi Loess Tableland is the largest tableland of the Loess Plateau, and has an important and strategic position in China [5]. Gully erosion-related ecological and resource crises in this region severely impede its social and economic development [6]. Gully erosion is considered a major challenge for land resource degradation and is becoming increasingly widespread [7]. In recent decades, gully erosion has been recognized as an important environmental issue [8].

A gully system consists of the gully line, gully head, and shoulder line. In a gully system, the gully head is one of the most active parts: It controls the conflux of the watershed, affects the development of the loess landforms, and causes the loss of land resources [8,9,10]. A gully channel extends, widens, and deepens from the gully head [11]. Gully expansion mainly occurs by gully head erosion [12]. Therefore, determining the location of a gully head and exploring a model to predict the gully head erosion are significant.

There are many methods of determining the location of a gully head. Field surveying was the earliest research method used [12]. Currently, the main gully head datasets in the region are extracted from aerial images and digital elevation models (DEMs) [13]. Visual interpolation and pixel-based image analysis are reliable techniques that apply spectral and textural heterogeneity to extract gully systems. However, these approaches are sensitive to the land cover type and atmospheric conditions. Object-oriented analysis (OOA) is an effective method used to monitor gully erosion over large areas [14]. However, the results of OOA methods rarely reveal basic geomorphological properties (e.g., slope angle, catchment area, plan and profile curvature, and slope aspect) or structural features of the terrain (e.g., gully channel line and gully head location). These features are important for understanding the development process of a gully [14]. Some DEM-based methods can derive gully terrain information, but cannot accurately extract the gully head location [13]. Under these conditions, DEM-based hydrologic methods are used to extract the gully head. The gully head is accurately determined by the point of intersection of the gully shoulder line (boundary between positive and negative (P-N) terrains) and gully network line (channel network) [13].

Determining the methods for monitoring and analyzing gully head erosion has been a research priority [12]. Traditional methods monitoring gully head erosion include field surveys, positioning observations, and GPS methods. Traditional methods are labor-intensive, time-consuming, and applicable to only small catchments [15,16]. Moreover, the precision of monitoring may not be high. Gully head erosion has also been determined using a DEM, and due to the inherent accuracy limitations, this method has not been commonly used [17]. A 3D scanner has successfully obtained a high-resolution DEM, which can be used to assess surface landforms and gully erosion. However, it is not suitable for assessing deformation and erosion at a considerable regional scale. Compared to the abovementioned gully head erosion evaluation technologies, satellite (e.g., synthetic aperture radar interferometry (InSAR)) and conventional aerial remote sensing monitoring methods (e.g., unmanned aerial vehicles (UAVs)) have significant advantages, including the ability to quickly monitor surface deformation, an insensitivity from climatic conditions, and the production of large-scale and high-resolution products [16,18]. In general, UAVs are more advantageous for short-term gully head erosion monitoring, but are also more sensitive to the environment. In addition, UAVs cannot be used to monitor continuous and dynamic deformation of gully surfaces due to their limited flight range. Gully head surface deformation contributes to gully head erosion and is thus considered a measure of gully head erosion. To predict the trend of a gully head, it is necessary to consider the gully head development over time, and one direct method used to understand this development is to measure the change in gully surface deformation. InSAR is a high-resolution technique mainly used for monitoring surface deformation over long periods and large areas [16,19,20]. Furthermore, multitemporal InSAR is suitable for analyzing time series images [21]. The small-baseline-subset InSAR (SBAS-InSAR) method [22] and the permanent-scatter InSAR (PS-InSAR) method [23,24] have been widely used in surface deformation and have the measurement accuracy of millimeters. Many studies indicate that the SBAS-InSAR method is more suitable than the other InSAR methods for addressing the low coherence of vegetation and gully areas [21]. The SBAS-InSAR method has been used in various landslide, urban subsidence, and volcano studies, for location detection and deformation monitoring [25,26]. However, it is rarely used to monitor gully deformation [21,26,27]. The SBAS-InSAR method can compensate for the lack of historical disaster data [21,27] and the gap in the evaluation of gullies at the regional scale [16]. Therefore, exploring the application of the SBAS-InSAR method in gully deformation may expand the spatial and temporal scale of gully head deformation monitoring, and promote the quick assessment of gully head deformation at the regional scale.

Samani [28] concluded that geomorphic threshold factors (e.g., catchment area and slope angle) were the major factors that affected gully erosion, and that soil attributes affected the gully initiation thresholds in arid and semiarid regions. Thompson [29] utilized the most important factors that affected gully erosion (e.g., catchment area and soil type) to statistically analyze data and obtained a logarithmic model for the United States. Based on previous studies, it is necessary to investigate the gully head erosion threshold and predict the gully erosion rate by examining the relationship among terrain attributes within the Loess Plateau of China.

The study area is one of the main gully head erosion areas, where gully head erosion has resulted in land resource degradation and infrastructure damage. In this study, we explore a new method to measure gully head erosion over large areas and long periods. The objective of this study is to measure and develop a model to assess gully head erosion on the Chinese Loess Plateau using a method based on SBAS-InSAR. The gully head locations are extracted from DEM data, the surface deformation of the gully head is monitored with ALOS PALSAR data, and the correctness of the results is verified.

2. Research Area

The DongZhi Loess Tableland, in the eastern reaches of the Malian River basin in northern China (106°14′–108°42′E, 34°50′–37°19′N), is covered with strongly erodible loess layers, which are more than 100 m thick on average (Figure 1a). This study is focused on Qingcheng County in the northeastern DongZhi Loess Tableland, eastern Malian River basin (Figure 1a). The geographical coordinates are 117°45′–117°52′E, 35°81′–36°18′N, and the research area is 29.3 km² (Figure 1a). The most common soil type in the research area is the thick wind-blown loess sediments deposited in central Asia during the Quaternary [30]. The lithological profile is divided into four units, and shown in Figure 1b. From top to bottom, these units are as follows: 1. The upper Pleistocene aeolian Malan loess (Q3) with an average thickness of 8 m; 2. the Lishi loess (Q2), which is more than 50 m thick; 3. an indurated layer of the Pliocene and Miocene Lantian red clay deposit (N2), 3–30 m tick; and 4. the lower Cretaceous mudstone bedrock with sparse sandstone partings. Q3 predominantly covers the inter-gully areas, whereas Q2 and red earth are mainly observed in the side slopes and bottoms of the gullies and short-lived river channels. Gully erosion in this small watershed is severe. The field survey found that the soil derived from the unconsolidated areas at the bottom of the gullies, and the slope surface that also contained Q3, was transported by the gully head soil erosion (Figure 2).

The study area is in a semiarid continental climate zone, which experiences both monsoons and the typical climate of the Loess Plateau. The average annual temperature is 13.7 °C (Figure 3). According to the data recorded by the local meteorological bureau, the annual precipitation was 451–650 mm in 2007–2011, with an average of 55.1 cm, and an average of half of the precipitation occurred during the monsoon season from June to September. On the Loess Plateau, extreme rainfall contributes to the expansion of gully erosion [31].

The terrain in the research area is sloping, rising from the northeast to the southwest. The elevation range is 1054–1384 m, the average altitude is 1250 m, and it is, densely dissected by gullies. The land use type is mainly cropland and grassland in the research area, which accounts for more than 77.6% of the watershed [32].

3. Materials and Methods

The methodology (Figure 4) in this study includes four steps: 1. Extracting the gully head location; 2. selecting factors that affect the spatial variation in gully head erosion; 3. calculating the deformation of the gully heads using SBAS-InSAR techniques; and 4. validating the accuracy of the gully head deformation model.

3.1. Terrain Attributes from a DEM

In this study, the DEM was created by digitizing the contour lines using topographic maps at a 1:10,000 scale. The contour interval of this topographic map is 10 m. The original resolution of the DEM is 10 m. The contour lines were interpolated using the Thin Plate Spline method (TIN) to generate the DEM with 5 × 5 m² cells. All the terrain attributes in this study are extracted from the DEM: The slope angle (Sl, °); catchment area (Ca, m²); topographical wetness index (TWI, m radian); stream power index (SPI, m radian⁻¹); slope-length factor (LS, m); plan curvature (Curv_pl, m m²); profile curvature (Curv_pr, m m²); and slope aspect (Sa). The terrain analysis extension of ArcGIS 10.2 was utilized in this study to estimate those terrain factors (Figure 5).

Sl is strongly correlated with the soil erosion and accumulation flow surface in an area [33]. For example, gentle slope areas due to overland flow accumulation are exposed to gully initiation. The range of slope in the research area is from 0 to 88° (Figure 5a). However, the slope of approach channel above the gully heads is between 0 to 40°.

Ca refers to the cumulative amount of flow directly or indirectly flowing into the grid at the gully head [34]. The range of Ca is from 0 to 600 m² (Figure 5b).

TWI indicates the tendency of gully heads to receive overland flows and could provide information on gully formation and evolution [35]. This factor is calculated by Equation (1):

$$\mathrm{TWI}=\mathrm{In} \left(\mathsf{\alpha }/\tan \beta \right)$$

(1)

where $\mathsf{\alpha }$ is the catchment area and $\beta$ is the slope gradient (in degrees). The range of TWI in the research area is from 0 to 22 m radian (Figure 5c). The range of TWI in the gully heads of the research area is from 0–14 m radian.

SPI is a relatively effective indicator for estimating the erosive power of surface runoff and surface erosion [36]. This factor is calculated by Equation (2):

$$\mathrm{SPI}=\mathsf{\alpha }\ast \tan \beta$$

(2)

where $\mathsf{\alpha }$ is the catchment area and $\beta$ is the slope gradient in degrees. The range of SPI of the research area is divided into four intervals (Figure 5d).

LS mainly affects the accumulation and acceleration of flow in the gully head erosion [33]. This factor is calculated by Equation (3):

$$\mathrm{LS}={\left(\mathsf{\alpha }\ast \mathrm{Cell} \mathrm{size}/22.3\right)}^{0.4}\ast {(\sin \beta /0.0896)}^{1.3}$$

(3)

where $\mathsf{\alpha }$ is the catchment area and $\beta$ is the slope gradient in degrees. The range of TWI in the research area is from 0 to 50 m (Figure 5e).

Curv_pl and Curv_pr can be used to consider the effect of the local terrain’s morphometric on overland flow distribution, runoff, and consequently gully head erosion [35]. Curvpl is tangent to a contour and is defined as the curvature on the corresponding normal interface. Curv_pr can be regarded as the second derivative of the slope angle (Figure 5f,g).

Sa through influences of the weathering mechanism, soil moisture, vegetation type, and geomorphological processes can indirectly affect gully head erosion processes. The range of Sa of the research area is divided into nine intervals (Figure 5h).

3.2. Gully Head Detection with a DEM

The gully head is a small area below the plateau, and its morphological characteristics are relatively stable during a certain stage of landform development. Specifically, the gully head is the catchment area at the intersection of the valley line and shoulder lines (Figure 6). The shoulder line is the boundary between gentle uphill and steep downhill terrain surfaces. The range of the gully head is from the intersection of the shoulder line and valley line to the point of maximal curvature changes along the valley line (the point corresponding to the maximum change in the slope) [37]. In this study, four steps are used to extract the gully head (Figure 7).

First, the criteria of extraction and validation are established through field surveys. Estimation of the typical topographic features of gullies is conducted based on laser rang-finders. For further investigation, the position of the gully head is measured by the total station. All of the spatial data obtained from the field surveys are digitized to known coordinates (WGS 84) in the GIS database and corresponding DEM. The data from field surveys contribute to verifying and improving the accuracy of the gully head extraction.

Second, an iterative channel deepening algorithm [13] is used to extract the gully system to eliminate parallel gullies and pseudo gullies. The ArcGIS Hydro analysis tools are used to extract a full gully network via the following approach: 1. Fill the sinks of the DEM; 2. set a larger flow accumulation threshold with 100 to generate the initial gully network; 3. use the modified DEM to calculate the flow direction and cumulative volume to extract a new gully network; 4. recodify the DEM according to the burn-in algorithm to eliminate parallel channels; 5. repeat steps 3 and 4 until all the small channels of the gully head are identified.

The third step is precisely positioning the gully head. Studies have shown that the location of the gully head is closely related to the shoulder line. The P-N terrain segmentation algorithm [38] is used to revise the shoulder line. The method to modify the gully head tracks the upstream unit with the largest flow accumulation (cell size: 5 × 5 m²) until one reaches the shoulder line [13].

Finally, the effectiveness of the method is assessed. The reference data for the assessment are described by the aerial image and field investigation, while the aerial image dataset may be affected by individual subjectivity. Therefore, the extracted results are mainly evaluated by field survey reference data in small sampling areas or manual identification. Digital orthophoto maps are used only as an auxiliary evaluation method.

3.3. Monitoring the Gully Head Erosion with SBAS-InSAR

To obtain a comprehensive and detailed description of the surface deformation to characterize and monitor the gully head erosion of the study area, 22 ascending L-band ALOS/PALSAR datasets (processed at Level 1.0) from 1 January 2007 to 27 February 2011 are utilized (

Table 1. Parameters of the ALOS PALSAR images in the research area.

Image No.	Incident Angle (°)	Acquisition Date	Polarization
1	38.708	1 January 2007	HH
2	38.705	16 February 2007	HH
3	38.702	4 July 2007	HH
4	38.709	19 August 2007	HH
5	38.698	4 October 2007	HH
6	38.727	19 November 2007	HH
7	38.693	4 January 2008	HH
8	38.721	19 February 2008	HH
9	38.685	5 April 2008	HH
10	38.694	21 May 2008	HH
11	38.705	6 July 2008	HH
12	38.714	6 January 2009	HH
13	38.714	21 February 2009	HH
14	38.706	9 July 2009	HH
15	38.701	24 August 2009	HH
16	38.696	9 October 2009	HH
17	38.720	9 January 2010	HH
18	38.723	24 February 2010	HH
19	38.709	12 July 2010	HH
20	38.714	12 October 2010	HH
21	38.699	12 January 2011	HH
22	38.703	27 February 2011	HH

). The frame number is 0710, and the tract number is 162. The L-band is the transmitting and receiving microwave frequency of the PALSAR, which operates at a wavelength of 23.6 cm with a temporal resolution of 46 days, and the incidence angle of the acquisitions is approximately 38.7°. The shuttle radar topography mission (SRTM) DEM with a resolution of 30 m was applied to remove topographic phase contributions. A multilook factor of 10 (two pixels in range and five pixels in azimuth directions) is applied to generate InSAR images at a spatial resolution of 15 m in both the range and azimuth directions. The location of reference area for processing is shown in Figure 12. The reason for choosing this area as a reference area is due to the fact that it is a stability area with no deformation and no atmospheric interference.

The topography of the gully head area is complex, and has low coherence [39]. The SBAS-InSAR is an effective method for measuring the surface deformation of low-coherence areas [40]. SBAS-InSAR is considered a multitemporal InSAR technique that uses the time series of multiple master InSAR images to effectively reduce atmospheric artifacts and residual terrain during the deformation retrieval [18]. The surface deformation of the gully head is determined by the small baseline interference pairs of SBAS-InSAR. Meanwhile, DEM data are combined with the SBAS-InSAR results to reliably identify gully head erosion.

The main errors of the SBAS-InSAR process arise from data phase unwrapping and atmospheric phase error. The data processing procedure is as follows:

The first step is to generate the connection graph. The spatial-temporal baseline threshold is determined to divide the data sets into different subsets (spatial baseline of 3000 m; temporal baseline of 360 days), and generate 76 high-quality interferometric pairs (Figure 8).

The process of generating the interferogram is the most critical step. Goldstein filters are used in this process to improve the phase unwrapping accuracy and measurement accuracy to maximize the signal-to-noise ratio of the interferograms. External DEM data are used to remove the topographical phase. The coherence threshold of 0.35 with the smallest cost flow is used for unwrapping.

The third step is to generate the surface deformation of a gully head. The residual topographic phase and the satellite orbit errors for the remaining pairs are removed by considering ground control points (GCPs). The GCPs should be in the best possible unwrapped phase area, and far from the deformation area. To remove the residual topography, a cubic inversion model is used. A spatiotemporal filter alleviates the effect of atmospheric delay (space filter window: 1800 m; time filter window: 365 days), and acquires the average deformation rate and the accumulated land deformation along the line of sight (LOS). The singular value decomposition is used to obtain the least squares solution, and then approximate the nonlinear time series deformation.

The fourth process is to produce the surface deformation along the slope angle. The values of the original surface deformation along the LOS obtained by SBAS-InSAR can be considered. However, the deformation rate in the direction of the LOS is insufficient to represent the true slope deformation in mountainous areas. The rate of deformation along the slope (V_slope) is defined in Equation (4), based on the deformation rate of the LOS direction (V_LOS) [27]:

$${\mathrm{V}}_{\mathrm{slope}}={\mathrm{V}}_{\mathrm{LOS}}/{\mathrm{n}}_{\mathrm{LOS}}.{\mathrm{n}}_{\mathrm{slope}},$$

$${\mathrm{n}}_{\mathrm{LOS}}=(-\sin \mathsf{\theta }\cos {\alpha }_s, \sin \mathsf{\theta }\cos {\alpha }_s, \cos \mathsf{\theta }{)}^{-1},$$

(4)

$${\mathrm{n}}_{\mathrm{slope}}=(-\sin \mathsf{\alpha }\cos \mathsf{\phi }, \sin \mathsf{\alpha }\cos \mathsf{\phi }, \cos \mathsf{\phi }{)}^{-1},$$

where $\mathsf{\phi }$ is the slope angle; $\mathsf{\alpha }$ is the slope aspect; $\mathsf{\theta }$ is the radar incident angle; ${\alpha }_s$ is the angle between the satellite’s orbital direction and true north (which is the flight direction of the radar satellite); the ascending image is negative; and the descending image is positive.

Finally, to achieve the result in the WGS 84 coordinate system, the surface displacement data are geocoded.

It is possible to measure the surface deformation of gully heads over large areas with SBAS-InSAR. The gully head erosion rate is calculated using the gully surface deformation and gully head area. Then, a regression model is formed with the gully head erosion rate results, and the gully head erosion can be assessed. Based on the obtained surface deformation along the slope angle ($d$), the gully head erosion ($\mathrm{g}$) is calculated (Figure 9). The gully head erosion in Figure 9 is highlighted in light green, which is defined by an integral value of every gully head surface deformation multiplied by the pixel size. In other words, the amount of gully head erosion is the cumulative sum of the deformation in the pixel area of gully head. The pixel area for every gully head can be determined (see Section 3.2), and each pixel area multiplied by the corresponding surface deformation along the slope direction (SBAS-InSAR values in every pixel) is the gully head erosion of the pixels. The amount of erosion for each pixel in the gully head area is summed to determine the gully head erosion (g). The gully head erosion rate (G) in this study is the ratio between the gully head erosion (g) and gully head area (A) (Equation (5)).

$$\mathrm{G}=\mathrm{g}/\mathrm{A}$$

(5)

3.4. Evaluating the Gully Head Erosion Rate (GHER)

To quantify the univariate relation between the gully head erosion and environmental characteristics (terrain attributes and soil type factors), a correlation matrix (

Table 2. Spearman’s rank correlation coefficient among terrain attributes (n = 415).

	GHER	Sl	Ca	TWI	SPI	LS	Curvpl	Curvpr	Sa
GHER	1
Sl	0.768 *	1
Ca	0.895 *	−0.670 *	1
TWI	0.096	0.273 *	−0.273 *	1
SPI	−0.071	0.029	0.035	−0.002	1
LS	0.901 *	−0.897 *	−0.900 *	−0.030	0.030	1
Curvpl	−0.016	0.001	0.049	0.120 *	−0.030	0.020	1
Curvpr	−0.034	0.039	0.044	−0.004	−0.053	0.039	−0.030	1
Sa	−0.242 *	0.112 *	0.233 *	0.118 *	0.066	0.190 *	0.001	−0.088	1

* Significant correlation at the p < 0.05 level.

) is generated. The relationship between the gully head erosion rate and the effective factors was determined by the regression analysis.

3.5. Evaluating the Validity and Efficiency of the Model

The coefficient of determination (R²) is used to evaluate the validity, goodness of fit, and efficiency of the model. R² reflects the proportion of dispersion described by the prediction (Equation (6); Nagelkerke et al., 1991) [41]. A large R² (value closer to 1) indicates that the model is a strong predictor [42], and a value greater than 0.5 is considered acceptable [9].

$${\mathrm{R}}^2={\left({{\displaystyle \sum }}_{i=1}^n(g_i^m-{\overline{g}}_i^m)(g_i^e-{\overline{g}}_i^e)/\sqrt{{{\displaystyle \sum }}_{i=1}^n(g_i^m-{\overline{g}}_i^m)}\sqrt{{{\displaystyle \sum }}_{i=1}^n(g_i^e-{\overline{g}}_i^e)}\right)}^2,$$

(6)

where $g_i^m$, $g_i^e$, ${\overline{g}}_i^m$, and ${\overline{g}}_i^e$ are the measured, estimated, mean measured, and mean estimated gully head erosion rates (m³/m²/year).

4. Results

4.1. Locating of the Gully Heads

A gully system map is composed of a shoulder line, a gully network, and a gully head (Figure 10). In total, 415 gully heads are extracted in the research area. The gully network without pseudo channels and parallel channels is extracted after 14 iterations. Extracting the gully heads is the first important step in analyzing the gully head erosion.

This macroscopic analysis is feasible, but an accuracy assessment is necessary. The high accuracy of the extracted gully system is confirmed by a visual assessment and image-based interpretation. Field surveys are the most fundamental method of validation. Therefore, some of the valley lines and gully heads were validated by a field survey (Figure 11). The results show that the gully heads are effectively extracted with the DEM. However, some small gullies are not extracted since the resolution of the DEM data is 5 m.

4.2. Gully Head Erosion Rates

4.2.1. SBAS-InSAR Results

The SBAS-InSAR results are illustrated in Figure 12. The negative values correspond to the gully head erosion. The positive values indicate that soil accumulates in the gully head.

The surface deformation (${\mathrm{V}}_{\mathrm{slope}}$) ranges from −73 to 68 mm/year, the mean surface deformation is −13.5 mm/year, and the variance of the mean surface deformation is 11.0 mm/year. Evidently, most parts of the river basin gully bottoms and flat terrain areas are stable (Figure 12). On the one hand, a flat terrain does not easily form gullies. On the other hand, human activities, such as road and building construction, interfere with low-altitude areas. The severe deformation is homogeneously distributed in the gully area. Approximately 33.04% of gullies have gully heads with an average deformation rate above 50 mm/year, and are defined as rapid gully head erosion areas. The white pixels represent a lack of data due to temporal decorrelation. Therefore, it may not be possible to detect some areas of deformation that suddenly collapse.

4.2.2. Accuracy Assessment of the SBAS-InSAR Results

In total, 27 points in flat areas (slope is approximately 0°) are used to validate the effects of SBAS-InSAR, as shown in Figure 13. The DEM values describe the vertical direction. The SBAS values of the LOS were projected in the vertical direction to be consistent with the DEM values to assess the SBAS results. The considered DEM values are from January 2007 to 2008, while the considered SBAS values are from January 2007 to 2008. Figure 14 compares two sets of stable points for the SBAS-InSAR and DEM in the same time interval. The results indicate that the stable points in both 2007 and 2008 have a reasonable consistency. Furthermore, R² and the root mean squared error (RMSE) of the two sets are 0.969 and 6.11 mm (the annual rate is 7.60 mm), respectively. The mean and maximum differences are 0.64 and 3.50 mm, respectively. The data are consistent in other years, which mean that the SBAS-InSAR results in other years are also suitable. Therefore, the monitoring of SBAS-InSAR in the research areas has produced relatively reliable results.

4.3. Spatial Factors That Influence the Gully Head Erosion Rate

Many simulations and studies of regional gully head erosion in the Loess Plateau have based their analyses and evaluations on the assumed topography, climate, soil type, and land use [31]. The increase in gully head size mainly depends on the catchment area, slope angle, topographical wetness index, slope aspect, stream power index, curvature, and precipitation [20]. We analyze the main factors that affect the change in gully erosion in the research area as follows.

4.3.1. Precipitation

Precipitation is one of the main forces that cause the gully head erosion [43]. Gully head erosion likely occurs due to continuous and torrential rainfall [44]. The statistics of precipitation and temperature data for the region in 2007–2011 are shown in Figure 3. The temperature and precipitation are basically positively correlated, with the most rainfall and highest temperature occurring in July-August and the lowest precipitation occurring in January–February each year, which is generally snow. Although there are variations in 2007–2011 and the rainfall data present seasonal and annual changes, the rainfall is generally stable. Therefore, the changes in rainfall may not account for the changes in the observed rates of gully head erosion and trends. In addition, there is only one rainfall monitoring station in the research area, so the spatial distribution of rainfall cannot be described.

4.3.2. Terrain Attributes

Topography is a significant factor that induces gully head erosion. The main topographic factors that we consider are Sl, Ca, TWI, SPI, LS, Curv_pl, Curv_pr, and Sa, which play very important roles in initiating and developing gully erosion. Spearman’s rank correlation was applied to describe the relationships among the gully morphological parameters (Table 2). A high gully head erosion rate is associated with large Sl, Ca, and LS values, as indicated by the very strong positive correlations. The correlation coefficients among the terrain attributes are significant (p < 0.05; for example, r = 0.895 among both the gully head erosion rate (GHER) and Ca). The most significant factors are LS, Ca, Sl, and Sa. The r value of 0.901 indicates an increase in GHER with increasing LS, while the r value of −0.242 indicates an increase in GHER with decreasing Sa. Thus, a higher GHER requires larger Sl, Ca, and LS values in the research area.

For a significant r among the Sl, Ca, TWI, LS, and GHER, a prediction model can be constructed by the fitting curve of the rate of gully head erosion. The least squares method of fitting is adopted, and a power function form (Equation (7)) is used. For the prediction of the trend term, the parameters of the power function are used and the results are shown in Figure 15 and

Table 3. Single factor parameters based on power function fitting.

Type	Parameter			Accuracy (R2)
	a	b	c
Sl	−0.256	1.324	0	0.665
Ca	−0.033	1.124	−1.946	0.758
TWI	−0.095	2.42	−5.44	0.519
LS	−0.857	1.053	0	0.812

.

$$\mathrm{D}\left(\mathrm{x}\right)={\mathrm{ax}}^{\mathrm{b}}+\mathrm{c},$$

(7)

Figure 16 shows the GHER as a function of class for Curv_pl, Curv_pr, SPI, and As. Excluding the lower class (SPI (<1); Curv_pl (<−2); Curv_pr (−0.5)), GHER when SPI, Curv_pl, and Curv_pr increased. Moreover, for the class of SPI (−1~1), Curv_pl (−2~−1), Curv_pr (−0.5~0.5), and Sa (E), the results indicated the broadest GHER range, which implies that the gully heads exhibited both high and low erosion rates.

4.3.3. Soil Type

GHER does not vary significantly among the soil types and parent materials in the research area (Figure 17). The highest GHER is observed for Q₃ (mean = −30 m³/m²/year), followed by Q₂ (mean = −24 m³/m²/year) and Q₁ (mean = −20 m³/m²/year). At p < 0.05, however, these differences are negligible. This result is identical to that of Frankl et al. (2013) [45]. Additionally, this result is probably due to the limited difference in soil types across the sample region [45].

4.4. Estimating the Erosion Rates of the Gully Head

In this study, the highest correlations of the fitted GHER curves with Sl, Ca, TWI, and LS are explored. The coefficient results also show that the high correlations represent a strong ability to predict the gully head erosion. Based on the correlation equation of every significant factor, the prediction model of the gully head erosion is developed by the following relationship function (Equation (8)).

$$\mathrm{G}\left(\mathrm{x}\right)=\mathrm{D}\left(\mathrm{Sl}\right)\ast \mathrm{D}\left(\mathrm{Ca}\right)\ast \mathrm{D}\left(\mathrm{TWI}\right)\ast \mathrm{D}\left(\mathrm{LS}\right)\ast \mathrm{D}\left({\mathrm{x}}_{cr}\right)$$

(8)

Here, the correction values ($\mathrm{D}\left({\mathrm{x}}_{cr}\right)$) are related to SPI, Curv_pl, Curv_pr, and Sa. The fitting relationship function of each factor (Sl, Ca, TWI, LS) are plugged into the model. Seventy percent of the gully heads are used to estimate the model coefficient, the correction value method is preliminarily tested, and the gully head erosion rate model is presented as follows (Equation (9)):

$$\mathrm{G}=1.5{\mathrm{Sl}}^{0.256}{\mathrm{Ca}}^{0.095}{\mathrm{TWI}}^{0.033}{\mathrm{LS}}^{0.857}$$

(9)

Based on the linear regression of the GHER values measured from the SBAS-InSAR results and gully head erosion model, we randomly selected 30% of the gully heads to calculate the error of the GHER results, and the linear fit is good (R² = 0.889) (Figure 18).

5. Discussion

5.1. Characterization of Gully Head Activities

The GHER distribution is half of the normal distribution (Figure 19). The average estimated GHER value is −7.5 m³/m²/year, and the distribution is concentrated in the range of -15~0 m³/m²/year, which is within the range of Xu et al. results [3], who reported GHER values as high as −3.8 to −52.87 m³/m²/year in the loess area of China. This is a relatively high gully head erosion rate, compared to those in the literature. Frankl et al. [46] recorded average gully head erosion rates of −5.2 ± 5 m³/m²/year in the semiarid highlands of northern Ethiopia, and Vandekerckhove et al. reported an overall rate of 4.0 m³/m²/year in southeastern Spain [47]. However, slightly higher gully head erosion rates have been observed at steep slopes in China, due to the influence of a more concentrated runoff caused by urban expansion and infrastructure development, more extreme storms, and thicker loess [3]. A potential reason for the high erosion rates of gully heads in the current study is the presence of large catchment areas, which cause the runoff to violently wash the gully head and cause a serious gully head erosion during the heavy rainfall.

Based on the rate of gully head erosion, we classified gullies as active gullies and inactive gullies [48]. Active gullies have average erosion rates above −15 m³/m²/year. Figure 17 shows that 152 gully heads have erosion rates above −15 m³/m²/year, which account for 33.04% of all gully heads. The locations of active gully heads are shown in Figure 12. This information can provide a reference for the project site selection in engineering construction (the construction of roads, high-speed railways, etc.) so that active gullies can be avoided.

5.2. Gully Volume Model vs. GHER Model

There are a few models that can be used to quantitatively evaluate gully erosion. Several studies have studied the relationship between the volume and length of a gully using a power equation of the form of $\mathrm{V}={\mathrm{aL}}^b$ [3,17,45,49,50,51]. Thompson (1964) [29,52] used a function focused on multiple factors (relationship between gully head and Sl, rainfall (P), Ca, and clay content (E)) based on short-term observations. However, in the study of the effects of topographical and environmental factors on the gully development, e.g., with the model $\mathrm{V}=0.15{\mathrm{Ca}}^{0.49}{\mathrm{Sl}}^{0.14}{\mathrm{P}}^{0.74}E$, the coefficient and exponents are uncertain for long-term estimates. The gully volume was measured according to this equation using the gully areas obtained from high-resolution maps, and the gully erosion volumes during the assessment period were summarized. Li et al. [17] observed that a strong relationship exists between the gully volume and area. Compared to the relationship between the volume and length of a gully, the relationship between the volume ($\mathrm{V}$) and area (${\mathrm{A}}_g$) of a gully was more prominent in the Loess Plateau of China, and the equation of $\mathrm{V}=0.2762{\mathrm{A}}_g{}^{1.3971}$ was proposed to estimate the gully volume. This model could better predict the gully volume in a larger area. However, this approach requires very high-resolution satellite data, and it is relatively labor intensive and time consuming to use in very large areas.

In this study, we explored the relationship among terrain attributes and the gully head erosion rate, and the gully head erosion rate was measured with SBAS-InSAR. The model is verified to be accurate for the Loess Plateau. The erosion rates of gully heads were directly calculated using our assessment model over a long period and over a large area, and it was found that the SBAS-InSAR and regression model considering terrain attributes can be used to accurately measure the gully head erosion.

6. Conclusions

Based on an application of the SBAS-InSAR technique in the Chinese Loess Plateau, a new model to assess the gully head erosion is proposed. The main conclusions are as follows:

(i)

First, 415 gully heads are extracted by the DEM, and this approach is proven to be feasible by comparing the results with field survey data and digital orthophoto maps. In addition, the ALOS PALSAR data and the SBAS-InSAR method are used to obtain favorable results in monitoring the surface deformation of gully heads, and the SBAS-InSAR results suggest that most of the gully bottoms and flat terrain areas are more stable than the gully heads.

(ii)

A simple regression analysis estimates the relation among the erosion rates of gully head erosion, terrain attributes, and soil type. Gully head erosion is strongly positively related to the topographical factors of the slope angle, catchment area, topographic wetness index, and slope length. In contrast, the soil type does not significantly affect the gully head erosion.

(iii)

A new framework based on spatial factors ($\mathrm{G}=1.5{\mathrm{Sl}}^{0.256}{\mathrm{Ca}}^{0.095}{\mathrm{TWI}}^{0.033}{\mathrm{LS}}^{0.857}$, R² = 0.889) is proposed to model the gully head erosion. One of the main advantages of combining the erosion rate of the gully head with terrain factors is to simplify the evaluation of the gully head erosion over large areas and long periods, particularly when only environmental information obtained from remote sensing data is available.

(iv)

Gully head protection should focus on controlling the terrain attributes of the slope angle and catchment area (e.g., reducing the catchment area of a gully head is an effective method to decrease the gully head erosion). The results of this study can be used for the control and risk assessment of the gully head erosion.