The Suitability of Remote Sensing Images at Different Resolutions for Mapping of Gullies in the Black Soil Region, Northeast China

Remote sensing images with different spatial resolutions have different performance capabilities for gully extraction, so it is very important to study the suitability of different spatial resolutions for this purpose. In this study, part of the black soil area in Northeast China with serious gully erosion was taken as the study area, and Google Earth images with seven spatial resolutions ranging from 0.51 to 32.64 m, commonly used in gully erosion research, were selected as data sources. Combined with auxiliary data, gullies were extracted by visual interpretation. The interpretation results of images of different spatial resolutions were analyzed qualitatively and quantitatively, and the interpretation suitability of images of different spatial resolutions for different types of gullies under different classification systems was emphatically explored. The results indicate that the image with a spatial resolution of 1.02 m has the best performance when not considering the types of gullies. However, the image with a spatial resolution of 2.04 m is the most cost-effective and, therefore, the most suitable for general research. When it is necessary to distinguish the type of gully, the image with a spatial resolution of 0.51 m can be adapted for all situations. However, research on ephemeral gullies is of little practical significance. Therefore, the image with a spatial resolution of 1.02 m is the most universally useful image, being cheaper and easier to obtain. When the spatial resolution is 2.04 m or lower, it is necessary to select the spatial resolution according to the gully type required for practical application. When the spatial resolution is 8.16 or lower, the interpretation of gullies becomes very difficult or even impossible.


Introduction
Gully erosion refers to the vertical erosion, lateral erosion and retrogressive erosion of gully flow, resulting in the formation of gullies [1]. The phenomenon is most obvious in areas with heavy rainfall, slopes, sparse vegetation and thick loose material coverage. Despite being a natural process, excessive human activities such as long-term cultivation and overgrazing are the primary inducing factors of gully erosion [2,3]. Several studies proposed that gully erosion is the main soil erosion type in the black soil region of Northeast China, causing serious soil erosion that erodes cultivated land, destroys roads and causes great soil and water loss [4,5]. Gully mapping is very useful for land use planning, conservation practices and soil erosion mitigation [6][7][8].
With the development of remote sensing techniques, various remote sensing data have become cheap and readily available. Compared with field investigation using tools such as micro-topographic profile meters [9], total stations [10], erosion pins [11] or global navigation satellite systems [12] and aerial photo interpretation by photogrammetric techniques [13,14], mapping gullies using satellite imagery interpretation has been considered to be a preferable method with the characteristics of lower labor and time consumption, availability to larger regions of interest and dynamic research. Thus, satellite imagery interpretation has been widely used [15][16][17].
Visual image interpretation, as the earliest way of mapping gullies by satellite imagery interpretation [18], has been applied for years due to its high precision and stability [8,19]. However, inefficiency and the requirement for background knowledge of the interpreter limit its application. In view of these flaws, automatic gully extraction methods have been developed [20][21][22], including pixel-based and object-oriented image analysis approaches. Compared with pixel-based extraction, object-oriented extraction takes the spatial and geometric information of an object into consideration, which typically makes it outperform pixel-based extraction [23,24]. However, object segmentation is not just a key to objectoriented image analysis-it is a very difficult process [25,26]. For any gully extraction method involving satellite imagery interpretation, selecting a suitable spatial resolution of remote sensing data is an initial and significant step [27]. Thus, Luoman et al. [28] compared the results of several remote sensing images with different resolutions when extracting gullies by visual interpretation. They reported that remote sensing images with a higher spatial resolution did not always provide the best information as they contained more noise. Some scholars have studied the selection of a suitable pixel size when using automatic gully extraction methods. Using five datasets with spatial resolutions ranging from 2 to 30 m, Younes et al. [29] explored the effects of different spatial resolutions on some controlling factors extracted from remote sensing data when drawing a gully susceptibility map based on pixel-based image analysis and concluded that 10 m was the optimal spatial resolution. Shruthi et al. [30] showed that the rule set developed for gully extraction based on object-oriented image analysis is sensitive to the spatial resolution of remote sensing images. However, these studies did not consider the impact of resolution on gully extraction under different classification systems based on different application purposes.
As one of the three major black soil regions in the world, the black soil region of Northeast China has been a major grain production area of China, producing rice, corn and soybeans. With topographical characteristics of rolling hills and some poor management practices, it has been suffering from serious soil erosion, threatening food security and agricultural sustainability [31], of which gully erosion is one of the main parts. According to the first national water conservancy census, gully erosion is widely distributed and mainly appears in sloping farmland [32]. Furthermore, gully erosion has been intensifying, as the number of gullies is increasing and the channels have been expanding in recent years [3].
The present study aimed to analyze the suitability of remote sensing data at different resolutions for different application purposes on the small watershed scale of the black soil region of Northeast China. This objective is addressed by (1) classifying the types of gullies according to application purposes based on field investigation, (2) establishing methods of visual interpretation and extracting gullies, (3) analyzing the suitability of different resolution data and selecting the optimal resolution based on performance analyses and (4) mapping gullies with the optimal pixel size.

Study Area
The study area is located in the southeast of Baiquan County, which is in the centralwestern part of Heilongjiang Province, belonging to the transition zone between Xiaoxing'an Mountains and Songnen Plain, at a longitude of 126 • 16 15 E to 47 • 18 52 E and latitude of 47 • 23 37 N to 47 • 26 47 N (Figure 1), covering an area of about 12.58 km 2 . It is situated in the middle of the black soil region and is covered by typical black soil, meadow soil, albic soil and dark brown soil. The climate is a warm, temperate, continental monsoon climate, with hot and wet summers and cold and dry winters (that is, four seasons and high seasonal precipitation). Additionally, the long-term average annual rainfall is 490 mm. The rolling hill topography is obvious, with an altitude of 240-280 m above sea level and a slope of 4-6 degrees. Due to the geographical environment and reclamation activities, gully erosion is widespread in this region [32]. In addition, ecological protection efforts that were started at the end of the last century improved the gully erosion to a certain extent [33]. Therefore, gully erosion in this area has gone through the natural development stage, the human intervention stage and the regulation stage, which is representative for studying the occurrence and development of gully erosion and its impact on the ecological environment in the black soil region in Northeast China.
Remote Sens. 2020, 17, x FOR PEER REVIEW 3 of 17 meadow soil, albic soil and dark brown soil. The climate is a warm, temperate, continental monsoon climate, with hot and wet summers and cold and dry winters (that is, four seasons and high seasonal precipitation). Additionally, the long-term average annual rainfall is 490 mm. The rolling hill topography is obvious, with an altitude of 240-280 m above sea level and a slope of 4-6 degrees. Due to the geographical environment and reclamation activities, gully erosion is widespread in this region [32]. In addition, ecological protection efforts that were started at the end of the last century improved the gully erosion to a certain extent [33]. Therefore, gully erosion in this area has gone through the natural development stage, the human intervention stage and the regulation stage, which is representative for studying the occurrence and development of gully erosion and its impact on the ecological environment in the black soil region in Northeast China.

Data Sources
The outline and internal details of gullies are the main factors determining their distribution and category. Whether they are covered or not and their clarity directly affect the interpretation accuracy of gullies. Time phase, range of spatial resolution and source of images were taken into account when selecting remote sensing images. Studies on gullies conducted in an environment similar to the study area have shown that images from late spring, autumn and early summer were preferable to those from other phases, as these images experienced minimal effects from vegetation cover or snowfall [34]. Additionally, images with a spatial resolution from sub-meter to 30 m (Landsat data) have been chosen for gully erosion research [18,26,28]. However, many kinds of images are not open-source and need to be purchased at a high price. Furthermore, it is quite difficult to find a data source with high-quality images from the desired time phase. Thus, Google Earth (GE) images at resolutions of 0.51, 1.02, 2.04, 4.08, 8.16, 16.32 and 32.64 m acquired on 9 May 2018 were used in this study (downloaded from the software 91weitu [35]).
The interpretation of remote sensing images needs to be supported and validated with ground data, which can be obtained from field investigations. Two field trips that lasted one week each were carried out in May 2018 and 2019, respectively. During the two field trips, as the main equipment, a portable GNSS with centimeter accuracy was used for precise positioning of the gullies. In the first field trip, some typical areas east of Baiquan County ( Figure 1c) were investigated, which provided a basis for the determination of the types of gullies and the establishment of interpretation marks. The second trip was to obtain the distribution of gullies in the validation areas (I, II and III). When gullies were measured in the field survey, we used a portable GNSS with centimeter accuracy and high precision to collect measurement points along the gully. The density of points increased when the terrain features changed sharply, and the density of points decreased at the place where the topographic features changed gradually without affecting the reflection of the topographic changes. In addition, there was at least one point at the head and tail of the gullies. Using the collected measurement points and GE images, gullies were drawn within the ArcGIS 10.5 software. Then, each gully was given a category attribute and area statistics were determined (Figure 1d).
In addition, auxiliary data (roads, residential sites, etc.) were also provided by the software 91weitu.

Visual Interpretation for Gully Mapping
According to the field investigation, there were different types of gullies based on different classification criteria (Table 1). Based on the development degree of the gullies, they could be divided into rill, ephemeral gully, permanent gully and modern incised valley [36]. Rills were not the object of this study since they can be eliminated by conventional tillage. Ephemeral gullies of about 30-50 cm in width cannot be completely eliminated by tillage, and in the next erosion process, ephemeral gullies would develop again in the same place. Through upstream erosion at the gully head, undercutting erosion at the gully bottom and lateral erosion at the gully slope, permanent gullies with a width and depth of more than 50 cm were formed. The profile of a modern incised valley was a concave arc different from that of a permanent gully. According to the activity of the gullies, they could be divided into active gullies and stable gullies. Gullies with lots of sediments, a large number of plants growing at the wall and bottom and shrubs, small arbors or artificial forests appearing were stable. According to whether gullies were treated or not, they could be divided into untreated gullies and treated gullies. Governors sorted out the boundaries of the gullies; flattened the wall and bottom; and planted willows, sylvestris or other shrubs at the edge, wall and bottom to treat gullies. Then, visual interpretation marks were established based on image features (shown in Table 1 with GE image at 0.5 m spatial resolution as an example).

Classification Criteria
Type of Gully Ground Photo GE Image (0.5 m) Image Feature

Severity and morphology
Ephemeral gully when the length was greater than or equal to 100 m. In view of the weak penetration of optical images, gullies in the forest were not included. In addition, different sections of the same gully might belong to different types. In this case, high-priority types were adopted. The priority order was as follows: modern incised valley > permanent gully > ephemeral gully; active gully > stable gully; untreated gully > treated gully. Clear and regular boundary, uniform texture inside when the length was greater than or equal to 100 m. In view of the weak penetration of optical images, gullies in the forest were not included. In addition, different sections of the same gully might belong to different types. In this case, high-priority types were adopted. The priority order was as follows: modern incised valley > permanent gully > ephemeral gully; active gully > stable gully; untreated gully > treated gully. the center line, and the length was calculated. The branch of a gully would be counted when the length was greater than or equal to 100 m. In view of the weak penetration of optical images, gullies in the forest were not included. In addition, different sections of the same gully might belong to different types. In this case, high-priority types were adopted. The priority order was as follows: modern incised valley > permanent gully > ephemeral gully; active gully > stable gully; untreated gully > treated gully. Clear and regular boundary, uniform texture inside the center line, and the length was calculated. The branch of a gully would be counted when the length was greater than or equal to 100 m. In view of the weak penetration of optical images, gullies in the forest were not included. In addition, different sections of the same gully might belong to different types. In this case, high-priority types were adopted. The priority order was as follows: modern incised valley > permanent gully > ephemeral gully; active gully > stable gully; untreated gully > treated gully. the center line, and the length was calculated. The branch of a gully would be counted when the length was greater than or equal to 100 m. In view of the weak penetration of optical images, gullies in the forest were not included. In addition, different sections of the same gully might belong to different types. In this case, high-priority types were adopted. The priority order was as follows: modern incised valley > permanent gully > ephemeral gully; active gully > stable gully; untreated gully > treated gully. Clear and regular boundary, uniform texture inside the center line, and the length was calculated. The branch of a gully would be counted when the length was greater than or equal to 100 m. In view of the weak penetration of optical images, gullies in the forest were not included. In addition, different sections of the same gully might belong to different types. In this case, high-priority types were adopted. The priority order was as follows: modern incised valley > permanent gully > ephemeral gully; active gully > stable gully; untreated gully > treated gully. the center line, and the length was calculated. The branch of a gully would be counted when the length was greater than or equal to 100 m. In view of the weak penetration of optical images, gullies in the forest were not included. In addition, different sections of the same gully might belong to different types. In this case, high-priority types were adopted. The priority order was as follows: modern incised valley > permanent gully > ephemeral gully; active gully > stable gully; untreated gully > treated gully. Clear and regular boundary, uniform texture inside the center line, and the length was calculated. The branch of a gully would be counted when the length was greater than or equal to 100 m. In view of the weak penetration of optical images, gullies in the forest were not included. In addition, different sections of the same gully might belong to different types. In this case, high-priority types were adopted. The priority order was as follows: modern incised valley > permanent gully > ephemeral gully; active gully > stable gully; untreated gully > treated gully. the center line, and the length was calculated. The branch of a gully would be counted when the length was greater than or equal to 100 m. In view of the weak penetration of optical images, gullies in the forest were not included. In addition, different sections of the same gully might belong to different types. In this case, high-priority types were adopted. The priority order was as follows: modern incised valley > permanent gully > ephemeral gully; active gully > stable gully; untreated gully > treated gully. Clear and regular boundary, uniform texture inside the center line, and the length was calculated. The branch of a gully would be counted when the length was greater than or equal to 100 m. In view of the weak penetration of optical images, gullies in the forest were not included. In addition, different sections of the same gully might belong to different types. In this case, high-priority types were adopted. The priority order was as follows: modern incised valley > permanent gully > ephemeral gully; active gully > stable gully; untreated gully > treated gully. ern incised valley, active gully and untreated gully. Gullies were linearized by extracting the center line, and the length was calculated. The branch of a gully would be counted when the length was greater than or equal to 100 m. In view of the weak penetration of optical images, gullies in the forest were not included. In addition, different sections of the same gully might belong to different types. In this case, high-priority types were adopted. The priority order was as follows: modern incised valley > permanent gully > ephemeral gully; active gully > stable gully; untreated gully > treated gully. Clear and regular boundary, uniform texture inside ern incised valley, active gully and untreated gully. Gullies were linearized by extracting the center line, and the length was calculated. The branch of a gully would be counted when the length was greater than or equal to 100 m. In view of the weak penetration of optical images, gullies in the forest were not included. In addition, different sections of the same gully might belong to different types. In this case, high-priority types were adopted. The priority order was as follows: modern incised valley > permanent gully > ephemeral gully; active gully > stable gully; untreated gully > treated gully. ern incised valley, active gully and untreated gully. Gullies were linearized by extracting the center line, and the length was calculated. The branch of a gully would be counted when the length was greater than or equal to 100 m. In view of the weak penetration of optical images, gullies in the forest were not included. In addition, different sections of the same gully might belong to different types. In this case, high-priority types were adopted. The priority order was as follows: modern incised valley > permanent gully > ephemeral gully; active gully > stable gully; untreated gully > treated gully. Clear and regular boundary, uniform texture inside ern incised valley, active gully and untreated gully. Gullies were linearized by extracting the center line, and the length was calculated. The branch of a gully would be counted when the length was greater than or equal to 100 m. In view of the weak penetration of optical images, gullies in the forest were not included. In addition, different sections of the same gully might belong to different types. In this case, high-priority types were adopted. The priority order was as follows: modern incised valley > permanent gully > ephemeral gully; active gully > stable gully; untreated gully > treated gully. Based on the interpretation marks and auxiliary data (road, residential sites, etc.), the processes of interpretation and vectorization of the gullies were completed in ArcGIS 10.5 software. Except for ephemeral gullies, the interpreted gullies included two attributes: area and category. Ephemeral gullies were only outlined for length, not area. When it was impossible to judge the type of gully, the type attributes of gullies were assigned as modern incised valley, active gully and untreated gully. Gullies were linearized by extracting the center line, and the length was calculated. The branch of a gully would be counted when the length was greater than or equal to 100 m. In view of the weak penetration of optical images, gullies in the forest were not included.
In addition, different sections of the same gully might belong to different types. In this case, high-priority types were adopted. The priority order was as follows: modern incised valley > permanent gully > ephemeral gully; active gully > stable gully; untreated gully > treated gully.

Selection of Optimal Spatial Resolution
By comparing the area of gullies from visual interpretation with that from the field survey data in three validation areas (I, II and III) and considering the difficulty of data Remote Sens. 2021, 13, 2367 6 of 16 acquisition and application requirements, the optimal resolution of visual interpretation was obtained.
Three indexes, namely precision, recall and F-score, were introduced to help evaluate the performance of visual interpretation. Precision indicates how many areas were correctly interpreted compared with the interpreted gully area. Recall indicates how many areas were correctly interpreted compared with the actual gully area. F-score is a comprehensive measure of precision and recall. Ideally, all three indexes should be 0. The calculation formulas are as follows: where A is the area of interpreted gullies and B is the corresponding actual gully area.

Overall Situation of Gully Interpretation
The distribution maps of gullies are shown in Figure 2. Here, the higher the spatial resolution of images is, the easier is it to obtain pure pixels. Thus, the contour, shape and internal details of gullies were clearer and more independent in the process of interpretation with images of a higher spatial resolution. On the contrary, with images of a lower spatial resolution, gullies were vague shapes that could only be identified simply or could not be identified at all. In the image with a spatial resolution of 8.16 m, the shape and texture of the gullies became very blurred, but they still could be roughly identified and outlined by the overall trend and location. The feature information of gullies almost disappeared when the resolution of the image was reduced to 16.32 m, making it impossible to interpret gullies.
The distinction between field paths and gullies needed to be paid attention to [8]. The field investigation showed that the field paths were flat, while gullies had a certain undercutting depth. In addition, the boundary of the field paths was smooth, while that of gullies was tortuous. In the process of interpretation, gullies and field paths showed significantly different characteristics in images with a resolution of 0.51 and 1.02 m, enabling them to be distinguished. However, gullies might show almost the same characteristics as field paths in images with a lower spatial resolution, preventing them from being correctly recognized.
The accuracy of gully interpretation was evaluated in three verification areas, shown in Figure 3 and Table 2. The interpretation in images with a resolution of 0.52, 1.02 and 2.04 m was in good agreement with the field investigation. The average F-score of the three verification areas was higher than 95%, and the F-score of each verification area was higher than 90%. Among them, images with the resolution of 1.02 m had the best performance. When the resolution reached 4.08 m, the performance became evidently worse.

Interpretation of Different Types of Gullies
Images of different resolutions could interpret different types of gullies. The number of different types of gullies in the three verification areas is shown in Figure 4. According to the interpretation process and verification of the number of gullies, all types of gullies could be well identified in images with a spatial resolution of 0.51 m. It began to be difficult to identify ephemeral gullies in images with a spatial resolution of 1.02 m. It was completely impossible to distinguish ephemeral gullies and was difficult to determine the types of active gully and stable gully in images with a spatial resolution of 2.04 m. When the spatial resolution reached 4.08 m, the classification of permanent gullies and modern incised valleys became difficult, as did the determination of untreated gullies and treated gullies. With images of 8.16 m spatial resolution, all types of gullies could not be determined, and The accuracy of gully interpretation was evaluated in three verification areas, shown in Figure 3 and Table 2. The interpretation in images with a resolution of 0.52, 1.02 and 2.04 m was in good agreement with the field investigation. The average F-score of the

Interpretation of Different Types of Gullies
Images of different resolutions could interpret different types of gullies. The number of different types of gullies in the three verification areas is shown in Figure 4. According to the interpretation process and verification of the number of gullies, all types of gullies could be well identified in images with a spatial resolution of 0.51 m. It began to be difficult to identify ephemeral gullies in images with a spatial resolution of 1.02 m. It was completely impossible to distinguish ephemeral gullies and was difficult to determine the types of active gully and stable gully in images with a spatial resolution of 2.04 m. When the spatial resolution reached 4.08 m, the classification of permanent gullies and modern incised valleys became difficult, as did the determination of untreated gullies and treated gullies. With images of 8.16 m spatial resolution, all types of gullies could not be determined, and gully identification was extremely difficult. Therefore, the interpretation with images of 8.16 m spatial resolution is not shown below.  In order to further study the ability to interpret different types of gullies from images with different spatial resolutions, the extracted areas of different types of gullies in the three verification areas are shown in Figure 5, and the accuracy evaluation is shown in Tables 3-5. The images with a spatial resolution of 0.51 and 1.02 m had excellent performance for the interpretation of different types of gullies under any classification system, and the average F-scores of the three validation areas were all above 98%. When the spatial resolution of the images began to decline from 2.04 m, the interpretation ability for different types of gullies under all classification systems showed a downward trend. Except for the influence of type misclassification, when the spatial resolution of the image was 4.08 m, the ability to interpret modern incised valleys was significantly lower than that for permanent gullies, and the ability to interpret untreated gullies was significantly higher than that for treated gullies. When the spatial resolution of the image was 2.04 m, the ability to interpret active gullies was significantly higher than that for stable gullies. Remote Sens. 2020, 17, x FOR PEER REVIEW 9 of 17 In order to further study the ability to interpret different types of gullies from images with different spatial resolutions, the extracted areas of different types of gullies in the three verification areas are shown in Figure 5, and the accuracy evaluation is shown in Tables 3-5. The images with a spatial resolution of 0.51 and 1.02 m had excellent performance for the interpretation of different types of gullies under any classification system, and the average F-scores of the three validation areas were all above 98%. When the spatial resolution of the images began to decline from 2.04 m, the interpretation ability for different types of gullies under all classification systems showed a downward trend. Except for the influence of type misclassification, when the spatial resolution of the image was 4.08 m, the ability to interpret modern incised valleys was significantly lower than that for permanent gullies, and the ability to interpret untreated gullies was significantly higher than that for treated gullies. When the spatial resolution of the image was 2.04 m, the ability to interpret active gullies was significantly higher than that for stable gullies.
Considering the different geographical environments in different regions, the interpretation performance of an image with one spatial resolution might be different in different verification areas. The stability of image interpretation in different regions is very important; therefore, relevant exploration was carried out. When the spatial resolution was 4.08 m, the performance of interpreting modern incised valleys in verification area II was quite poor compared with that in the other two verification areas. Therefore, it was considered that the image with a spatial resolution of 4.08 m was not stable in the interpretation of gullies classified by their development degree. Similarly, images with a spatial resolution of 2.04 m or lower were not stable enough in the interpretation of gullies classified by activity and treatment.      Considering the different geographical environments in different regions, the interpretation performance of an image with one spatial resolution might be different in different verification areas. The stability of image interpretation in different regions is very important; therefore, relevant exploration was carried out. When the spatial resolution was 4.08 m, the performance of interpreting modern incised valleys in verification area II was quite poor compared with that in the other two verification areas. Therefore, it was considered that the image with a spatial resolution of 4.08 m was not stable in the interpretation of gullies classified by their development degree. Similarly, images with a spatial resolution of 2.04 m or lower were not stable enough in the interpretation of gullies classified by activity and treatment.

Discussion
A high spatial resolution was helpful to interpret gullies, but an excessively high resolution increased noise. Although the shadow caused by gully erosion, which could help determine the boundaries of gullies, was clear in the high-resolution images [8], the shadow of vegetation growing in gullies was too clear, which covered up the outline of the gullies and made it difficult to draw them [28]. Therefore, the interpretation accuracy of the image with a spatial resolution of 1.02 m was the highest when not considering the type of gully.
When choosing the spatial resolution of the image, one should consider not only the gully interpretation accuracy but also the price of the image and the application purpose of the research. Generally speaking, the higher the spatial resolution of an optical image is, the higher the price is. Taking the company Beijing Lanyu Fangyuan Information Technology Co., Ltd. [37], which specializes in providing satellite images and other image processing services in China, as an example, the prices of some satellite images widely used and currently in service are shown in Table 6. It should be noted that the interpretation accuracy of gullies became significantly worse in general when the spatial resolution was 4.08 m or lower, so only the prices of images with a spatial resolution of 0.5-2.1 m are listed. Based on the performance of gully interpretation and the prices of the images, images with a resolution of 2.04 m are the best choice when not considering the types of gullies interpreted, which could meet the general requirements of relevant professionals. In fact, many studies have used SPOT images with a spatial resolution of 2.5 m to interpret gullies [28,38]. When finer monitoring in a large area was carried out, images with a resolution of 1.02 m could be selected. When studying soil loss caused by gullies at different development stages, the types of gullies that need to be interpreted are ephemeral gully, permanent gully and modern incised valley. To extract the information of these three types of gullies completely, images with a spatial resolution of 0.51 m are the only choice. However, research on ephemeral gullies has little practical significance, as they have narrow width and shallow depth, do not hinder general tillage-or sometimes disappear due to tillage-and often appear repeatedly at the same position. In this situation, images with a spatial resolution of 1.02 and 2.04 m could be considered, and images with a spatial resolution of 2.04 m are better in that they are cheaper and still meet the requirements. Pullman showed that ZY-3 satellite images with a spatial resolution of 2 m were an ideal choice for interpreting permanent gullies and modern incised gullies [28]. When studying the activity of gullies to formulate relevant governance policies, the types of gullies that need to be interpreted are active gullies and stable gullies. Images with a spatial resolution of 0.51 and 1.02 m can both be used, but considering the practical significance, price factors and interpretation performance, images with a spatial resolution of 1.02 m are undoubtedly better. When it is necessary to monitor the treatment of gullies, gullies first need to be identified as untreated gullies or treated gullies. Similar to the selection of the spatial resolution of images used to interpret active gullies and stable gullies, images with a spatial resolution of 0.51 and 1.02 m could be listed as candidates, but the images with a spatial resolution of 1.02 m are more suitable. It should be noted that the abundance of auxiliary data might reduce the requirement for spatial resolution. Auxiliary data include high-precision DEMs, fine road network data and vegetation fractions. Studies showed that SPOT images with a spatial resolution of 2.5 m and NDVI data could be used to interpret active gullies, semiactive gullies and stable gullies, and the result was ideal [39]. Moreover, the spectral and temporal resolutions of images are important factors affecting the characteristics of surface features in remote sensing images, besides spatial resolution [40,41]. The improvement of spectral resolution could provide a more precise surface spectrum, which could be used as auxiliary data to better interpret gullies. With the improvement of temporal resolution, more dynamic information could be obtained. If images before and after a rainstorm could be obtained, the rainwater accumulation shown in the image after the rainstorm could be used as auxiliary information to interpret gullies. The interpretation of gullies was only studied from the perspective of spatial resolution with limited auxiliary data in this paper. Future research needs to take spectral resolution, temporal resolution and enrichment of auxiliary data into account.
GE images, a type of optical image, were the only satellite image data used for mapping gullies in this study. Some other studies have used GE images to map gullies and confirmed the availability of GE images [42][43][44]. Compared with other high-resolution optical images, GE images can be used to evaluate gullies in a very large area without any cost of image acquisition. However, the dates of GE images are often not homogeneous across a region provided for some platforms, which remains a limitation of the application of GE images [45]. In addition, gully mapping was restricted to agricultural land, as forest and other surface cover precluded interpretation of gullies beneath [46].
Comparing the GE images of the two time phases used in the study, it was found that there were almost no changes in the types and distribution of ground objects, which shows that the states of the gullies in the study area were basically the same at the two time points. Therefore, the lag of verification data acquisition had little effect on the results in this study. However, studies have shown that gullies can be greatly developed in a short time if conditions are favorable due to rainstorms and freeze-thaw erosion [47,48]. Therefore, the acquisition time of validation data should be as close as possible to the acquisition time of images used to interpret gullies, and it is better not to use an interval including a rainy season or a freeze-thaw period when relevant research of gullies is carried out.
This study was carried out in a typical black soil region of Northeast China. The method and results could be extended to other black soil areas in Northeast China for the following reasons: (1) there were similarly distributed rules for gully erosion across the entire black soil region, and (2) the study area had all types of gullies in the black soil area due to a long history of gully development and management. However, some adjustments need to be made to this method to obtain results suitable for local conditions when attempting to extend the present approach to areas other than the black soil region. For example, it is necessary to reselect the time phase of images and reclassify the types of gullies according to field investigation.
In response to the problem of gullies swallowing cultivated land, land managers are looking for effective land management measures to control gully erosion, such as the terracing of hillslopes [49]. Land use planning and land degradation are interrelated and interact with each other [50]. Further research is needed to clarify the relationship between gully development and land use planning, which can provide suggestions for gully control. Indicators are often selected and calculated with remote sensing and GIS techniques to reflect the relationship between humans and land [51]. Based on this study, the spatial and temporal distribution of different types of gullies could be obtained by interpreting remote sensing images with optimal spatial resolution. Then, the indicators of gully severity such as gully intensity and density can be used to quantitatively study the relationship between land use planning and gully development.

Conclusions
The main purpose of this study was to reveal the applicability of images of different spatial resolutions for the interpretation of different types of gullies so as to select the optimal spatial resolution under different application requirements. The study draws the following conclusions: Firstly, improvement of the spatial resolution will make the contour and internal details of gullies clearer. However, at the same time, it will increase the noise information. Based on this, the interpretation performance of images with a spatial resolution of 1.02 m is the most satisfactory, without considering the type of gully. Additionally, the interpretation result of the image with a spatial resolution of 2.04 m remained excellent, with all F-scores above 95%. Therefore, images with a spatial resolution of 2.04 m are most suitable for general research, considering acquisition cost. When the spatial resolution falls to 16.32 m, gullies cannot be interpreted completely. Secondly, images with different spatial resolutions have different interpretation abilities for different types of gullies. The interpretation of all types of gullies from images with a spatial resolution of 0.51 m is outstanding. Different from images with a 0.51 m spatial resolution, it is difficult to identify ephemeral gullies from images with a spatial resolution of 1.02 m. However, ephemeral gully research has little practical significance due to minor hazards and temporary occurrences. Therefore, images with a spatial resolution of 1.02 m are the most universally appropriate choice. Thirdly, when the spatial resolution is 2.04 m or lower, the ability to interpret some other types of gullies in addition to ephemeral gullies is weakened. Therefore, the spatial resolution should be selected according to the types of gullies needed in practical application. When the spatial resolution of an image is 8.16 m, it is impossible to identify any types of gullies, and it is also very difficult to identify gullies at all.
There are some problems in this study that require further research. Besides the spatial resolution, the spectral resolution and temporal resolution also have an impact on the interpretation of gullies. In addition, the selection and enrichment of auxiliary data is another factor that changes the interpretation performance. Furthermore, different sections of the same gully may belong to different types. In this study, a simple classification was carried out, but a finer classification method is needed.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy.