Topographic Spatial Variation Analysis of Loess Shoulder Lines in the Loess Plateau of China Based on MF-DFA

The Loess Plateau in China is internationally known for its unique geographical features and has therefore been studied by many researchers. This research exploits the regional differences in the spatial morphological characteristics of Loess shoulder lines in different landform types as an important basis for geomorphological regionalization. In this study, we used ensemble empirical mode decomposition (EEMD), multi-fractal detrended fluctuation analysis (MF-DFA), and detrended cross-correlation analysis (DCCA) to analyze topographic data series extracted from shoulder lines. Loess shoulder-line land variations series data from the Suide, Ganquan, and Chunhua areas on the Loess Plateau were selected and a combination of the two above-mentioned methods was used to study land variations at these three sample sites. The results revealed differences in the topographic variations of the multi-fractal characteristics and the topographic spatial variation in the Loess shoulder line of the three sample sites. Furthermore, the extent to which the results were affected by noise and the analysis scale differed among the three areas.


Introduction
The Loess Plateau of China has drawn the attention of many researchers worldwide for its unique geographical condition, geomorphologic research values, loess landform feature, and its special natural and cultural landscapes [1,2].Many previous studies about the Loess Plateau have been carried out to solve scientific problems regarding its formation process [3,4], landform shape [5][6][7], gully development [8][9][10], soil erosion [11,12], environmental evolution [13,14] etc., and many outstanding research results have been achieved.In recent years, the development of Digital Elevation Model (DEM) applications and Geographic Information System (GIS) spatial analysis techniques has greatly promoted the study of the land surface analysis of the Loess Plateau.Digital land surface analysis has been used to extract the basic terrain derivatives (e.g., slope and aspect) [15,16] and to study complicated applications (e.g., paleo-topographic reconstruction and landform classification) [17][18][19][20][21] of the Loess Plateau.The terrain points [22][23][24], lines [25][26][27], and surfaces [28][29][30] of the Loess Plateau have been well studied in terms of extraction methods based on DEMs.However, these studies mainly focused on research with the view of developing and evaluating technologies and methods, and there is a lack of research and discussion on the role of terrain data in the study of the topography and geomorphology of the Loess Plateau.Tang et al. [1] established the theory of the slope topo-sequence, which they successfully applied to study the Loess landform.The corresponding slope topo-sequence, which showed an obvious spatial distribution, of any of the different Loess landforms could be found [31].The topo-sequence of hypsometric integrals (HI) was applied to the loess watershed landform in an area of the Loess Plateau affected by severe soil erosion, based on a digital land surface analysis and geo-informatics method named Tupu [32].Furthermore, an existing method has been put forward to study the spatial distribution based on the landform planation index (LPI), which is the terrain derivative extracted from DEMs [33].In spite of the fact that these methods realized a quantitative analysis of the macroscopic differentiation rules of the Loess Plateau, these methods are not sensitive to the local topographic changes in the Loess Plateau, which cover up the local variation of the terrain.
The Loess shoulder-line is one of the most significant topographic structures in the Loess Plateau, which separates the complicated physiognomy of the plateau distinctly into positive (inter-gully area) and negative terrains (inner-gully area) (P-N terrains) in terms of the formation and topographic feature of the landforms.This shoulder-line plays important roles in the research of soil and water conservation.As the most vital boundaries, extraction methods [34] and efficient calculations [35,36] have mainly been applied to the Loess shoulder lines in many studies.Few researchers have applied the quantitative index of the shoulder lines to study the spatial distribution of the Loess Plateau land at present.The shoulder lines formed and developed at the slope scale, meandering in the whole watershed, and carved out the spatial pattern of Loess gullies and Loess tableland, Loess ridges, and Loess hills.The regional difference in the spatial morphological characteristics of Loess shoulder lines in different landform types has become an important basis for geomorphological regionalization.
Fractal theory has been used to address the limitation of traditional quantitative methods on geomorphic features, and has opened up a new approach to generally describe the geomorphic features of a drainage basin and its geomorphic evolution [37,38].More specifically, fractal theory has provided a way to quantitatively describe self-similar or self-affine landforms [39].The classic application of fractals to geomorphology focused on line elements such as coastlines, rivers, and faults [18].Based on grid DEMs, the fractal dimension was shown to characterize more than mere terrain "irregularity" or "roughness" [40]; instead, it can also be employed to describe the changes in the variability of the topography in relation to the distance [41].The relationship between the fractal parameters and the physical processes needs further research.For a loess landform, the fine fractal structure of a single watershed varies for different landform types and developmental processes.Thus, an analysis of the geometric singularity distribution of loess shoulder lines in different loess landforms to obtain the fractal parameters of these shoulder lines could provide theoretical and methodological support for further research on scientific and reasonable classification of loess geomorphic regions.
The research described in this paper is based on elevation data, which was obtained from the actual direction of the curve at a certain distance along shoulder lines extracted from DEMs.The elevation data, as with most time series data, also consist of nonlinear and non-stationary data capable of reflecting the topographic spatial variation and relief of the area surrounding the shoulder lines.The main difference between the time series data and elevation data is that the two kinds of data are based on the time and sampling distance as independent variables, respectively.Due to a high degree of similarity within time series, it is feasible to use the elevation sampling data of shoulder lines to apply to land surface analysis.
This research mainly focuses on the introduction of the multi-fractal spectrum exponent and Hurst exponent in order to describe the topographic variation law of different landform types that can reflect the overall macro-and local micro-pattern of the land surface of the Loess Plateau.The research results provide quantitative information in support of the loess landform classification.

Study Area
The hilly area of the Loess Plateau of China was taken as the study area.This area was developed on a base of ancient landform consisting of Mesozoic bedrock.The Cenozoic laterite and Loess layers are cut by the flow of water, with soil erosion forming the complex gully systems.The selected area contains three types of landform: loess tableland, loess ridge, and loess hill, which are studied by our proposed method to understand their distributions and differences.Distinguishing these differences required the study sites to be confined to those with the most typical geomorphologic characteristics and without a mixture of landforms.Therefore, three sample areas were selected: (1) loess hill; (2) loess ridge; and (3) loess tableland.Geographic descriptions of each of the sample sites are provided in Table 1.

Topographical Data
All data used in this study are from the revised shoulder lines of these three sites and were extracted by DEMs with a resolution of 5 m (Figure 1).The DEMs used in this study were obtained from the National Administration of Surveying, Mapping, and Geoinformation of China, and were generated from the 1:10,000 topographic maps with a 1 m contour interval.Firstly, the paper topographic maps were drawn by ground-based surveying and then scanned into the computer with a geometric correction.Then the contour lines were digitalized and interpolated into a triangulated irregular network (TIN).Finally, the grid DEMs were interpolated by the TIN after ordinary manual editing in TIN to correct the mistakes.The shoulder line delineation could be achieved by a manual interpretation process with the automatic result derived by DEMs since the error of automatic results may lead to wrong analytic conclusions.Automatic extraction could be performed by using the GVF Snake algorithm [25].Among the direction of the shoulder lines, elevation series were sampled with 5-m distance intervals.The number of elevation series in the three sites are 16,000 (Suide), 10,000 (Ganquan), and 7000 (Chunhua).

EEMD
Ensemble empirical mode decomposition (EEMD) is an adaptive analysis method for time series proposed by Wu [62], and is performed with the assistance of data noise.The addition of white noise to the signal to be decomposed allows the combination of signal and noise to be regarded as a whole, and then the Intrinsic Mode Functions (IMFs) with modal consistency can be obtained.Since the distribution of the white noise spectrum is uniform, the signals of different scales are automatically distributed to the appropriate reference scale when the signals are added to the white noise, which has a background with the same time-frequency distribution.Because the mean value of this noise is zero, noise of different frequencies will cancel each other after multiple averaging operations, and the result of the integrated mean can be directly used as the final result.The algorithm can be described as follows: Step 1: Add white noise to the signal to be decomposed where x(t) is the true signal and n(t) is the noise.
Step 2: Perform the EMD decomposition into a series of IMFs j c using the signals with noise

EEMD
Ensemble empirical mode decomposition (EEMD) is an adaptive analysis method for time series proposed by Wu [62], and is performed with the assistance of data noise.The addition of white noise to the signal to be decomposed allows the combination of signal and noise to be regarded as a whole, and then the Intrinsic Mode Functions (IMFs) with modal consistency can be obtained.Since the distribution of the white noise spectrum is uniform, the signals of different scales are automatically distributed to the appropriate reference scale when the signals are added to the white noise, which has a background with the same time-frequency distribution.Because the mean value of this noise is zero, noise of different frequencies will cancel each other after multiple averaging operations, and the result of the integrated mean can be directly used as the final result.The algorithm can be described as follows: Step 1: Add white noise to the signal to be decomposed where x(t) is the true signal and n(t) is the noise.
Step 2: Perform the EMD decomposition into a series of IMFs c j using the signals with noise where r m is the residue of data X(t), after m number of IMFs are extracted.
Step 3: Repeat Steps 1 and 2 using a different type of noise each time.
Then decompose Step 4: Calculate the mean value of the IMFs The final result is obtained Each IMF component obtained by the decomposition, irrespective of whether it is pure noise or the original sequence component with a clear physical significance, can be judged by a Monte-Carlo significance test proposed by Wu et al. [63].IMFs that cannot pass the test with a confidence limit level of 99% can be seen as noise.Here, it is important to note that the noise in this particular situation (referring to the local errors in DEMs caused by sampling or some abnormal data) is totally different from the white noise used in EEMD (which is simply a mathematic support for data decomposition).Finally, all remaining IMFs are reconstructed into an entirely new data series and the fractal characteristics analyzed by the following MF-DFA method.

MF-DFA
The MF-DFA is an effective measure to reflect the developmental process of a system [42].For the sum standardization of the fixed length time series, the MF-DFA method consists of six steps.The first three steps are essentially the same as in conventional DFA methods.Suppose that {x i } is a numeral series of length N, and this series obeys some compact support, i.e., x i = 0 means an insignificant fraction of the values only.The procedure of the MF-DFA method is as follows: Step 1: Calculate the cumulative difference sequence of {x i } where <x> is the mean value of {x}.
Step 2: Apply backward division of the series Y(i) into N s = int(N/s) individual segments with an equal length scale s.Generally, there will remain a short part at the end of the series since the length N of the series is often not a multiple of the given scale s.Therefore, a forward division procedure should be repeated so that 2Ns segments can be obtained altogether.
Step 3: Calculate the local trend for each 2Ns segments by the least-squares procedure.Then determine the variance where y p (i) is the fitting polynomial of segment p, p is the number of segments, s is the length scale and, and F is the variance.
Step 4: Calculate the qth order fluctuation function by averaging over all segments.
(when q is any real number that is not zero) (10) Step 5: Use several different q to determine the scaling exponent of fluctuation function F q (s).If there is long range correlation in series {x i }, F q (s) will exhibit an impact of power law relationship on length scale s, i.e., log where C is a constant, Hq is the qth scaling exponent of fluctuation function, i.e., the Hurst index.
Since the Hurst exponent Hq can only measure the fluctuation of an elevation series from a global aspect, the local fluctuation trend could be indicated by the local Hurst exponent Ht.Ht could be estimated by performing the Hq calculation procedure in a local series segment around the given series location.Compared with Hq, Ht takes advantage of identifying the instant structural changes in the elevation series.
Step 6: Calculate the qth mass exponent tq, the qth singularity exponent hq and qth singularity dimension Dq.In Equation ( 7), tq is another measurement for scaling exponents and Hq is the Hurst index, and the relationship between these two parameters is as below, The mass exponent tq could be applied for estimating the qth singularity exponent hq and the qth singularity dimension Dq as below: where hq is the singularity strength, used to describe the singularity degree of each interval in the time series.
Another way to characterize a multi-fractal series is ∆hq, which can be defined as below, ∆hq = hq max − hq min (16) where the maximum and minimum values of hq are expressed by hq max and hq min , respectively.Since the multi-fractal spectrum has either a left or a right truncation when using a pair of the qth-order of negative and positive values to generate the Hurst exponent, the width of the multi-fractal spectrum, ∆hq can indicate the homogeneity of the time series variation.A long tail to the left of the multi-fractal spectrum indicates that the series has a multi-fractal structure insensitive to local fluctuations with small magnitudes, whereas a long tail to the right indicates that the series has a multi-fractal structure insensitive to local fluctuations with large magnitudes.The larger ∆hq, the more uneven the time series distribution is reflected, and the more obvious the multi-fractal feature.The above statistical and analysis work was carried out by using the software MATLAB2014 [64].

DCCA
In the MF-DFA method, the length scale s is an important parameter for analysis.The scale effect assessment between the local Hurst exponent Ht and the elevation sequence is achieved by introducing a well-used method for long-term cross-correlations investigation [65], detrended cross-correlation analysis (DCCA).The procedure of the DCCA method consists of the following four steps: Step 1: Determine the new profiles from two series, Ht series {x i } and elevation sequence {y i }.
Step 2: Similar to Step 2 in the MF-DFA method, divide the profiles {x k } and {y k } into n = 2Ns non-overlapping segments with equal length n.
Step 3: Calculate the detrended covariance as below: where x k,i and y k,i are local fitted (least-squares fit) trends of each of the divided segments.If F(n) behaves as a power-law function of n, i.e., The cross-correlation exponent α can be measured by the slope of the log-log plot between F(n) and n.If α > 0.5, the long-term cross-correlations between two series are persistent, and vice versa.If α = 0.5, there is a non-cross-correlated relationship between the two series, namely a change in one series has no effect on the behavior of the other series.

Multi-Fractal Analysis of Loess Shoulder Line
Figures 2a, 3a and 4a show the q-order scaling function F(q) of the study sites.According to Equations ( 7)-( 9), the tendency by the least-squares method between the scaling function and the scale can be simulated by MATLAB software.Here, q is set to −5, 0, and 5; the scale is 16 to 1024; and the scaling function uses 2 as the basis of the logarithm.From the figures, we can see an increase in the trend between F(q) and q, and a decrease in the trend between H(q) and q.Therefore, we can conclude that the elevation series of the Loess shoulder line has a multi-fractal characteristic.
Figures 2b, 3b and 4b show the q-order Hurst exponents of each of the study sites according to Equation (10).We see a clear decrease in the trend between Hq and q, which confirms our argument that the elevation series of the Loess shoulder line has a multi-fractal characteristic.For the sample series, when q increases from −5 to 5, Hq of Suide decreases from 0.9729 to 0.6069, with the extent 0.366; Hq of Ganquan decreases from 1.2078 to 0.6375, with the extent 0.570; and Hq of Ganquan decreases from 1.0717 to 0.6669, with the extent 0.4048.Thus, Hq of all three areas are significantly not constant values, which reflect the strong multi-fractal characteristic.The range of Hq has a maximum value in Ganquan and a minimum value in Suide, which shows that the multi-fractal strengths of the three sites are different (Ganquan > Chunhua > Suide).In addition, in Suide and Chunhua, when q is at (−5, 5), Hq is at (0.5, 1), which indicates that the elevation series of Loess shoulder-lines has long-range correlation.This is also valid for Ganquan when q is at (−2.3, 5).
Figures 2c, 3c and 4c show the t-tq relationship of the three areas.We can see a clear increasing trend between tq and q, where tq is the upper convex curve.Thus, there is a nonlinear relationship between tq and q, which further confirms the multi-fractal characteristics of the Loess shoulder ISPRS Int.J. Geo-Inf.2017, 6, 141 of 19 line [66,67].The higher nonlinearity of the tq curves is translated into a wider f(a) dispersion with additional pronounced multi-fractal characteristics.
ISPRS Int.J. Geo-Inf.2017, 6, 141 of 19 line [66,67].The higher nonlinearity of the tq curves is translated into a wider f(a) dispersion with additional pronounced multi-fractal characteristics.q =-5 q =0 q =5 q =-5 q =0 q =5 data tq( 5 q =-5 q =0 q =5 q =-5 q =0 q =5 data tq( 5  line [66,67].The higher nonlinearity of the tq curves is translated into a wider f(a) dispersion with additional pronounced multi-fractal characteristics.q =-5 q =0 q =5 q =-5 q =0 q =5 data tq( 5 q =-5 q =0 q =5 q =-5 q =0 q =5 data tq( 5  Figures 2d, 3d, and 4d shows the Dq-hq relationship of the three sites according to Equation (12).All of the bell curves in Figures 2d, 3d, and 4d indicate strong multi-fractal characteristics.The top of the multi-fractal spectrum of Ganquan in Figure 3d and Table 2 is relatively flat (the spectrum width ∆ℎ is 0.8945), which shows that the elevation sampling data exhibit considerable fluctuation and non-uniform distribution.Contrary to this, the plots of Suide and Chunhua are relatively sharp (the spectrum widths ∆ℎ are 0.6479 and 0.6852), which shows a gentle fluctuation and uniform distribution.The Suide site is characterized by a long tail to the left in its multi-fractal spectrum (Figure 2d), which is indicative of strong soil erosion in Loess positive and negative terrains.For the Chunhua site, the multi-fractal spectrum (Figure 4d) displays a long tail to the right, which shows that the elevation sampling data sequence changes gently, and that the changes in the land surface along the shoulder line are relatively small due to its relatively weak soil erosion.In Ganquan, Figure 3d shows that the multi-fractal spectrum has an approximately symmetrical shape without obvious tails to the left or right, which indicates that the elevation sampling data sequence of the q =-5 q =0 q =5 q =-5 q =0 q =5 data tq( 5  Figures 2d, 3d and 4d shows the Dq-hq relationship of the three sites according to Equation (12).All of the bell curves in Figures 2d, 3d and 4d indicate strong multi-fractal characteristics.The top of the multi-fractal spectrum of Ganquan in Figure 3d and Table 2 is relatively flat (the spectrum width ∆hq is 0.8945), which shows that the elevation sampling data exhibit considerable fluctuation and non-uniform distribution.Contrary to this, the plots of Suide and Chunhua are relatively sharp (the spectrum widths ∆hq are 0.6479 and 0.6852), which shows a gentle fluctuation and uniform distribution.The Suide site is characterized by a long tail to the left in its multi-fractal spectrum (Figure 2d), which is indicative of strong soil erosion in Loess positive and negative terrains.For the Chunhua site, the multi-fractal spectrum (Figure 4d) displays a long tail to the right, which shows that the elevation sampling data sequence changes gently, and that the changes in the land surface along the shoulder line are relatively small due to its relatively weak soil erosion.In Ganquan, Figure 3d shows that the multi-fractal spectrum has an approximately symmetrical shape without obvious tails to the left or right, which indicates that the elevation sampling data sequence of the shoulder line in Ganquan varies, and the land surface changes are uneven in some areas along the shoulder line and gentle in other areas.

Spatial Variation of Loess Shoulder Line
Two scales were selected to analyze the spatial distribution of the Hurst exponent.Figure 5a is the normalization result of the elevation sampling data sequence, Figure 5b is Ht with scale = 7, and Figure 5c is Ht with scale = 17.In Figure 5b,c we can see that, when scale = 7, the maximum Ht is 1.3165 (Pn = 8020), the minimum is 0.5467 (Pn = 11800), and the mode value is 0.8292; when scale = 17, the maximum Ht is 1.1736 (Pn = 8020), the minimum is 0.5082 (Pn = 11,800), and the mode value is 0.8292.The minimum value indicates that the land there is fractured while the maximum is gentle.It should be noticed that the extreme point of Ht on both scales is the same, because the scale has little effect on the analysis in the case of Suide, or probably because, at Suide, the land surface of the Loess shoulder line is fractured but has an undulating land topography distribution.shoulder line in Ganquan varies, and the land surface changes are uneven in some areas along the shoulder line and gentle in other areas.

Spatial Variation of Loess Shoulder Line
Two scales were selected to analyze the spatial distribution of the Hurst exponent.Figure 5a is the normalization result of the elevation sampling data sequence, Figure 5b is Ht with scale = 7, and Figure 5c is Ht with scale = 17.In Figure 5b,c we can see that, when scale = 7, the maximum Ht is 1.3165 (Pn = 8020), the minimum is 0.5467 (Pn = 11800), and the mode value is 0.8292; when scale = 17, the maximum Ht is 1.1736 (Pn = 8020), the minimum is 0.5082 (Pn = 11,800), and the mode value is 0.8292.The minimum value indicates that the land there is fractured while the maximum is gentle.It should be noticed that the extreme point of Ht on both scales is the same, because the scale has little effect on the analysis in the case of Suide, or probably because, at Suide, the land surface of the Loess shoulder line is fractured but has an undulating land topography distribution.Figure 6a-c shows the results for Ganquan.In Figure 6b,c we can see that, when scale = 7, the maximum Ht is 1.8004 (Pn = 3900), the minimum is 0.6799 (Pn = 10,200), and the mode value is 1.035; and when scale = 17, the maximum Ht is 1.5420 (Pn = 2100), the minimum is 0.6838 (Pn = 1000), and the mode value is 1.035.The extreme point of Ht on both scales is totally different.This may be due to the fact that the scale has a strong effect on the analysis in Ganquan, which indicates that land morphology features vary unevenly.Figure 6a-c shows the results for Ganquan.In Figure 6b,c we can see that, when scale = 7, the maximum Ht is 1.8004 (Pn = 3900), the minimum is 0.6799 (Pn = 10,200), and the mode value is 1.035; and when scale = 17, the maximum Ht is 1.5420 (Pn = 2100), the minimum is 0.6838 (Pn = 1000), and the mode value is 1.035.The extreme point of Ht on both scales is totally different.This may be due to the fact that the scale has a strong effect on the analysis in Ganquan, which indicates that land morphology features vary unevenly.
Figure 7a-c shows the result for the Chunhua site.In Figure 7b,c we can see that, when scale = 7, the maximum Ht is 1.5751 (Pn = 1500), the minimum is 0.5807 (Pn = 6000), and the mode value is 0.8963; and when scale = 17, the maximum Ht is 1.3978 (Pn = 2700), the minimum is 0.5996 (Pn = 3400), and the mode value is 0.8963.The local Hurst exponent is mostly greater than the mode, while the extreme point of Ht of both scales is similar in that it is totally different, this may be due to the strong effect of the scale on the analysis of the Chunhua site, which indicates that the land morphology features are distributed unevenly.Figure 7a-c shows the result for the Chunhua site.In Figure 7b,c we can see that, when scale = 7, the maximum Ht is 1.5751 (Pn = 1500), the minimum is 0.5807 (Pn = 6000), and the mode value is 0.8963; and when scale = 17, the maximum Ht is 1.3978 (Pn = 2700), the minimum is 0.5996 (Pn = 3400), and the mode value is 0.8963.The local Hurst exponent is mostly greater than the mode, while the extreme point of Ht of both scales is similar in that it is totally different, this may be due to the strong effect of the scale on the analysis of the Chunhua site, which indicates that the land morphology features are distributed unevenly.

Scale Effect of Loess Shoulder-Line Topographic Variation
Figure 8a,b shows that the correlation coefficient α of the loess shoulder-line topography and local Hurst exponent of topographic variation of the Suide study site for the analysis of scale 7 and 17 are relatively small, i.e., 0.4163 and 0.3283, respectively.The DCCA parameter α decreases with increasing analysis scale, which indicates the analysis scale to loess shoulder-line topographic variation is not affected.In the different analysis scales, the DCCA correlation coefficient α of the Suide site is less than that of the Ganquan and Chunhua study sites.In all the analysis scales, the highest correlation coefficient is 0.4163, and the best analysis scale is 7 for the Suide study site.All

Scale Effect of Loess Shoulder-Line Topographic Variation
Figure 8a,b shows that the correlation coefficient α of the loess shoulder-line topography and local Hurst exponent of topographic variation of the Suide study site for the analysis of scale 7 and 17 are relatively small, i.e., 0.4163 and 0.3283, respectively.The DCCA parameter α decreases with increasing analysis scale, which indicates the analysis scale to loess shoulder-line topographic variation is not affected.In the different analysis scales, the DCCA correlation coefficient α of the Suide site is less than that of the Ganquan and Chunhua study sites.In all the analysis scales, the highest correlation coefficient is 0.4163, and the best analysis scale is 7 for the Suide study site.All of the DCCA correlation coefficients α are very low for the Suide site and this suggests that the local Hurst exponent is not reliable for representing the loess hill land surface variation.Figure 8c,d shows that the DCCA correlation coefficient α in the analysis scale 7 and 17 are 0.4246 and 0.4512, respectively, for the Ganquan study site.In this case, the DCCA parameter α increases when the analysis scale increases.This indicates the reliability of the local Hurst exponent that reflects the loess ridge land surface variation.For the Ganquan site, the optimal analysis scale is scale 17 when using the local Hurst exponent to describe the topographic variation.Figure 8e,f shows that the correlation coefficient α is more than 0.5 for the Chunhua site, at the analysis scale 7 and 17, the values of the correlation coefficient α are 0.5229 and 0.5055, respectively.The local Hurst exponent and loess shoulder-line land variation has a relatively strong correlation in the Chuanhua site and this indicates that the local Hurst exponent is very suitable to reflect the loess tableland topographic variation.According to Figure 9, at the different analysis scales, the DCCA correlation coefficient α is the highest for the Chunhua site, for which the correlation coefficient α is different from that of the other two sites.An increase in the analysis scale first increases and then decreases the correlation coefficient α, which reaches its maximum when the analysis scale is 11, which is the best analysis scale of the Chunhua site.In the loess shoulder-line land surface analysis of loess tableland, the selection of a suitable analysis scale has a great influence on the results of land surface analysis.Figure 10 shows the Ph-Ht probability distribution for all study sites.As demonstrated in Figure 11, the plot of the local Hurst exponent conforms to an approximate normal distribution and the mean Ht increases in the order: Suide < Chunhua < Ganquan.This result indicates that the topographic variations are relatively smoother at the Ganquan site, and relatively steeper at the Suide site.As demonstrated in Figure 11, the standard deviation of Ht decreases in the order: Chunhua > Ganquan > Suide.This result also indicates that Ht is concentrated around the mean of the Suide site, Ht is scattered around the mean of the Ganquan site, and Ht is the most disperse around the mean of the Chunhua site.Figure 10 shows the Ph-Ht probability distribution for all study sites.As demonstrated in Figure 11, the plot of the local Hurst exponent conforms to an approximate normal distribution and the mean Ht increases in the order: Suide < Chunhua < Ganquan.This result indicates that the topographic variations are relatively smoother at the Ganquan site, and relatively steeper at the Suide site.As demonstrated in Figure 11, the standard deviation of Ht decreases in the order: Chunhua > Ganquan > Suide.This result also indicates that Ht is concentrated around the mean of the Suide site, Ht is scattered around the mean of the Ganquan site, and Ht is the most disperse around the mean of the Chunhua site.

Influence Factors of Loess Shoulder-Line
The analytic result shows the MF-DFA methods achieve a much deeper understanding of the data structure than other data analysis methods.Since the multi-fractal analysis method can be successfully applied to many complex systems, it has great potential to be applied to the land surface structure of loess shoulder lines in the Loess Plateau.The influence factors and the physical meaning of the multi-fractal characteristics of the loess shoulder-line land elevation sampled must be considered to be important.
Figure 11a shows the multi-fractal spectrum obtained by MF-DFA analysis of loess shoulder-line land elevation sampling at the Chunhua, Ganquan, and Suide study sites.This figure shows that the three curves of the multi-fractal spectrum are very similar to the curve with a long tail to the right.Determining the factors causing these results is a problem worthy of further discussion with respect to the method of land surface analysis in the Loess gully areas.First, the

Influence Factors of Loess Shoulder-Line
The analytic result shows the MF-DFA methods achieve a much deeper understanding of the data structure than other data analysis methods.Since the multi-fractal analysis method can be successfully applied to many complex systems, it has great potential to be applied to the land surface structure of loess shoulder lines in the Loess Plateau.The influence factors and the physical meaning of the multi-fractal characteristics of the loess shoulder-line land elevation sampled must be considered to be important.
Figure 11a shows the multi-fractal spectrum obtained by MF-DFA analysis of loess shoulder-line land elevation sampling at the Chunhua, Ganquan, and Suide study sites.This figure shows that the three curves of the multi-fractal spectrum are very similar to the curve with a long tail to the right.Determining the factors causing these results is a problem worthy of further discussion with respect to the method of land surface analysis in the Loess gully areas.First, the extraction accuracy of the Loess shoulder line is one of the important factors to influence the study.Many researchers [25,26,34] have performed an in-depth study on the extraction methods and the accuracy of the Loess shoulder line.This study is based on previous researchers' high-accuracy extraction of the Loess shoulder line in the study area and to ensure that Loess shoulder-line data can accurately reflect the topography change of the study site.The Loess shoulder-line land elevation sample data in our study is a one-dimensional elevation series, as in other research fields, such as the time series of climate change, hydrological time series, and seismic data series.There must be noise in the series data, which may be another significant factor to cause the deviation (error) in the result of the topographic analysis of the Loess shoulder line.Thus, noise filtering of the data series is very important when using the MF-DFA method to analyze the topographic changes of the Loess shoulder-line.After the noise removal process described in Section 2.3.1, the Loess shoulder-line land elevation sample series data is analyzed by the MF-DFA method.Figure 11 shows a comparison between the multi-fractal spectrum obtained before and after noise filtering, which differs completely.The results obtained with raw data without any noise removal process, displayed in Figure 11a, show the similarity of the multi-fractal characteristics of the three study sites.Figure 11b shows that the changes in the land surface of the Loess shoulder line are different for the three study sites, and that the topographic change reflected in the multi-fractal spectrum of the Loess shoulder-line is consistent with the actual terrain.A detailed description of the loess shoulder-line topographic variation at the three sites is presented in the Results Section.This is reflected in the land surface analysis of the loess shoulder line, where the noise in the data is a non-negligible impact factor.

Correlation Analysis with Other Geographical Agents
The results of Suide show that the width of the multifractal spectrum of the Suide site is 0.6479, which indicates that the internal variation of the terrain of this site is relatively uniform.This may be related to the land use type of the Suide site, which is mainly artificial terraced fields.The Suide site is characterized by a large number of artificial terraces of various types, which have a great effect on reducing the soil erosion [68].In addition, the site also has a large number of check-dams, which also has a great effect on weakening soil erosion [69][70][71].The check-dam intercepts a large amount of sediment from an upstream channel, and the slope in combination with the siltation of the check-dam, gradually forms a new equilibrium profile upstream, and elevates the erosion datum of its control site.The erosion energy of the upper reaches of the dam and the slope on both sides of the dam gradually decrease; therefore, the erosion effect is weakened, forcing the evolution of the small watershed to accelerate the transition into "old age." The results of the Ganquan site show that the multifractal spectral width is 0.8945, which indicates that the variation in terrain changes along the Ganquan site is very different and the terrain change is more intense.This may also be related to the land use type, rainfall intensity, and so on.The main land use type in the Ganquan site is grassland, but the coverage is very low, and many places are bare loess, with a weak effect on soil erosion.The average annual precipitation at the Ganquan site is ~500 mm, which is mainly concentrated in the form of heavy rain in summer, and soil erosion is very serious.Unlike the Suide site, Ganquan has few such soil and water conservation facilities.Heavy rain in summer could cause serious soil erosion, so the gullies develop very actively, resulting in the intensive terrain variation of the Ganquan site.
The results of the Chunhua site show that the multifractal spectral width is 0.6852, which indicates that the difference in terrain variation along the shoulder line is similar.The landform type of the Chunhua site is dominated by loess tableland, and the ratio of hilly to gully sites is about 8:2.Although the annual precipitation is 600.6 mm, the land use type is mainly garden and woodland, and the vegetation cover is more pronounced.In addition, the loess soil layer at the Chunhua site is thinner and the bedrock is exposed, leading to a high erosion base as well as weak soil erosion.

Conclusions
A multifractal, also known as a multiple fractal measure, reveals a class of complexities and singularity.A whole site with a complex fractal characteristic can be divided into many small sites with different singularities, so that the internal fine structure of the fractal can be understood more thoroughly.
Loess shoulder-line land surface variations series data from the Suide, Ganquan, and Chunhua sites were selected and the MF-DFA method was used to study the variations in these three land sites.The following conclusions were drawn.The loess shoulder-line topographic variations have strong multi-fractal characteristics at the three sites.The strongest multi-fractal characteristics were found at the Ganquan site, followed by the Chunhua site, with the weakest multi-fractal characteristics at the Suide site.The noise and analysis scale of sampling data have a strong influence on the analysis of the loess shoulder line, and the best scale values of the three areas are different, namely 7, 17, and 11 for the Suide, Ganquan, and Chunhua sites, respectively.The topographic spatial variations of the loess shoulder-line are mainly related to other geographic agents, such as lithology, precipitation, and land use.At the Suide site, artificial landforms (such as artificial terraces and check-dams) are the main factor; at the Ganquan site, precipitation; and at the Chunhua site, lithology.The results showed that the combination of EEMD, MF-DFA, and DCCA achieved success in revealing the terrain variations on a watershed.The results could provide theoretical and methodological support for further research on scientific and reasonable classification of loess geomorphic regions.We are optimistic about the potential use of our methods.

Figure 1 .
Figure 1.Data used in this study: (a) Suide; (b) Ganquan; (c) Chunhua; and (d) comparison of the automatic extraction and manual revision (enlargement of part of the Suide site).

Figure 1 .
Figure 1.Data used in this study: (a) Suide; (b) Ganquan; (c) Chunhua; and (d) comparison of the automatic extraction and manual revision (enlargement of part of the Suide site).

Figure 2 .
Figure 2. Multi-fractal analysis of the Loess shoulder line in the Suide site.(a) q-order scale function; (b) q-order Hurst exponent; (c) q-order mass exponent; (d) multifractal spectrum.

Figure 3 .
Figure 3. Multi-fractal analysis of the Loess shoulder line in the Ganquan site.(a) q-order scale function; (b) q-order Hurst exponent; (c) q-order mass exponent; (d) multifractal spectrum.

Figure 2 .
Figure 2. Multi-fractal analysis of the Loess shoulder line in the Suide site.(a) q-order scale function; (b) q-order Hurst exponent; (c) q-order mass exponent; (d) multifractal spectrum.

Figure 2 .
Figure 2. Multi-fractal analysis of the Loess shoulder line in the Suide site.(a) q-order scale function; (b) q-order Hurst exponent; (c) q-order mass exponent; (d) multifractal spectrum.

Figure 3 .
Figure 3. Multi-fractal analysis of the Loess shoulder line in the Ganquan site.(a) q-order scale function; (b) q-order Hurst exponent; (c) q-order mass exponent; (d) multifractal spectrum.

Figure 3 .
Figure 3. Multi-fractal analysis of the Loess shoulder line in the Ganquan site.(a) q-order scale function; (b) q-order Hurst exponent; (c) q-order mass exponent; (d) multifractal spectrum.

Figure 4 .
Figure 4. Multi-fractal analysis of the Loess shoulder line in the Chunhua site.(a) q-order scale function; (b) q-order Hurst exponent; (c) q-order mass exponent; (d) multifractal spectrum.

Figure 4 .
Figure 4. Multi-fractal analysis of the Loess shoulder line in the Chunhua site.(a) q-order scale function; (b) q-order Hurst exponent; (c) q-order mass exponent; (d) multifractal spectrum.

Figure 5 .
Figure 5. Hurst exponent of sampling data of shoulder line with different scales for the Suide site.(a) the standard of profile of shoulder line; (b) Ht variation in analysis scale 7; (c) Ht variation in analysis scale 17.

Figure 5 .
Figure 5. Hurst exponent of sampling data of shoulder line with different scales for the Suide site.(a) the standard of profile of shoulder line; (b) Ht variation in analysis scale 7; (c) Ht variation in analysis scale 17.

Figure 6 .
Figure 6.Hurst exponent of sampling data of shoulder line with different scales for the Ganquan site.(a) the standard of profile of shoulder line; (b) Ht variation in analysis scale 7; (c) Ht variation in analysis scale 17.

Figure 6 .
Figure 6.Hurst exponent of sampling data of shoulder line with different scales for the Ganquan site.(a) the standard of profile of shoulder line; (b) Ht variation in analysis scale 7; (c) Ht variation in analysis scale 17.ISPRS Int.J. Geo-Inf.2017, 6, 141 12 of 19

Figure 7 .
Figure 7. Hurst exponent of sampling data of shoulder line with different scales for the Chunhua site.(a) the standard of profile of shoulder line; (b) Ht variation in analysis scale 7; (c) Ht variation in analysis scale 17.

Figure 7 .
Figure 7. Hurst exponent of sampling data of shoulder line with different scales for the Chunhua site.(a) the standard of profile of shoulder line; (b) Ht variation in analysis scale 7; (c) Ht variation in analysis scale 17.

19 Figure 8 .
Figure 8. DCCA analysis of the local Hurst exponent of the topographic variation and loess shoulder-line elevation sample data at different scales.(a) The α of analysis scale 7 in Suide site; (b) the α of analysis scale 17 in Suide site; (c) the α of analysis scale 7 in Ganquan site; (d) the α of analysis scale 17 in Ganquan site; (e) the α of analysis scale 7 in Chunhua site; (f) the α of analysis scale 17 in Chunhua site.

Figure 8 .
Figure 8. DCCA analysis of the local Hurst exponent of the topographic variation and loess shoulder-line elevation sample data at different scales.(a) The α of analysis scale 7 in Suide site; (b) the α of analysis scale 17 in Suide site; (c) the α of analysis scale 7 in Ganquan site; (d) the α of analysis scale 17 in Ganquan site; (e) the α of analysis scale 7 in Chunhua site; (f) the α of analysis scale 17 in Chunhua site.

Figure 8 .
Figure 8. DCCA analysis of the local Hurst exponent of the topographic variation and loess shoulder-line elevation sample data at different scales.(a) The α of analysis scale 7 in Suide site; (b) the α of analysis scale 17 in Suide site; (c) the α of analysis scale 7 in Ganquan site; (d) the α of analysis scale 17 in Ganquan site; (e) the α of analysis scale 7 in Chunhua site; (f) the α of analysis scale 17 in Chunhua site.

Figure 9 .
Figure 9. Scale effects of DCCA parameters α on the local Hurst exponent of loess shoulder-line and topographic variation.

Figure 9 .
Figure 9. Scale effects of DCCA parameters α on the local Hurst exponent of loess shoulder-line and topographic variation.

19 Figure 10 .
Figure 10.Probability distribution of the local Hurst exponent of loess shoulder-line land surface variation.

Figure 10 .
Figure 10.Probability distribution of the local Hurst exponent of loess shoulder-line land surface variation.

Figure 10 .
Figure 10.Probability distribution of the local Hurst exponent of loess shoulder-line land surface variation.

Figure 11 .
Figure 11.Multi-fractal spectrum of loess shoulder-line elevation sample data: (a) multi-fractal spectrum of raw data; and (b) multi-fractal spectrum of noise filtering data.

Figure 11 .
Figure 11.Multi-fractal spectrum of loess shoulder-line elevation sample data: (a) multi-fractal spectrum of raw data; and (b) multi-fractal spectrum of noise filtering data.

Table 1 .
Geographic description of study sites.
•15 00 -110 • 22 30 E; 37 • 32 30 -37 • 37 30 N Located in Suide county, Shaanxi province with an elevation between 814 and 1188 m, an average slope inclination of 29 • and 9.39 km −1 of gully destiny.The intensive soil erosion in this area characterizes the fragmented land, which have a strong self-like morphology.•3300 -36•35 00 N Located in Ganquan, Shaanxi province, with an elevation between 990 and 1404 m, an average slope inclination of 27 • and 7.74 km −1 of gully density.This area is characterized by an irregular bank-shaped distribution in morphology, which is almost a loess ridge with 200-300 m in width.• 50 00 -34 • 55 00 N Located in Chunhua, Shaanxi province, with an average elevation between 768 and 1188 m, an average slope inclination of 12, and 1.79 km −1 of gully density.The gentle slopes of land surface decreased the total soil loss by erosion.

Table 2 .
Parameters of multi-fractal analysis of the study sites.

Table 2 .
Parameters of multi-fractal analysis of the study sites.