Soil erosion is the combined result of natural factors and human factors, and has become one of the most significant environmental problems worldwide, threatening agricultural productivity and food security (DeGraffenried & Shepherd 2009; Zhang et al., 2010; Verachtert et al., 2011; Zhao et al., 2018) . Efficient intervention to control soil erosion requires assessment models for predicting the spatial location and intensity of degradation (Cohen et al., 2005; Zhang et al., 2012) . Among the factors affecting soil erosion (rainfall, soil erodibility, vegetation cover, topographic factors, land cover, etc.), rainfall and vegetation cover are most active factors. Rainfall is the main driving factor (Zhang & Fu, 2003) . Raindrop splashing and separating soil particles, and runoff eroding and transporting lead to soil erosion. However, not all rainfall events can cause soil erosion, and only those being able to generate sufficient runoff to transport sediment are erosive rainfall (Xie et al., 2002; Liu et al., 2007) . Hancock (2009) demonstrated that annual sediment output was variable and non-linear with increased rainfall amount and intensity. Vegetation is a positive factor in preventing soil erosion, and is one of the most active factors (Zhan & Wang, 2001) . Inappropriate land use and destruction of vegetation will lead to increasing soil erosion (Huang et al., 2004) . Vegetation can improve the ability of soil resistance to erosion, such as lowering rainfall kinetic energy by vegetation stem, slowing runoff velocity by plant stems and litter, increasing ability of soil resistance to erosion by plant roots (Zhang et al., 2000) . Vegetation weaken raindrops hitting the ground, increase ground roughness and disperse the raindrops power in the canopy; carious vegetation also increase soil organic matter content and further improve the physical and chemical properties of soil (Gilly & Risse, 2000; Zhang & Liang, 1996) . However, vegetation obviously changes over time, and vegetation cover will also be significantly different in different seasons, thus the ability preventing soil erosion varies during the year.
Vegetation cover is usually extracted from satellite images in current erosion risk assessment methods (Mutekanga et al., 2010) . It holds an important role in protecting soil (Vrieling et al., 2008; Menéndez-Duarte et al., 2009) . As vegetation grows, vegetation cover which is an important impact factor in soil erosion assessment significantly varies. Consequently, the temporal phase of images impacts accuracy of extracted soil erosion risk information (Gao & Wang, 2004) . However, there is not a clear guideline to choose the most appropriate temporal image data for soil erosion risk remote sensing assessment. Usually, the images in growing season, especially periods of high intensity rainfall (Mutekanga et al., 2010; Zhang et al., 2010; Zhang et al., 2012) , are selected and used in the studies. In addition, after summarizing the experience of previous soil erosion remote sensing investigation, some scholars put forward some methods used to identify best temporal phase of images. For example, it is the best temporal phase when the tree and shrub species have begun to vigorously grow and the crops have not yet covered the surface in the loess plateau (Cai et al., 2002; Qiao, 2003) . Because the significant chlorophyll effect of vegetation results in prominent vegetation information, while the loess information is also very distinct, these two kinds of objects can form a strong contrast. In this time, the effect is best for visual interpretation. This method will allow researchers to distinguish more vegetation types, resulting in more effective data interpretation. However, in soil erosion remote sensing study, the erosion risk or intensity is not intuitive information and cannot be directly interpreted from the digital image. That is, the best temporal phase for visual interpretation is not necessarily suitable for soil erosion assessment. Therefore, the vegetation of this period may not be necessarily representative for erosion risk assessment.
For a certain study area, in addition to rainfall and vegetation, the other impact factors can be considered as constants within a certain period of time. Therefore, accurate analysis of temporal matching relationship of rainfall and vegetation during the year is of great significance for assessing soil erosion risk and optimizing soil and water conservation. Langbein & Schumm (1958) reveal a complex relationship that the rainfall directly result in soil erosion and simultaneously facilitate vegetation growth, as an indirect way to reduce the effects of surface erosion. Vrieling et al. (2008) analyzed the severe erosion periods of Cerrados regions in Brazil using moderate resolution imaging spectroradiometer (MODIS) normalized difference vegetation index (NDVI) time series data and high-intensity rainfall showed by Tropical Rainfall Measuring Mission (TRMM) three-hourly data. The purpose of this paper is to find a good way for choosing remote sensing images which is used to assess soil erosion risk, that is, to identify optimal vegetation cover for annual soil erosion risk assessment. TRMM-3B43 monthly rainfall data and MODIS-NDVI data are used to study the temporal matching relationship of rainfall and vegetation in the Huangfuchuan basin and an indicator RV is developed to aid remote sensing image selection in practical applications.
2. Materials and Methods
2.1. Study Area
The Huangfuchuan Basin is located in the middle reaches of the Yellow River, and the main source of coarser sediment in the Yellow river basin (Figure 1). The average annual sediment transport volume is 0.505 billion tons. The soil erosion area is 3215 km2 (39˚11' - 40˚3'N, 110˚15' - 111˚19'E), accounting for 99% of the total area.
Figure 1. Location of study area.
The temporal variation of rainfall is determined using time-series TRMM 3B43 products. These data are available from 1998 and were obtained through the TRMM Online (http://daac.gsfc.nasa.gov/data/datapool/TRMM). TRMM 3B43 is a monthly rainfall product and available at a 0.25˚ grid. Its spatial resolution is low, but it is considered that a good indication for temporal rainfall distribution (Vrieling et al., 2008; Zhang et al., 2012; Zhao et al., 2018) . This dataset is selected because rain gauge measurements do not always exit or are often hard to obtain, while TRMM data is available for any region between latitudes 50˚N and 50˚S, furthermore, its spatial distribution is more authentic relative to the interpolated results of point measured data. A point-measurement value is hard to represent a broader range of rainfall distribution, although it is more accurate at this point location (Wilheit et al., 1977) . In the observation of convective precipitation, this point data is more unrepresentative. Radar is also an alternative, however, it can only observe limited area, and its disadvantage is uncertainty, low coverage, high cost (Cheng et al., 1993) . The meteorological satellites can provide macro and large-scale global data, and it has become an important observation means of global precipitation (WMO, 1993) . TRMM project is the first active remote sensing of earth atmosphere by satellite and can collect more accurate information on the spatial structure of rain clouds.
NDVI has been widely used in the study of vegetation remote sensing and phenology. It is a best indicator of plant growth status and spatial distribution density, and can be easily extracted from remote sensing images with a variety of spatial resolution. Although NDVI value is disturbed by soil background reflectance (Escadafal, 1994) and vegetation growth (De Jong, 1994) , it can still reflect the spatial and temporal variation of vegetation (Carlson & Ripley, 1997; Zhang et al., 2019) . It often is used to assess protective vegetation cover (De Jong et al., 1999; Jain & Goel, 2002; Symeonakis & Drake, 2004) . In this paper, MODIS imagery from the Terra platform is used to determine the temporal variation of NDVI during the year. The MOD13Q1 product, representing a 16-day vegetation index composite, is used and it is calculated on a pixel-by-pixel basis (https://wist.echo.nasa.gov/api/), in which images are selected or excluded based on quality, cloud cover and viewing geometry (Huete et al., 2002) . This product always provides a well-timed NDVI value for each pixel even during poor conditions, such as persistent cloud cover. And it is widely used to monitor ecological environment.
The sediment transport quantity is the sediment mass passing through a section of the river in a certain period of time, and the unit is ton. This paper uses it to indicate the erosion and sediment yield of the controlled area. The data include monthly and annual sediment transport data of hydrological station monitoring the
3.1. The Definition of Optimal Vegetation Cover
In the arid and semi-arid areas where low and concentrated rainfall, the time high intensity rainfall occurs is not fixed in a certain period. While the vegetation cover varies, resulting in the protective ability of vegetation changes for soil. If the high intensity rainfall occurs when vegetation protective ability is maximal, the most server erosion may not appear, thereby, the image of this period is not necessary representative for assessing soil erosion.
Similarly, in northern
3.2. Analysis of Time-Series Rainfall and NDVI
Rainfall and vegetation are the two factors with the greatest changes for annul soil erosion risk assessment. In order to study the matching relationship of variation of rainfall and vegetation in the year, several parameters are chosen to describe the time series data, that is, kurtosis, skewness and peak position. Compared to the normal distribution, kurtosis reflects the sharpness or flatness of a distribution. A positive value indicates a more sharp distribution. A negative value indicates a more flat distribution. Calculation formula (1) as follows:
where, KU is kurtosis; S is the standard deviation of time series data.
Skewness reflects the asymmetry degree relative to average value in a distribution. A positive value indicates the asymmetric part has a positive trend. A negative value indicates the asymmetric part has a negative trend. Calculated formula (2) as follows:
where SK is skewness; S is the standard deviation.
In this study area, rainfall and vegetation time-series are single peak (Zhang & Qin, 2016) , the maximum value are considered as the peak position.
3.3. Design of Optimal Vegetation Cover Indicator
In practical applications, it is impossible and unnecessary to firstly purchase all remote sensing images of rainy season for soil erosion risk assessment. The correct approach is to firstly identify the optimal temporal phase, and then purchase the image data. Therefore, it is necessary to develop an indicator of the optimal temporal phase for remote sensing images selection. For a specific study area, the erosion impact factors, such as soil type, topography, tillage practices, etc. can be considered as constant for annual soil erosion assessment. Rainfall and vegetation significantly change in the year, and selected to build the indicator, assisting remote sensing image selection. NDVI (Normalized Difference Vegetation Index) is an important parameter (see Equation (3)), often used to detect vegetation growth status and vegetation cover. Here, it is used to be on behalf of vegetation condition. Rainfall data is the TRMM 3B43 data set.
where NIR is the near infrared band of the image, and RED is red band.
Some studies have proved that rainfall which less than
where, R is the original rainfall, is the processed rainfall, i is month number.
Since the units and magnitude of rainfall and NDVI data are inconsistent, in order to express better the matching relationship, they are converted to relative values by standardizing with the sum of the year. Then an index RV is developed to help indicate the relative soil erosion risk.
where, is the standardized rainfall. The values reflect the percentage of erosive rainfall per month in the year; is the standardized
NDVI, the values indicate the relative size of vegetation cover for each moth in the year. Therefore, RV can indicate the relative erosion risk during the year. For example, the greater RV values indicate that relatively large rainfall encounters relatively small vegetation cover during the year, and the erosion risk will be greater. When the RV is negative, it still reflects the relative magnitude of erosion risk, However it indicates that erosion is almost impossible.
It should be noted that, the data is standardized with the sum of the year, and it is not comparable year to year. In addition, our aim is just to indicate the probability of erosion occurrence in each period, only the quantity and not intensity of rainfall is considered.
4. Results and Discussions
4.1. Comparison of Parameters
The kurtosis of 2004-2013 rainfall and NDVI time series data are calculated and shown in Figure 2. It can be seen, the kurtosis of NDVI is relatively stable, and ranges between −2 and 1. Rainfall kurtosis values are all greater than NDVI in each year, and the fluctuations are more obvious. That illustrates that rainfall is more concentrated and the concentration degree changes each year.
The skewness is shown in Figure 3. The changing trend of NDVI is similar with rainfall, but the amplitude is less than rainfall. That illustrates that rainfall may promote vegetation growth, but the change of NDVI is relatively slow.
Peak position is shown in Figure 4. Although both are located in rainy season, the changes of rainfall are still greater than NDVI.
The analysis of three parameters indicates that the temporal matching relationships of rainfall and NDVI vary in each year, and changing degree of rainfall is more evident than NDVI. So, the best temporal phase of image should be different in each year for soil erosion remote sensing assessment. It is necessary to analyze the specific circumstances of each year in order to find the most appropriate image.
Figure 2. Kurtosis of precipitation and NDVI time series.
Figure 3. Skewness of precipitation and NDVI time series.
Figure 4. Peak positions of precipitation and NDVI time series.
4.2. The Matching Relationship of Rainfall and NDVI
The curves of 2004-2013 rainfall and NDVI are shown in Figure 5. Due to the impact of climate, rainfall and other natural conditions, variation of NDVI is slightly different in each year. Basically, vegetation starts to rapidly grow from April, and reaches the peak between July and September, then began to decline. However, the rainfall is very different. Although it is mainly concentrated between April and September, there is significant randomness not only from the concentration degree but peak position. For example, the rainfall is most concentrative in 2013, and the rainfall between June-September accounts for 92.47% of total rainfall of the year; while that of the rainy season between June and August making up only 63.34% in 2011.
The RV values are calculated and shown in Figure 6. The protection ability of vegetation has certain differences even in three months of rainy season. Therefore, maximum of rainfall and RV did not appear in the same month. In other words, erosion risk is not necessarily most severe even if the rainfall is maximal
Figure 5. Temporal matching relations of precipitation and vegetation time series during the year of 2004-2013.
Figure 6. RV distribution of the year after avoiding erosive rainfall.
in this month, because the protection ability of vegetation also increased. According to the definition of optimal vegetation cover, the best temporal phase will be different due to the matching relationships. Because RV indicates the relative probability of erosion occurrence and erosion is very concentrated in the study area, the RV maximum will be corresponding to the phase of the optimal vegetation cover. Therefore, the best phase are in the rainy season, but not fixed in a certain month. From the calculation results of many years, RV can intuitively assist to select the best phase for remote sensing assessment of soil erosion.
4.3. Comparison of RV and Average Annual Sediment Transport
The sediment transport data monitored in river sections can be used to illustrate the amount of soil erosion. In this study, the average annual sediment transport data derived from the main hydrological stations in the lower reaches of the study area is used to compare with RV. The result is shown in Figure 7(a). We can see that the synergistic change between them is obvious, especially during the rainy season. During the period when rainfall can cause erosion, that is, from May to September, the correlation of RV and average annual sediment transport of years (ASTY) is shown in Figure 7(b). It can be seen that their correlation coefficient reaches 0.89.
The related studies have shown that the annual rainfall causing erosion in the study area is mostly a small part of flood season (Zhang & Qin, 2016) . During this period, it is dominated by water erosion. The sediment transport of rainy season accounts for 80% of the annual sediment transport in this study area (Wang et al., 1983). Therefore, RV can well reflect the relatively size of erosion risk and sediment.
There is almost erosive rainfall during January-April and October-December. Therefore, the sediment transport is caused by other factors during this period. RV can not indicate its changes. Related studies have shown that wind erosion dominates the region from March to May, especially the relatively small particle size is greatly affected by wind, and its erosion modulus exceeds water erosion and gravity erosion (Tang et al., 2001) .
For soil erosion by water, the vegetation cover corresponding to the most severe erosion is very vital in this study area due to the concentrated precipitation and erosion. Based on the sediment transport data of years, it is enough to demonstrate that the importance of the period corresponding to the maximum RV which is in best temporal phase for water erosion study.
Among many factors affecting soil erosion, the rainfall is the main driving factor
Figure 7. (a) Comparison of RV and average annual sediment transport of years (ASTY); (b) A scatter plot to show the correlation of RV and ASTY.
causing soil erosion (Zhang & Fu, 2003) , and the vegetation cover is the most positive factor protecting soil (Zhan & Wang, 2001) . The other factors can be considered as invariables for annual soil erosion assessment. So the temporal matching relationship between these two factors can indicate the relative erosion risk of different periods during the year. We analyze the matching relationship aiming to identify the optimal period of image which will be used in annual soil erosion assessment based on remote sensing. Since the rainfall of this study area is not abundant, and heavy rainfall is concentrated, the vegetation cover corresponding to maximal relative erosion risk is of great significance for studying soil erosion by water. As the greater variation of annual rainfall distribution in this region, therefore the temporal matching relationship with vegetation cover varies each year, resulting in the best phase is not fixed in a certain month.
The proposed RV index based on the matching relationships between temporal rainfall and NDVI can indicate the relative erosion risk during the year. The maximal RV value has a good correlation with the optimal vegetation cover, and will effectively assist selecting the images for soil erosion remote sensing assessment. It is a significant contribution to soil erosion research and facilitates erosion research in the future, showing good potential for successful application in other places where heavy rainfall is concentrated.
We acknowledge the financial support provided by the National Science and Technology Platform Construction Project (2005DKA32300), the Major Research Projects of the Ministry of Education (16JJD770019), the Science and Technology Development Plan Project of Henan Province, China (121100111300) and the cooperation base open fund of the Key Laboratory of Geospatial Technology for the middle and lower Yellow River regions and CPGIS (JOF 201602). We especially thank the anonymous reviewers for their constructive comments and insightful suggestions that greatly improved the quality of this manuscript.
 Cai, J. Q., Ren, Z. Y., & Li, Y. C. (2002). Discussion on Related Technical Issues of Remote Sensing Investigation on Soil Erosion (in Chinese, with English Abstract). Bulletin of Soil and Water Conservation, 22, 45-47.
 Carlson, T. N., & Ripley, D. A. (1997). On the Relation between NDVI, Fractional Vegetation Cover, and Leaf Area Index. Remote Sensing of Environment, 62, 241-252.
 Cheng, M., Brown, R., & Collier, C. G. (1993). Delineation of Precipitation Areas Using Meteosat Infrared and Visible Data in the Region of the United Kingdom. Journal of Applied Meteorology, 32, 884-898.
 Cohen, M. J., Shepherd, K. D., & Walsh, M. G. (2005). Empirical Reformulation of the Universal Soil Loss Equation for Erosion Risk Assessment in a Tropical Watershed. Geoderma, 124, 235-252.
 De Graffenried, J. B., & Shepherd, K. D. (2009). Rapid Erosion Modeling in a Western Kenya Watershed Using Visible near Infrared Reflectance, Classification Tree Analysis and 137Cesium. Geoderma, 154, 93-100.
 De Jong, S. M. (1994). Derivation of Vegetative Variables from a Landsat TM Image for Modelling Soil Erosion. Earth Surface Processes and Landforms, 19, 165-178.
 De Jong, S. M., Paracchini, M. L., Bertolo, F., Folving, S., Megier, J., & De Roo, A. P. J. (1999). Regional Assessment of Soil Erosion Using the Distributed Model SEMMED and Remotely Sensed Data. Catena, 37, 291-308.
 Escadafal, R. (1994). Soil Spectral Properties and Their Relationships with Environmental Parameters: Examples from Arid Regions. In: J. Hill, & J. Mégier (Eds.), Imaging Spectrometry—A Tool for Environmental Observations (pp. 71-87), Dordrecht: Kluwer Academic Publishers,.
 Gao, Y. H., & Wang, H. Q. (2004). Application of Remote Sensing Technology in Dynamic Monitoring of Soil Erosion (in Chinese, with English abstract). Soil and Water Conservation in China, 1, 33-34.
 Hancock, G. R. (2009). A Catchment Scale Assessment of Increased Rainfall and Storm Intensity on Erosion and Sediment Transport for Northern Australia. Geoderma, 152, 350-360.
 Huang, Z. L., Chen, L. D., Fu, B. J., & Wu, X. L. (2004). Effect of Soil Erosion Reduction and Time Changing of Different Vegetation Types in Semi-Arid Loess Hilly and Gully Region (in Chinese, with English Abstract). China Water Resources, 20, 38-40.
 Huete, A., Didan, K., Miura, T., Rodriguez, E. P., Gao, X., & Ferreira, L. G. (2002). Overview of the Radiometric and Biophysical Performance of the MODIS Vegetation Indices. Remote Sensing of Environment, 83, 195-213.
 Jain, S. K., & Goel, M. K. (2002). Assessing the Vulnerability to Soil Erosion of the Ukai Dam Catchments Using Remote Sensing and GIS. Hydrological Sciences Journal, 47, 31-40.
 Langbein, L. B., & Schumm, S. A. (1958). Yield of Sediment in Relation to Mean Annual Precipitation. Transactions, American Geophysical Union, 39, 1076-1084.
 Mutekanga, F. P., Visser, S. M., & Stroosnijder, L. (2010). A Tool for Rapid Assessment of Erosion Risk to Support Decision-Making and Policy Development at the Ngenge Watershed in Uganda. Geoderma, 160, 165-174.
 Symeonakis, E., & Drake, N. (2004). Monitoring Desertification and Land Degradation over Sub-Saharan Africa. International Journal of Remote Sensing, 25, 573-592.
 Tang, Z. H., Cai, Q. G., Li, Z. W. et al. (2001). Study on Interaction among Wind Erosion, Hydraulic Erosion and Gravity Erosion in Sediment-Rock Region of Inner Mongolia. Journal of Soil and Water Conservation, 15, 25-29.
 Verachtert, E., Maetens, W., Van Den Eeckhaut, M., Poesen, J., & Deckers, J. (2011). Soil Loss Rates Due to Piping Erosion. Earth Surface Processes and Landforms, 36, 1715-1725.
 Vrieling, A., De Jong, S. M., Sterk, G., & Rodrigues, S. C. (2008). Timing of Erosion and Satellite Data: A Multi-Resolution Approach to Soil Erosion Risk Mapping. International Journal of Applied Earth Observation and Geoinformation, 10, 267-281.
 Wilheit, T. T., Chang, A. T. C., Rao, M. S. V. et al. (1977). A Satellite Technique for Quantitatively Mapping Rainfall Rates over the Oceans. Journal of Applied Meteorology, 16, 551-560.
 Zhan, X. G., & Wang, P. (2001). Research on Dynamical Supervision Design of Water and Soil Loss in Three Gorges Reservoir Area Based on RS and GIS (in Chinese, with English Abstract). Journal of Yangtze River Scientific Research institute, 18, 41-44.
 Zhang, G. H., & Liang, Y. M. (1996). A Summary of Impact of Vegetation Coverage on Soil and Water Conservation Benefit (in Chinese, with English Abstract). Research of Soil and Water Conservation, 3, 104-110.
 Zhang, X. C., Shao, M. A., Huang, Z. B., & Lu, Z. F. (2000). An Experimental Research on Soil Erosion and Nitrogen Loss under Different Vegetation Cover (in Chinese, with English Abstract). Acta ecologica Sinica, 20, 1038-1044.
 Zhang, X. W., Qiu, F., & Qin, F. (2019). Identification and Mapping of Winter Wheat by Integrating Temporal Change Information and Kullback-Leibler Divergence. International Journal of Applied Earth Observation and Geoinformation, 76, 26-39.
 Zhang, X. W., Wu, B. F. et al. (2012). Soil Erosion Risk and Its Spatial Pattern in Upstream Area of Guanting Reservoir. Environmental Earth Sciences, 65, 221-229.
 Zhao, H. G., Tang, Y. Y., & Yang, S. T. (2018). Dynamic Identification of Soil Erosion Risk in the Middle Reaches of the Yellow River Basin in China from 1978 to 2010. Journal of Geographical Sciences, 28, 175-192.