Hyper-temporal remote sensing data in bare soil period and terrain attributes for digital soil mapping in the Black soil regions of China

Catena 184 (2020) 104259

Contents lists available at ScienceDirect

Catena journal homepage: www.elsevier.com/locate/catena

Hyper-temporal remote sensing data in bare soil period and terrain attributes for digital soil mapping in the Black soil regions of China

T

Haoxuan Yanga, Xiaokang Zhangc, Mengyuan Xua, Shuai Shaoa, Xiang Wanga, Wenqi Liua, ⁎ Danqian Wua, Yuyang Maa, Yilin Baoa, Xinle Zhanga, Huanjun Liua,b, a

College of Resources and Environmental Sciences, Northeast Agricultural University, Harbin 150030, China Northeast Institute of Geography and Agroecology, Chinese Academy of Sciences, Changchun 130102, China c International Institute for Earth System Sciences, Jiangsu Provincial Key Laboratory of Geographic Information Science and Technology, Nanjing University, Nanjing 210023, China b

A R T I C LE I N FO

A B S T R A C T

Keywords: Hyper-temporal Digital soil mapping Landsat Soil type Terrain attribute

Remote sensing image data are often used as input in digital soil mapping (DSM). However, it is difficult to distinguish and identify soil types with less difference in reflectance spectral characteristics, because a small amount of input is not enough to provide enough common features. We consider that the hyper-temporal remote sensing data can be used to extract more common features of soil. The accuracy of DSM is improved by using the common features of soil or effective terrain attributes. We took Mingshui County of the Songnen Plain in northeast China as study area, which is known as a Black soil region. STRM DEM, legacy soil data, and 20 scenes Landsat images of bare soil period from 1984 to 2018 (April and May are considered a period of cultivated soil exposure in the study area), were used, with a maximum likelihood method classifier. A digital soil mapping model was constructed based on hyper-temporal data. Results from the study show that the accuracy of mapping with hyper-temporal classification characteristics, with an overall accuracy of 85.18% and a Kappa coefficient of 0.772, is higher than that of mono-temporal classification characteristics, with an average overall accuracy of 64.35% and an average Kappa coefficient of 0.467. After the introduction of relief degree of land surface (RDLS), the overall accuracy and Kappa coefficient of hyper-temporal mapping were 88.22% and 0.818, higher than the accuracy of other terrain factors. The research results signal the advantages of hyper-temporal remote sensing data in DSM, and the common features were able to improve the accuracy of DSM extracted from hyper-temporal data. This paper provided new insight to explain the impact of diverse terrain on DSM of Black soil region, and the mapping of soil type level could be accomplished more easily by the combination of the two characteristics.

1. Introduction

five factors and their derivatives as environmental covariates to complete DSM (Zhang et al., 2017; Minasny and McBratney, 2016). However, the soil landscape model requires soil sampling at a large scale, which is relatively difficult and costly to implement. Therefore, a fast and convenient DSM method is proposed, which with the combination of remote sensing technology and legacy soil data is taken as a covariate (Stumpf et al., 2016; Kempen et al., 2009; Pahlavan-Rad et al., 2016; Boettinger et al., 2008). At present, most legacy soil data are obtained by the traditional field survey method. Although the accuracy of legacy soil data is relatively low in some countries and regions, this still gives some information about the study area. If DSM is accomplished using a combination of RS and legacy soil data, it can reduce costs and improve speed and efficiency; therefore, it is applicable for DSM at a large scale (Sulaeman et al., 2013; Ben-Dor et al., 2008).

The spatial distribution of soil types is an important input for research into precision agriculture and land resource management. At present, the main data source is still the legacy soil data with low accuracy, which is limited by backward mapping technology and fuzzy data. Therefore, it is important to update the soil data with digital soil mapping (DSM) technology (Carré et al., 2007). DSM is a practical technology, and its core principle is the synergistic relationship between soil spatial differences and environmental factors. The spatial changes of soil can be predicted and mapped by the synergistic relationship (McBratney et al., 2003; Zhu et al., 2001; Zhu et al., 2011). As soil is the product of five major soil-forming factors (climate, biology, parent material, terrain, and time), many scholars regard the

⁎

Corresponding author at: College of Resources and Environmental Sciences, Northeast Agricultural University, Harbin 150030, China. E-mail address: [email protected] (H. Liu).

https://doi.org/10.1016/j.catena.2019.104259 Received 23 February 2019; Received in revised form 2 September 2019; Accepted 11 September 2019 0341-8162/ © 2019 Elsevier B.V. All rights reserved.

Catena 184 (2020) 104259

H. Yang, et al.

Fig. 1. Overview of the study area.

temporal remote sensing images. This also means that the hyper-temporal data will reflect more complete soil surface information, which will be conducive to the development of soil classification mapping, as well as ensuring the accuracy of soil classification mapping. Generally, we could extract information from image data and process a large amount of data into one or several datasets by data compression or dimensionality reduction technology (PCA, K-T transformation, etc.). These processed data could represent some characteristics of the study area in a new form (Kunkel et al., 2011; Dobos et al., 2001). Previous studies show that the combination of terrain factors and RS image data could provide more accurate soil type information (Ehsani and Quiel, 2009; Dobos et al., 2000). Terrain is one of the five major soil-forming factors, and is also one of the most commonly used covariates for predicting soil types and properties, especially for DSM (McBratney et al., 2003). Researchers often combine the use of monotemporal or multi-temporal remote sensing data and other data for DSM (Liu et al., 2014; Kang et al., 2008; Boettinger et al., 2008). The soil surface information could be expressed more completely by using hyper-temporal remote sensing data compared with the mono-temporal or multi-temporal remote sensing data. The main objectives of the study are as follows: (1) study the methods of extracting common features of soil types from hyper-temporal remote sensing images for digital soil mapping; (2) find a suitable terrain classification characteristic for DSM and combine it with hyper-temporal classification characteristics to improve the accuracy of soil type mapping in Black soil regions. The results of this study can provide new ideas for DSM by using hypertemporal remote sensing images and a new method for soil types mapping in Black soil regions with similar spectral curves. The study can also be used as technical support for soil protection, monitoring, and sustainable use of cultivated land in Black soil regions.

Mono-temporal remote sensing data, such as vegetation type or terrain, are suitable for extracting static information as one of the most common data types (Metelka et al., 2018; Buis et al., 2009). However, some static information tends to change over time. Therefore, multitemporal data or hyper-temporal data are used in many studies for extracting the characteristics of spatiotemporal dynamic changes, which amplify the difference of data. For example, NDVI changes constantly over time (Wang et al., 2018; Maynard and Levi, 2017), and the parameter extraction and model construction are accomplished by using the change trend of remote sensing image characteristics. By contrast, there are some substances that do not change within a relatively short time period, such as soil type, genus, and trace elements; the differences in these substances, due to errors in temperature, precipitation, season, or image data processing, are hard to distinguish in remote sensing images (Tsui et al., 2004; Liu et al., 2018). The difference between hyper-temporal and multi-temporal can also be interpreted as similar to the difference between hyperspectral and multispectral. Multi-temporal data capture changes mainly at coarse temporal resolutions, and hyper-temporal data are capable of describing long-term variability at a high temporal resolution and capturing the subtle changes in spectral characteristics. The high temporal resolution depends on research objects, for crop growth monitoring or phenological information for example, one or several weeks are usually taken as the temporal resolution of image selection (De Bie et al., 2008; Maynard and Levi, 2017; Khan et al., 2010). For objects that change slowly, however, such as land use change or environmental change, one or several years are taken as the temporal resolution of image selection (Kleynhans et al., 2015; Piwowar and LeDrew, 1995). Hyper-temporal data have more informative detail compared with multi-temporal data. In the existing studies using hyper-temporal remote sensing data, the research objects were mainly residential areas, cultivated land and vegetation, etc. These objects either undergo continuous change or have seasonal rules, which can be reflected by the high temporal resolution characteristics of hyper-temporal data. The variation of characteristics is identified by the differences in hyper-temporal data (Azzari et al., 2019; Ali et al., 2013; De Bie et al., 2008; Kleynhans et al., 2015; Piwowar and LeDrew, 1995). In this paper, we studied the soil type of bare soil using hypertemporal data. Considering that the variation of soil type is slow, it is difficult to reflect the variation of soil type through the temporal differences. It is necessary to extract the common features from hyper-

2. Materials and methods 2.1. Study area The study area is MingShui county (124°18′–125°21′E, 46°44′–47°29′N), which is located in the north of Songnen plain, Heilongjiang province, northeast China (Fig. 1). The county has an area of about 2400 km2. The study area belongs to the monsoon climate of medium latitudes or Dwa in the Köppen-Geiger climate classification, meaning it is warm and rainy in the summer and cold and dry in the 2

Catena 184 (2020) 104259

H. Yang, et al.

winter. The annual average temperature is about 2.9 °C and the annual average precipitation is about 480 mm. The eastern part of the study area is composed of undulating plains. The central region has a higher altitude and the western region has a relatively flat terrain. The region is called a Black soil region by local people because the surface of soil appears black and the topsoil of the region is covered with black or dark humus. The typical soil type in Black soil regions is also directly defined as Black soil in the Genetic Soil Classification of China (Gu et al., 2018); it is also named Phaeozems in the World Reference Base for Soil Resources. The study area belongs to the Black soil region, which includes Chernozems, Meadow soil, and Black soil in the Genetic Soil Classification of China. Chernozems and Meadow soil are called Chernozems and Cambisols, respectively, in WRB. Humus accumulation is significant in Black soil regions; the tillage has high cultivated value with rich organic matter. The area is an important grain production base (Jain et al., 2018). There are some large areas around the world similar to the Black soil region of northeast China, including the Mississippi River basin of the United States, the Great Plains of Ukraine, and the Caco-Pampa region of South America. Agriculture is an important form of land use in these areas (Liu et al., 2017). However, the black soil layer has become thinner and soil erosion intensified with the excessive use of tillage in recent years, so the protection of black soil resources has become an urgent issue (Liu et al., 2010; Zhu et al., 2013; Shen et al., 2019). Most land that could be planted has been turned into cultivated land by farmers in northeast China, and only some land that is prone to waterlogging or salinization has not been cultivated. The study area is dominated by annual vegetation. Most farmers in the study area still use the traditional way of burning to remove crop residue, which hardly has left any residue in the cultivated land. The cultivated land is ploughed from the end of March to the beginning of April every year, and then the soil of the cultivated land is directly exposed to the surface after that (Fig. 2). There was neither a large area of vegetation nor snow cover, as April and May belong to the period of bare soil for cultivated land. The land use types are mainly cultivated land and grassland, and the area ratio of cultivated land and grassland is approximately 3:1 in the study area. Therefore, the bare soil period is used to define this period. We also found that the change rate of land use types did not exceed 5% after comparing land use data in 1984 and 2018. The satellite images capture the surface information, the bare soil is captured in the cultivated land, and some grass and soil are also captured in the grassland. The cultivated land is mainly distributed in the range of Phaeozems, Chernozems, and a small part of Cambisols, Grassland is distributed only in the uncultivated Cambisols. In China, the data on soil types come from the second national soil survey completed in the 1980s. At that time, the technology for soil surveying was backward, which led to the low accuracy of the second soil survey. It is difficult to apply in precision agriculture. In addition, because of the relatively similar spectral reflection curves of disparate soil types in the Black soil region, it is more difficult to use remote

Table 1 Types and dates of remote sensing data in bare soil period for DSM. Image type

Date (mmdd-yyyy)

Image type

Date (mmdd-yyyy)

Image type

Date (mmdd-yyyy)

Landsat-5 Landsat-5 Landsat-5 Landsat-5 Landsat-5 Landsat-5 Landsat-5

05-16-1984 04-17-1985 05-06-1986 04-12-1989 05-17-1990 05-04-1991 04-20-1992

Landsat-5 Landsat-5 Landsat-5 Landsat-5 Landsat-5 Landsat-5 Landsat-5

05-12-1994 05-28-2000 04-29-2001 04-03-2003 05-07-2004 04-08-2005 04-16-2008

Landsat-5 Landsat-5 Landsat-5 Landsat-8 Landsat-8 Landsat-8

05-05-2009 04-22-2010 04-09-2011 04-30-2013 05-03-2014 04-28-2018

sensing images for DSM (Zhang et al., 2018; Liu et al., 2018; Wang et al., 2019). 2.2. The data source Table 1 shows a total of 20 Landsat remote sensing images in bare soil period from 1984 to 2018 in MingShui county (Path 119/Row 27). The Landsat images were acquired from USGS (https://glovis.usgs.gov). The selected images were not covered by clouds. Landsat images of 6 bands were used, including Blue, Green, Red, NIR, SWIR1 and SWIR2. April and May belong to the period of bare soil in the study area (Table 1). So it met the requirements of our study on digital soil mapping. The legacy soil data came from the second national soil survey; the second national soil survey was based on the theory of soil genesis, and soil types were obtained by field surveys and a soil profile analysis (Han and Li, 2018). At the same time, soil evolution and change are a long process, so there was no obvious change in short-term natural and human factors. Disparate soil types have different spectral reflection curves, which is the reason we used remote sensing images for digital soil mapping. Owing to different soil types being developed by different topographic factors, we also used STRM-DEM data, with a spatial resolution of 90 m. 2.3. Data processing method 2.3.1. Data pre-processing ENVI5.3 was used to radiometric calibration and FLAASH atmospheric correction. Since the spatial resolution of Landsat remote sensing image we used was 30 m, we resampled DEM data with a resolution of 90 m, the NEAREST method was used in the resampling method, and the DEM data resolution after resampling was 30 m. Arcgis10.2 was used to clip remote sensing images and DEM data according to the vector boundary of the study area were used for the extraction of subsequent classification characteristics. 2.3.2. The method of classification characteristics extraction (1) Mono-temporal classification characteristics extraction

Fig. 2. Photographs of the soil surface in cultivated land after ploughed. 3

Catena 184 (2020) 104259

H. Yang, et al.

20 scenes of remote sensing images, we tried to construct hyper-temporal classification characteristics by using the three indexes' classification characteristics of 20 images after K–T transformation. An attempt was also made to reduce each exponential feature from twentydimensional band spatial data to three-dimensional band spatial data. Brightness indexes extracted from 20 remote sensing images were stacked, and were reduced dimensionally by PCA. Three principal components were retained as hyper-temporal classification characteristics of the brightness index and written as PC3-brightness. The cumulative contribution rate of PC3-brightness was 68.42%. The greenness and wetness indexes were processed from 20 scene images in the same way, and the hyper-temporal classification characteristics of the greenness index and multi-temporal wetness index were obtained and written as PC3-greenness and PC3-wetness; the cumulative contribution was 76.71% and 63.88%, respectively. We did not use the remaining principal components, because each of them expressed a smaller fraction and contained more redundant information. Hyper-temporal classification characteristics were stacked by the above three indexes classification characteristics and written as PC9.

Kauth–Thomas transformation (K–T) is a special principal component analysis (PCA), which is also named tasseled cap transformation. It can not only reduce the noise of remote sensing data but also enable the processed remote sensing image to achieve the effect of data compression and image enhancement (Zanchetta et al., 2016; Liang, 2008; Luo and Tao, 2017). K–T transformation is performed on 20 remote sensing images, and three indexes are obtained for each remote sensing image, brightness, greenness and wetness. The brightness index is able to reflect the overall reflection effect of surface characteristics, the greenness index reflects the situation of surface vegetation, and the wetness index is able to reflect the water condition of the surface (Serra and Pons, 2008; Chen and Jiang, 2009). A corresponding extension tool in ENVI5.3, called “Tasseled Cap for Landsat8 OLI”, was used to accomplish the K-T transformation for Landsat8 data, and the tool followed the approach of Muhammad Hasan Ali Baig et al. (Baig et al., 2014) to accomplish K-T transformation in Landsat8. We denoted the classification characteristics obtained after K–T transformation of each image as KT3. Remote sensing images of bare soil period could reduce the disturbance of vegetation to soil and describe the reflectance characteristics of soil in cultivated land as much as possible. The classification characteristics represented the brightness, greenness, and wetness of soil after K–T transformation.

(3) Terrain classification characteristics extraction As one of the five major soil-forming factors, terrain is a factor that has to be considered in DSM. The photothermal and moisture conditions of soil were diverse under different terrain attributes (Jenny, 1941). This plays a significant role in the formation of soil, even in the plains. However, it is difficult to directly identify the differences in soil types in the plains from satellite images. Fig. 4 shows the spectral reflection curve characteristics of different soil types and different land use types extracted from the training samples. Three soil types have similar spectral curve characteristics in cultivated land. Therefore, the terrain factor is more important in DSM. Slope (SLO), aspect (ASP), curvature (CUR), relief degree of land surface (RDLS), and elevation (ELE) were considered as five separate topographic classification characteristics in the paper.

(2) Hyper-temporal classification characteristics extraction The spectral characteristic information of each image is a little different, because the shooting time of each image is disparate in the hyper-temporal remote sensing data. The variation of soil type is slow as long as there is no sudden climate change or geological disaster in the short term (Krumbein, 1994; Jenny, 1941). It was difficult to represent the differences in soil types by processing 20 scenes of remote sensing image data into time series data. We further compressed the K-T transformation products from the 20 scenes into a single image to make the processed data smaller and more representative; on the other hand, the method further avoided abnormal situations of partial remote sensing image data (Arias et al., 2018). The detailed process is shown in Fig. 3. In order to accomplish dimensionality reduction and compression of

2.3.3. Soil separability analysis method Jeffries-Matusita distance (J-M distance) was used for analyzing the

Fig. 3. Remote sensing image data processing flow chart. 4

Catena 184 (2020) 104259

H. Yang, et al.

Fig. 4. Soil spectral reflection curve characteristics extracted from Landsat images. (a, Landsat 5 TM with 09 April 2011; b, Landsat 8 OLI with 28 April 2018).

temporal remote sensing data on accuracy of DSM, we built two types of datasets. The mono-temporal remote sensing dataset was based on the classification characteristics of each remote sensing image after K–T transformation (KT3), the hyper-temporal remote sensing dataset was based on the constructed hyper-temporal classification characteristics (PC9), shown in Table 2. Datasets were constructed with different terrain attributes by introducing different terrain classification characteristics. We compared the accuracy to judge the effect of different terrain classification characteristics on mapping in Black soil regions. Table 2 shows the classification characteristics dataset under different terrain attributes. Numbers 2–6 added diverse terrain attributes on the basis of Number 1 in Table 2.

separability of the training samples. The distance values of J-M are in the range 0–2. The closer the distance values of these two soil types are to 2, the higher the separability (Bruzzone and Roli, 1995; Padma and Sanjeevi, 2014). The formulas of J-M distance are listed below.

(J − M distance)i j = 2(1 − e−Bij ) Bij =

|(Vi + Vj )/2| Vi + Vj −1 1 ⎤ (Mi − Mj ) + 1 ln (Mi − Mj )T ⎡ ⎢ 2 ⎥ 8 2 |Vi ||Vj | ⎦ ⎣

(1)

(2)

where (J − M distance)ij is the J-M distance between soil types i and j, Bij is the Bhattacharya distance between soil types i and j, Mi and Mj are the sample mean of soil types i and j, Vi and Vj are the variances of soil types i and j.

2.6. Classification and validation methods 2.4. Training samples and verification samples Machine learning technology has developed rapidly in recent decades, and a large number of models have been made into classifiers and applied into DSM (Lary et al., 2016). Common models include maximum likelihood method, support vector machine, neural network model, etc. Disparate classifiers have different features. For example, Support Vector Machines are suitable for using in a small range of data samples (Jain et al., 2018). Neural Network Model is easy to overfit (Aitkenhead and Coull, 2016). We were aiming for a more convenient and rapid approach of DSM. Due to the large amount of raster data in our study area, we also need to avoid overfitting, so we used the Maximum Likelihood Classifier. This is a common and widely used classifier. It can quickly and conveniently classify most image data, and also effectively avoid overfitting compared with other classifiers. A confusion matrix was used to verify the accuracy of mapping, and the overall accuracy and Kappa coefficient were used to evaluate our mapping results in the confusion matrix. The formulas are listed below.

Both training samples and verification samples were collected from legacy soil data, and the selection of samples was mainly concentrated in typical soil regions (shown in Fig. 1). There were 35 points in Phaeozems, 35 points in Chernozems, and 55 points in Cambisols. The area of Phaeozems was similar to that of Chernozems and the area of Cambisols was a little larger than the above two types in the study area. There were more training sample points in Cambisols. The training sample points were buffer processed outwards in each soil type; the buffer distance was 100 m (Grinand et al., 2008). The final training samples were obtained, including 1763 pixels of Phaeozems, 1698 pixels of Chernozems, and 2635 pixels of Cambisols. The number of verification sample points was double that of the training sample points, so there were 70 points in Phaeozems, 70 points in Chernozems and 110 points in Cambisols. After the same buffer processing, a final verification sample was obtained, including 3368 pixels of Phaeozems, 3357 pixels of Chernozems, and 5320 pixels of Cambisols. The following principles are followed in training sample collection: (1) The same soil type should be sampled separately in disparate regions; because the spectral reflection characteristics of the same soil type are different in disparate regions, some errors could be avoided in mapping. (2) In order to avoid classification errors caused by fuzzy boundary of soil types, training samples should be selected in core areas of soil types as far as possible. (3) The number ratio of sample points of different soil types is roughly the same as the area ratio between the actual soil types in the study area, to ensure that the verification sample distribution is uniform (Qi, 2004).

Table 2 Groups of classification characteristic dataset. Number

Mono-temporal remote sensing data

Hyper-temporal remote sensing data

1 2 3 4 5 6

KT3 KT3 + SLO KT3 + ASP KT3 + CUR KT3 + RDLS KT3 + ELE

PC9 PC9 + SLO PC9 + ASP PC9 + CUR PC9 + RDLS PC9 + ELE

Note: KT3 and PC3 represent mono-temporal and hyper-temporal classification characteristics, respectively; SLO represents slope; ASP represents aspect, CUR represents curvature; RDLS represents relief degree of land surface; ELE represents elevation.

2.5. Classification characteristic datasets In order to compare the effect of mono-temporal and hyper5

Catena 184 (2020) 104259

H. Yang, et al. n

Overall Accuracy =

∑i = 1 Pii N

n

Kappa =

× 100%

the mono-temporal classification characteristic (KT3). Even though the shooting time of each image was diverse, the most common features of all remote sensing images were included in the extracted mono-temporal classification characteristics, which identified the surface characteristics of disparate soil types. As a result, the mapping accuracy of hyper-temporal characteristics was higher than that of mono-temporal characteristics, which indicates that the mono-temporal characteristics were extremely unstable in DSM.

(3)

n

N ∑i = 1 Pii − ∑i = 1 (Pi + × P+i ) N2 −

n ∑i = 1

(Pi + × P+i )

(4)

The column numbers of the confusion matrix are n; Pii is the number of pixels on row i and column i in the confusion matrix, which represents the number of correctly distributed pixels. Pi+ and P+i are the total number of pixels on row i and column i, respectively. N is the total number of pixels used for verification. We evaluate the effect of diverse terrain attributes with the User's accuracy and Producer's accuracy on the mapping accuracy of disparate soil types and the formulas are listed below.

User′s accuracy =

Pii × 100% Pi +

Producer′s Accuracy =

Pii × 100% P+i

3.4. Effects of terrain factors on digital soil mapping Table 3 also lists the accuracy of DSM with different terrain factors as input. It was found that the DSM accuracy of mono-temporal and hyper-temporal characteristics was significantly improved after adding terrain classification characteristics. Compared with the mapping accuracy of mono-temporal characteristics without terrain classification characteristics, the addition of RDLS could increase the overall accuracy by 11.47% and the Kappa coefficient by 15.9% on average; the addition of ELE could increase the overall accuracy by 7.77% and the Kappa coefficient by 11.1% on average. Under the hyper-temporal classification characteristic conditions, the addition of RDLS could increase the overall accuracy by 3.53% and the Kappa coefficient by 5.3%; the addition of ELE could increase the overall accuracy by 3.04% and the Kappa coefficient by 4.6%. The DSM accuracy of Black soil regions could not be significantly improved by the introduction of ASP and CUR. Fig. 6 shows the results of DSM under different classification characteristics. In order to further discuss the effect of terrain classification characteristics on the accurate determination of soil types, we accomplished the mapping and obtained the annual User's accuracy and Producer's accuracy of diverse soil types for each mono-temporal classification characteristics. The annual User's accuracy and Producer's accuracy were averaged for the same soil type and the same terrain attribute. The average User's accuracy and Producer's accuracy of mono-temporal were obtained under diverse soil types and terrain attributes as shown in Table 4, and the average mapping accuracy without terrain (KT3) was subtracted from KT3 plus different terrain attributes. Through the difference value, we were able to directly identify the effect of different terrain attributes on the mapping accuracy of disparate soil types in the Black soil region. As shown in Fig. 7a (Mono-temporal), RDLS improves the mapping accuracy of Phaeozems and Chernozems: the User's accuracy and Producer's accuracy of Phaeozems improved by 25.87% and 21.99% respectively. The User's accuracy and Producer's accuracy of Chernozems improved by 8.10% and 8.06%, respectively. RDLS was able to distinguish between Phaeozems and Chernozems in mapping Black soil regions. ELE was able to effectively improve the mapping accuracy of Cambisols: the User's accuracy and Producer's accuracy of Cambisols were improved by 3.38% and 12.13%, respectively, and the amplitude of improvement was larger than that of other terrain classification characteristics for the mapping accuracy of Cambisols. However, the addition of ASP had little impact on the mapping accuracy of the above three soil types, and the variation in mapping accuracy for the three soil types did not exceed 0.4%. SLO and CUR could also improve the mapping accuracy of Phaeozems and Chernozems, but the amplitude of improvement was smaller than for RDLS. After adding all of the terrain attributes that we extracted, the mapping accuracy was generally higher than that after adding independent terrain attributes. The User's and Producer's accuracy of hyper-temporal was obtained under diverse soil types and terrain attributes. As shown in Table 5, the mapping accuracy without terrain (PC9) was also subtracted from PC9 plus different terrain attributes. As shown in Fig. 7b (hyper-temporal), we found that terrain attributes still improve mapping accuracy, but the amplitude of improvement was smaller than that of using the average of mono-temporal image data. Because the mapping accuracy when using hyper-temporal data without terrain attributes was already very high, it was difficult to improve the mapping accuracy more remarkably after

(5)

(6)

3. Results and analysis 3.1. Spectral characteristic analysis Fig. 1 shows that the land use types are mainly cultivated land and grassland in the study area. The soil types of the cultivated land include Chernozems, Phaeozems and Cambisols, but the soil types of the grassland include only Cambisols. According to Fig. 4 (Cultivated landCambisols and Grassland-Cambisols), there is a certain difference in the spectral reflection curve of Cambisols between cultivated land and grassland during the bared soil period due to the impact of vegetation. Cambisols in cultivated land is basically all exposed to the surface, while Cambisols in grassland is still covered by grass. Since the temperature is generally low in the bare soil period, the grass covering the surface did not turn significantly green. 3.2. Analysis of classification characteristics In order to quantitatively study the separability between soil types under diverse classification characteristics and compare the separability of classification characteristics of mono-temporal and hyper-temporal data, we calculated the J-M distance of training samples in disparate temporal classification characteristics. Fig. 5 shows that the separability of hyper-temporal classification characteristic (PC9) was better than that of mono-temporal classification characteristics (KT3). Mono-temporal classification characteristics lacked stability because they are easily interfered with by other factors, which results in large error. However, the hyper-temporal classification characteristics after compression processing can effectively reduce the interference of image error with the characteristics, and the separability and stability are improved. We also calculated the J-M distance of the training samples under each mono-temporal classification characteristics, which was useful for comparing the classification degree of each terrain classification characteristic for disparate soil types. Fig. 5 shows that RDLS was helpful to distinguish Phaeozems from Chernozems and Phaeozems from Cambisols. ELE had a good effect on distinguishing Cambisols from Chernozems. 3.3. Digital soil mapping accuracy of mono-temporal and hyper-temporal Table 3 shows the digital soil mapping accuracy of mono-temporal and hyper-temporal. The overall accuracy when using hyper-temporal classification characteristics (PC9) was improved by 11.12–30.97%, and the Kappa coefficient was improved by 16.4–44.6% compared with 6

Catena 184 (2020) 104259

H. Yang, et al.

Fig. 5. J-M distance between soil types under disparate terrain attributes. (a, Cambisols-Chernozems; b, Cambisols-Phaeozems; c, Phaeozems-Chernozems).

uses the advantages of hyper-temporal data for common feature extraction. Because the research objects have slow transformation (soil types), the advantages of high temporal resolution are difficult to pinpoint. From the verification results, we found that the OA and Kappa coefficient of using hyper-temporal data without terrain factors could reach 85.18% and 0.772, respectively; the average OA and Kappa coefficient of using mono-temporal data without terrain factors were 64.35% and 0.467, respectively. Therefore, the accuracy of mapping with hyper-temporal data could be greatly improved, and the stability of mapping could be guaranteed.

the introduction of terrain attributes. In addition, the mapping accuracy of the three soil types studied in this paper could be improved by using hyper-temporal data. Generally, terrain attributes were conducive to the improvement of soil mapping accuracy, but the sensitivity of diverse terrain factors to disparate soil types was different. 4. Discussion 4.1. Advantages of DSM based on hyper-temporal remote sensing data

4.2. The differences and causes of mapping accuracy between each scene of mono-temporal remote sensing data

In this paper, 20 scenes of remote sensing images are used in DSM. The results show that the accuracy and stability of models using hypertemporal remote sensing images were much higher than models using mono-temporal remote sensing images (Ma et al., 2017; Massawe et al., 2018). The OA and Kappa coefficient using different mono-temporal images' data are still unstable. In general, using hyper-temporal remote sensing images has advantages in terms of quickly obtaining data needed for DSM mapping. Moreover, the mapping process is simple and low-cost. Previous studies used hyper-temporal data to study the surface objects, which included the detection of settlement expansion (Kleynhans et al., 2015), monitoring of environmental change and land cover (Piwowar et al., 1998; Cable et al., 2014), crop statistics (Khan et al., 2010), carbon storage (Ren et al., 2013), etc. Their research mainly used the advantage of time series analysis of hyper-temporal data in high temporal resolution to accomplish the research. Our study mainly

The accuracies of DSM were different when using diverse monotemporal classification characteristics, which indicate that the soil surface reflectance, expressed by each Landsat satellite images, is also different. On one hand, due to the different state of the surface at different shooting times, the image of remote sensing would be affected by some factors such as precipitation and temperature, which would cause further differences in the surface reflectance. On the other hand, the interannual vegetation reflectance was also different, especially for Cambisols covered with grass, which would influence the mapping accuracy. The surface coverings have a certain impact on the mapping accuracy at the pixel scale. Because those data are used to accomplish the mapping, the overall accuracy is low. Moreover, the OA and Kappa 7

Catena 184 (2020) 104259

H. Yang, et al.

Table 3 Overall accuracy and Kappa coefficient of DSM by using different classification characteristic dataset. Temporal type

Mono-temporal

Hyper-temporal

Date (mm-dd-yyyy)

05-16-1984 04-17-1985 05-06-1986 04-12-1989 05-17-1990 05-04-1991 04-20-1992 05-12-1994 05-28-2000 04-29-2001 04-03-2003 05-07-2004 04-08-2005 04-16-2008 05-05-2009 04-22-2010 04-09-2011 04-30-2013 05-03-2014 04-28-2018 max avg min

KT3

KT3 + SLO

KT3 + ASP

KT3 + CUR

KT3 + RDLS

KT3 + ELE

KT3 + ALL

OA

OA

OA

OA

OA

OA

Kappa

OA

68.40% 63.12% 67.23% 66.10% 74.06% 66.12% 69.77% 69.25% 66.75% 58.86% 54.21% 70.16% 58.22% 64.33% 62.21% 55.08% 69.22% 57.78% 61.17% 65.05% 74.06% 64.35% 54.21% PC9 85.18%

0.518 0.446 0.505 0.489 0.608 0.489 0.548 0.536 0.502 0.383 0.326 0.550 0.382 0.465 0.433 0.354 0.540 0.371 0.415 0.470 0.608 0.467 0.326

74.63% 0.612 70.44% 0.555 72.38% 0.581 71.42% 0.568 75.20% 0.625 71.40% 0.567 72.09% 0.581 72.45% 0.582 70.30% 0.552 66.98% 0.501 60.93% 0.421 72.11% 0.578 63.77% 0.459 69.38% 0.536 68.78% 0.526 57.21% 0.382 71.41% 0.570 61.02% 0.415 66.95% 0.500 71.24% 0.562 75.20% 0.625 69.00% 0.534 57.21% 0.382 PC9 + SLO 86.09% 0.785

0.772

Kappa

Kappa

68.53% 0.520 63.14% 0.446 66.95% 0.501 66.33% 0.492 73.55% 0.600 66.34% 0.493 69.64% 0.546 69.27% 0.536 66.64% 0.500 58.85% 0.382 54.33% 0.328 70.12% 0.550 58.66% 0.388 64.30% 0.465 62.32% 0.435 55.10% 0.354 69.03% 0.537 57.45% 0.366 61.37% 0.418 64.80% 0.467 73.55% 0.600 64.34% 0.466 54.33% 0.328 PC9 + ASP 85.24% 0.773

Kappa

71.61% 0.567 67.13% 0.506 69.88% 0.544 69.20% 0.535 74.22% 0.610 68.97% 0.531 70.76% 0.562 70.53% 0.554 68.52% 0.526 63.01% 0.444 57.85% 0.378 71.24% 0.566 61.16% 0.423 66.77% 0.499 65.55% 0.480 55.78% 0.363 69.78% 0.547 57.87% 0.371 64.23% 0.460 68.41% 0.520 74.22% 0.610 66.62% 0.499 55.78% 0.363 PC9 + CUR 85.52% 0.777

Kappa

81.03% 0.710 77.62% 0.662 78.34% 0.672 77.26% 0.656 79.14% 0.684 78.37% 0.672 76.37% 0.645 79.35% 0.686 77.39% 0.658 76.01% 0.637 69.82% 0.549 76.59% 0.646 71.27% 0.570 76.55% 0.644 76.71% 0.645 63.20% 0.467 77.25% 0.656 70.00% 0.548 74.84% 0.619 79.29% 0.683 81.03% 0.710 75.82% 0.635 63.20% 0.467 PC9 + RDLS 88.22% 0.818

Kappa

73.66% 0.597 71.07% 0.562 75.04% 0.620 73.33% 0.595 80.51% 0.703 73.43% 0.595 75.87% 0.636 75.92% 0.633 73.11% 0.592 69.13% 0.531 62.72% 0.444 78.08% 0.665 68.34% 0.524 72.92% 0.588 71.72% 0.568 63.42% 0.464 76.94% 0.650 67.15% 0.504 69.00% 0.529 71.09% 0.559 80.51% 0.703 72.12% 0.578 62.72% 0.444 PC9 + ELE 87.71% 0.810

Kappa

82.91% 0.737 79.10% 0.683 81.07% 0.711 79.63% 0.690 82.32% 0.729 81.06% 0.710 78.69% 0.676 82.08% 0.725 80.66% 0.705 78.70% 0.675 74.99% 0.622 81.95% 0.723 76.20% 0.639 78.90% 0.677 80.45% 0.699 68.18% 0.530 81.42% 0.715 75.13% 0.621 76.78% 0.646 79.93% 0.692 82.91% 0.737 79.01% 0.680 68.18% 0.530 PC9 + ALL 89.95% 0.844

Note: Max, avg, and min are the maximum, mean, and minimum values of all mono-temporal image classification characteristics. OA is the overall accuracy and kappa is the Kappa coefficient.

drought in 40 years in 2003, accompanied by large-scale sand and dust weather, which affected the surface reflectance (Choobari et al., 2014; Chen et al., 2018). Irregular variations of surface reflectance reduce the mapping accuracy. By contrast, the image of 17 May 1990 could better reflect the situation of bare soil, because there was no abnormal weather when the image was taken. These are shown in Fig. 8(a and b).

coefficient of mapping with 17 May 1990 remote sensing images are 74.06% and 0.608, respectively, and the OA and Kappa coefficient of mapping with 3 April 2003 remote sensing images are 54.21% and 0.326, respectively; the OA and Kappa coefficient differences between the two images could reach 19.84% and 28.20%. respectively. The main reason is that Heilongjiang Province experienced the worst spring

Fig. 6. Digital soil mapping of different temporal data under diverse terrain classification attributes. 8

Catena 184 (2020) 104259

H. Yang, et al.

Table 4 Average mapping accuracy of soil types by using mono-temporal classification characteristic dataset. TYPE (monotemporal)

Mapping accuracy (average)

Phaeozems

Meadows

Chernozems

KT3

prod.acc user.acc prod.acc user.acc prod.acc user.acc prod.acc user.acc prod.acc user.acc prod.acc user.acc prod.acc user.acc

58.56% 53.82% 65.12% 66.18% 58.16% 53.89% 62.50% 59.75% 80.54% 79.69% 65.86% 66.09% 80.17% 81.20%

64.73% 86.68% 69.09% 86.54% 64.86% 86.57% 66.82% 86.70% 71.69% 87.45% 76.86% 90.06% 81.06% 89.34%

69.57% 53.30% 72.77% 54.91% 69.71% 53.20% 70.45% 53.89% 77.63% 61.40% 70.90% 57.99% 74.60% 64.79%

KT3 + SLO KT3 + ASP KT3 + CUR KT3 + RDLS KT3 + ELE KT3 + ALL

Table 5 Mapping accuracy of soil types by using hyper-temporal classification characteristic dataset. TYPE (hypertemporal)

Mapping accuracy

Phaeozems

Meadows

Chernozems

PC9

prod.acc user.acc prod.acc user.acc prod.acc user.acc prod.acc user.acc prod.acc user.acc prod.acc user.acc prod.acc user.acc

87.44% 80.95% 89.31% 83.12% 87.38% 81.19% 88.45% 81.48% 93.32% 87.45% 89.99% 87.15% 93.50% 90.20%

88.36% 88.38% 89.81% 88.70% 88.57% 88.27% 88.98% 88.50% 90.51% 88.95% 94.36% 90.30% 96.05% 90.11%

77.87% 84.65% 76.94% 85.00% 77.81% 84.75% 77.09% 85.13% 79.48% 87.82% 78.46% 87.57% 79.71% 89.32%

PC9 + SLO PC9 + ASP PC9 + CUR PC9 + RDLS PC9 + ELE PC9 + ALL

Note: prod.acc represents User's accuracy and user.acc represents Producer's accuracy.

Note: prod.acc represents User's accuracy and user.acc represents Producer's accuracy.

mono-temporal data for improving the effects of DSM: the accuracy of mapping was improved, and the stability of mapping was guaranteed. The main reason was that the hyper-temporal classification characteristics extracted in this paper could remove the redundant information from the image, and the common features formed with the compression and dimensionality reduction of multiple remote image data could effectively reveal the surface information of disparate soil types, even some vegetation can also show partial soil information, especially for the Cambisols covered with grass. The method of extracting common features was easier for transforming remote sensing data into research inputs. The layer dimension of the input was controlled based on demand, which promoted the integration of remote sensing technology and DSM. At present, remote sensing images have been used more widely, so the method in this paper has the potential to be used to extract or monitor surface information with slow-changing features. In

Therefore, compared with mono-temporal data, hyper-temporal data could provide common soil characteristics for DSM, without considering random factors such as interannual differences in imaging conditions and surface soil moisture. At the same time, it would be more stable to use the classification characteristics of hyper-temporal data because small errors in some image pixels would not lead to errors of mapping. As shown in Fig. 9, there was still some instability in the mapping of mono-temporal classification characteristics, even after the introduction of the terrain factor. 4.3. The significance and insufficiency of the study In this research we showed that hyper-temporal data are better than

Fig. 7. The soil types mapping accuracy difference under disparate terrain attributes. 9

Catena 184 (2020) 104259

H. Yang, et al.

Fig. 8. Landsat images of bare soil in different years.

techniques in the laboratory reached the level of soil genus within the Chinese Soil Genetic Classification System (e.g., Sticky calcareous meadow soils, Loess calcareous chernozem, Loess black soils et al.), with an OA of 83.3% and a Kappa coefficient of 0.8 (Wang et al., 2019), we will try to use hyper-temporal remote sensing data to DSM at the soil genus level in black soil regions. Because only Cambisols is covered by grass in the study area, Cambisols characteristics that we extracted actually contain some vegetation information. Regardless, it was still enough to accomplish the mapping. In the future, we will try to use the approach of Demattê et al. (2018) to deal with the vegetation issues on the soil surface.

addition, RDLS and ELE were effective terrain factors to improve the mapping accuracy of DSM in Black soil regions, and were also conducive to the development of research on the evolution and variation of soil types. There are still some defects in this study: since the accuracy of mapping is determined by the selection of classification characteristics, it remains to be verified whether the classification features selected in this study can be applied successfully in other study areas. Moreover, there are some differences of terrain attributes between the Black soil region in northeast China and other countries; the universality of the DSM method combined with terrain features still needs to be proven. Moreover, the common features of hyper-temporal remote sensing data might ignore some important information in the extraction process. These aspects still need attention and exploration.

5. Conclusions This study demonstrated the effectiveness of using hyper-temporal remote sensing data and terrain information to improve the mapping accuracy of soil types. We also emphasized another advantage of hypertemporal remote sensing data: they are conducive to extracting common feature information, especially for features similar to soil types that cannot change in the short term. The common features of surface soil information were extracted from a mass of image data. We accomplished the DSM of soil types by using this approach. This avoided the adverse impact of errors in mono-temporal remote sensing data and improved the stability of DSM by using remote sensing images. In addition, the importance of terrain was confirmed: RDLS could improve the mapping accuracy of Phaeozems and Chernozems, and ELE could help improve the mapping accuracy of Cambisols. In this paper, the hyper-temporal remote sensing data in the bare soil period and RDLS were used in the mapping of soil types in the typical Black soil region of Songnen Plain, with an OA of 88.22% and a kappa coefficient of 0.818. In addition to the three types of soil mapped in this paper, there are other types of soil in the Black soil region of Songnen plain in China,

4.4. Future research perspectives on the combination of DSM and hyperspectral technologies In the laboratory, a spectrometer is used to classify soil samples within the wavelength range of 500–2500 nm, with the highest accuracy reaching 97% (Zhang et al., 2018). The soil samples, which are collected in the field, will be used as training samples or verification samples before classifying or allocating by hyperspectral technology in the future. In order to enrich the acquisition methods and quantities of training samples or verification samples, we could collect spectral data of soil samples and predict the soil type based on the models, and then using that modeled point data to increase training or verification data in the hyper-temporal DSM model. With the rapid application of hyperspectral satellite remote sensing data in recent years, we can use hyperspectral satellite remote sensing data to collect more accurate data on soil types in the next study. Although the soil classification or allocation by hyperspectral

Fig. 9. The Overall Accuracy and Kappa coefficient of all mono-temporal remote sensing data (a, Overall Accuracy; b, Kappa coefficient). 10

Catena 184 (2020) 104259

H. Yang, et al.

including boggy soil and aeolian soil. It is also necessary to discuss whether DSM can be accomplished in other areas in the future by combining hyper-temporal remote sensing data or hyperspectral satellite image data with terrain attributes. Moreover, we will attempt to implement mapping of soil type or genus by using hyperspectral data due to the increasing popularity of hyperspectral satellite remote sensing images.

Han, X., Li, N., 2018. Research progress of black soil in northeast China. Sci. Geogr. Sin. 38, 1032–1041. Jain, D.K., Dubey, S.B., Choubey, R.K., Sinhal, A., Arjaria, S.K., Jain, A., Wang, H., 2018. An approach for hyperspectral image classification by optimizing SVM using self organizing map. Journal of Computational Science 25, 252–259. Jenny, H., 1941. Factors of soil formation: a system of quantitative pedology/Hans Jenny. Soil Sci. 42, 415. Kang, Q., Zhang, Z.X., Zhao, X.L., 2008. A study of soil classification based on remote sensing in arid area. Journal of Remote Sensing 12, 159–167. Kempen, B., Brus, D.J., Heuvelink, G.B.M., Stoorvogel, J.J., 2009. Updating the 1:50,000 Dutch soil map using legacy soil data: a multinomial logistic regression approach. Geoderma 151, 311–326. Khan, M.R., de Bie, C.A.J.M., van Keulen, H., Smaling, E.M.A., Real, R., 2010. Disaggregating and mapping crop statistics using hypertemporal remote sensing. Int. J. Appl. Earth Obs. Geoinf. 12, 36–46. Kleynhans, W., Salmon, B.P., Olivier, J.C., 2015. Detecting settlement expansion in South Africa using a hyper-temporal SAR change detection approach. Int. J. Appl. Earth Obs. Geoinf. 42, 142–149. Krumbein, W.C., 1994. Factors of soil formation: a system of quantitative pedology by Hans Jenny. Geoderma 68, 336–337. Kunkel, Melvin L., Flores, Alejandro N., Smith, Toni J., McNamara, James P., Benner, Shawn G., 2011. A simplified approach for estimating soil carbon and nitrogen stocks in semi-arid complex terrain. Geoderma 165, 1–11. Lary, D.J., Alavi, A.H., Gandomi, A.H., Walker, A.L., 2016. Machine learning in geosciences and remote sensing. Geosci. Front. 7, 3–10. Liang, Zhang, 2008. Study on the hyperspectral remote sensed image classify based on PCA and SVM. In: Optical Technique. Liu, X.B., Zhang, X.Y., Wang, Y.X., Sui, Y.Y., Zhang, S.L., Herbert, S.J., Ding, G., 2010. Soil degradation: a problem threatening the sustainable development of agriculture in Northeast China. Plant Soil Environment 56, 87–97. Liu, Juan, Cai, Yanjun, Wang, Jin, 2014. Soil classification of Qinghai Lake basin based on remote sensing. Remote Sensing for Land Resources 26, 57–62. Liu, Xiaobing, Charles Lee Burras, Yuri S. Kravchenko, Artigas Duran, Ted Huffman, Hector Morras, Guillermo Studdert, Xingyi Zhang, Richard M. Cruse, and Xiaohui %J Revue Canadienne De La Science Du Sol Yuan. 2017. 'Overview of Mollisols in the World: Distribution, Land Use and Management: Revue Canadienne de Phytotechnie', 92: 383–402. Liu, H., Yang, H., Xu, M., Zhang, X., Zhang, X., Yu, Z., Shao, S., Li, H., 2018. Soil classification based on maximum likelihood method and features of multi-temporal Landsat 8 remote sensing images in bare soil period. Transactions of the Chinese Society of Agricultural Engineering 34, 132–139. Luo, Kaisheng, Tao, Fulu, 2017. Method for wetland type extraction using remote sensing combing object-oriented and tasseled cap transformation. Transactions of the Chinese Society of Agricultural Engineering 33, 198–203. Ma, Yuxin, Minasny, Budiman, Wu, Chunfa, 2017. Modelling and mapping of key soil properties to support agricultural production in Eastern China. Geoderma Regional 10. Massawe, B.H.J., Subburayalu, S.K., Kaaya, A.K., Winowiecki, L., Slater, B.K., 2018. Mapping numerically classified soil taxa in Kilombero Valley, Tanzania using machine learning. Geoderma 311, 143–148. Maynard, J.J., Levi, M.R., 2017. Hyper-temporal remote sensing for digital soil mapping: characterizing soil-vegetation response to climatic variability. Geoderma 285, 94–109. McBratney, A.B., Mendonça Santos, M.L., Minasny, B., 2003. On digital soil mapping. Geoderma 3–52. Metelka, V., Baratoux, L., Jessell, M.W., Barth, A., Ježek, J., Naba, S., 2018. Automated regolith landform mapping using airborne geophysics and remote sensing data, Burkina Faso, West Africa. Remote Sens. Environ. 204, 964–978. Minasny, B., McBratney, A.B., 2016. Digital soil mapping: a brief history and some lessons. Geoderma 264, 301–311. Padma, S., Sanjeevi, S., 2014. Jeffries Matusita based mixed-measure for improved spectral matching in hyperspectral image analysis. Int. J. Appl. Earth Obs. Geoinf. 32, 138–151. Pahlavan-Rad, M.R., Khormali, F., Toomanian, N., Brungard, C.W., Kiani, F., Komaki, C.B., Bogaert, P., 2016. Legacy soil maps as a covariate in digital soil mapping: a case study from Northern Iran. Geoderma 279, 141–148. Piwowar, J.M., LeDrew, E.F., 1995. Hypertemporal analysis of remotely sensed sea-ice data for climate change studies. Prog. Phys. Geogr. 19, 216–242. Piwowar, J.M., Peddle, D.R., LeDrew, E.F., 1998. Temporal mixture analysis of Arctic sea ice imagery: a new approach for monitoring environmental change. Remote Sens. Environ. 63, 195–207. Qi, F., 2004. Knowledge discovery from area-class resource maps: data preprocessing for noise reduction. Trans. GIS 8, 297–308. Ren, H., Chen, H., Li, L., Li, P., Hou, C., Wan, H., Zhang, Q., Zhang, P., 2013. Spatial and temporal patterns of carbon storage from 1992 to 2002 in forest ecosystems in Guangdong, Southern China. J Plant, and Soil 363, 123–138. Serra, P., Pons, X., 2008. Monitoring farmers’ decisions on Mediterranean irrigated crops using satellite image time series. Int. J. Remote Sens. 29, 2293–2316. Shen, Qingsong, Wang, Yao, Wang, Xinrui, Liu, Xu, Zhang, Xingyi, Zhang, Shaoliang, 2019. Comparing interpolation methods to predict soil total phosphorus in the Mollisol area of Northeast China. CATENA 174, 59–72. Stumpf, Felix, Schmidt, Karsten, Behrens, Thorsten, Schönbrodt-Stitt, Sarah, Buzzo, Giovanni, Dumperth, Christian, Wadoux, Alexandre, Xiang, Wei, 2016. Incorporating limited field operability and legacy soil samples in a hypercube sampling design for digital soil mapping. Thomas %J Journal of Plant Nutrition Scholten, and Soil Science = Zeitschrift fuer Pflanzenernaehrung und Bodenkunde 179, 499–509.

Acknowledgments This research was supported by National Natural Science Foundation of China (41671438) and Natural Science Foundation of Heilongjiang Province of China (D2017001). We thank LetPub (www. letpub.com) for its linguistic assistance during the preparation of this manuscript. References Aitkenhead, M.J., Coull, M.C., 2016. Mapping soil carbon stocks across Scotland using a neural network model. Geoderma 262, 187–198. Ali, A., de Bie, C.A.J.M., Skidmore, A.K., Scarrott, R.G., Hamad, A., Venus, V., Lymberakis, P., 2013. Mapping land cover gradients through analysis of hyper-temporal NDVI imagery. Int. J. Appl. Earth Obs. Geoinf. 23, 301–312. Arias, O.V., Garrido, A., Villeta, M., Tarquis, A.M., 2018. Homogenisation of a soil properties map by principal component analysis to define index agricultural insurance policies. Geoderma 311, 149–158. Azzari, G., Grassini, P., Edreira, J.I.R., Conley, S., Mourtzinis, S., Lobell, D.B., 2019. Satellite mapping of tillage practices in the North Central US region from 2005 to 2016. Remote Sens. Environ. 221, 417–429. Baig, M.H.A., Zhang, L., Shuai, T., Tong, Q., 2014. Derivation of a tasselled cap transformation based on Landsat 8 at-satellite reflectance. Remote Sensing Letters 5, 423–431. Ben-Dor, E., Taylor, R.G., Hill, J., Demattê, J.A.M., Whiting, M.L., Chabrillat, S., Sommer, S., 2008. Imaging Spectrometry for Soil Applications. Advances in Agronomy Academic Press. Boettinger, J.L., Ramsey, R.D., Bodily, J.M., Cole, N.J., Kienast-Brown, S., Nield, S.J., Saunders, A.M., Stum, A.K., 2008. Landsat spectral data for digital soil mapping. In: Digital Soil Mapping with Limited Data, pp. 193–202. Bruzzone, L., Roli, F., 1995. An extension of the Jeffreys-Matusita distance to multiclass cases for feature selection. IEEE Trans. Geosci. Remote Sens. 33, 1318–1321 %@ 0196-2892. Buis, E., Veldkamp, A., Boeken, B., van Breemen, N., 2009. Controls on plant functional surface cover types along a precipitation gradient in the Negev Desert of Israel. J. Arid Environ. 73, 82–90. Cable, Jeffrey W., Kovacs, John M., Shang, Jiali, Jiao, Xianfeng, 2014. Multi-temporal Polarimetric RADARSAT-2 for Land Cover Monitoring in Northeastern Ontario, Canada. 6. pp. 2372. Carré, F., McBratney, A.B., Mayr, T., Montanarella, L., 2007. Digital soil assessments: beyond DSM. Geoderma 142, 69–79. Chen, Chao, Jiang, Tao, 2009. Research on remote sensing image fusion methods based on tasseled cap transformation. Science of Surveying Mapping 4, S155–S156. Chen, S., Yuan, T., Zhang, X., Zhang, G., Feng, T., Zhao, D., Zang, Z., Liao, S., Ma, X., Jiang, N., Zhang, J., Yang, F., Lu, H., 2018. Dust modeling over East Asia during the summer of 2010 using the WRF-Chem model. J. Quant. Spectrosc. Radiat. Transf. 213, 1–12. Choobari, O.A., Zawar-Reza, P., Sturman, A., 2014. The global distribution of mineral dust and its impacts on the climate system: a review. Atmos. Res. 138, 152–165. De Bie, C.A.J.M., Khan, M.R., Toxopeus, A.G., Venus, V., Skidmore, A.K., 2008. Hypertemporal image analysis for crop mapping and change detection. In: ISPRS 2008. Proceedings of the XXI Congress: Silk Road for Information From Imagery: The International Society for Photogrammetry and Remote Sensing, Beijing, China, pp. 803–812. Demattê, J.A.M., Fongaro, C.T., Rizzo, R., Safanelli, J.L., 2018. Geospatial Soil Sensing System (GEOS3): a powerful data mining procedure to retrieve soil spectral reflectance from satellite images. Remote Sens. Environ. 212, 161–175. Dobos, E., Micheli, E., Baumgardner, M.F., Biehl, L., Helt, T., 2000. Use of combined digital elevation model and satellite radiometric data for regional soil mapping. Geoderma 97, 367–391. Dobos, E., Montanarella, L., Nègre, T., Micheli, E., 2001. A regional scale soil mapping approach using integrated AVHRR and DEM data. Int. J. Appl. Earth Obs. Geoinf. 3, 30–42. Ehsani, A.H., Quiel, F., 2009. A semi-automatic method for analysis of landscape elements using Shuttle Radar Topography Mission and Landsat ETM+ data. Comput. Geosci. 35, 373–389. Grinand, C., Arrouays, D., Martin, M.P., Laroche, B., 2008. Extrapolating regional soil landscapes from an existing soil map: sampling intensity, validation procedures, and integration of spatial context. Geoderma 143, 180–190. Gu, Z., Xie, Y., Gao, Y., Ren, X., Cheng, C., Wang, S., 2018. Quantitative assessment of soil productivity and predicted impacts of water erosion in the black soil region of northeastern China. Sci. Total Environ. 637–638, 706–716.

11

Catena 184 (2020) 104259

H. Yang, et al.

1–15. Zhang, G., Liu, F., Song, X., 2017. Recent progress and future prospect of digital soil mapping: a review. J. Integr. Agric. 16, 2871–2885. Zhang, Xiaokang, Liu, Huanjun, Zhang, Xinle, Yu, Shengnan, Dou, Xin, Xie, Yahui, Wang, Nan, 2018. Allocate soil individuals to soil classes with topsoil spectral characteristics and decision trees. Geoderma 320, 12–22. Zhu, A.X., Hudson, B., Burt, J., Lubich, K., Simonson, D., 2001. Soil mapping using GIS, expert knowledge, and fuzzy logic. Soil Sci. Soc. Am. J. 65, 1463–1472. Zhu, A Xing, James E. Burt, Amanda C. Moore, and Michael P. Smith. 2011. ‘SoLIM: A New Technology for Soil Mapping Using GIS, Expert Knowledge & Fuzzy Logic’. Zhu, X., Han, B., Zhao, J., 2013. Spatial and temporal variability of soil nutrients in the black soil area of northeast China. Journal of Food Agriculture Environment 11, 1386–1389.

Sulaeman, Yiyi, Minasny, Budiman, McBratney, Alex B., Sarwani, Muhrizal, Sutandi, Atang, 2013. Harmonizing legacy soil data for digital soil mapping in Indonesia. Geoderma 192, 77–85. Tsui, C.-C., Chen, Z.-S., Hsieh, C.-F., 2004. Relationships between soil properties and slope position in a lowland rain forest of southern Taiwan. Geoderma 123, 131–142. Wang, J., Xiao, X., Qin, Y., Doughty, R.B., Dong, J., Zou, Z., 2018. Characterizing the encroachment of juniper forests into sub-humid and semi-arid prairies from 1984 to 2010 using PALSAR and Landsat data. Remote Sens. Environ. 205, 166–179. Wang, X., Zhang, X., Li, H., Zhang, X., Liu, H., Dou, X., Yu, Z., 2019. The minimum level for soil allocation using topsoil reflectance spectra: genus or species? CATENA 174, 36–47. Zanchetta, A., Bitelli, G., Karnieli, A., 2016. Monitoring desertification by remote sensing using the Tasselled Cap transform for long-term change detection. Nat. Hazards 83,

12