Assessment of Soil Erosion at Multiple Spatial Scales Following Land Use Changes in 1980–2017 in the Black Soil Region, (NE) China

Impact of land use and land cover change on soil erosion is still imperfectly understood, especially in northeastern China where severe soil erosion has occurred since the 1950s. It is important to identify temporal changes of soil erosion for the black soil region at different spatial scales. In the present study, potential soil erosion in northeastern China was estimated based on the Revised Universal Loss Equation by integrating satellite images, and the variability of soil erosion at different spatial scales following land use changes in 1980, 1990, 2000, 2010, and 2017 was analyzed. The regionally spatial patterns of soil loss coincided with the topography, rainfall erosivity, soil erodibility, and use patterns, and around 45% of soil loss came from arable land. Regionally, soil erosion rates increased from 1980 to 2010 and decreased from 2010 to 2017, ranging from 3.91 to 4.45 Mg ha−1 yr−1 with an average of 4.22 Mg ha−1 yr−1 in 1980–2017. Areas with a rate of soil erosion less than 1.41 Mg ha−1 yr−1 decreased from 1980 to 2010 and increased from 2010 to 2017, and the opposite changing patterns occurred in higher erosion classes. Arable land continuously increased at the expense of forest in the high-elevation and steep-slope areas from 1980 to 2010, and decreased from 2010 to 2017, resulting in increased areas with erosion rates higher than 7.05 Mg ha−1 yr−1. At a provincial scale, Liaoning Province experienced the highest soil erosion rate of 9.43 Mg ha−1 yr−1, followed by Jilin Province, the eastern Inner Mongolia Autonomous Region, and Heilongjiang Province. At a county scale, around 75% of the counties had a soil erosion rate higher than the tolerance level. The county numbers with higher erosion rate increased in 1980–2010 and decreased in 2010–2017, resulting from the sprawl and withdrawal of arable land.


Introduction
Soil determines the biological and geochemical cycles of the earth system and provides basic services such as food, clean water, and goods for humankind. Severe soil erosion can thus influence soil fertility, sustainable agriculture, and the functions of soil systems [1,2]. In recent decades, the rapid expansion of arable land at the expense of natural vegetation induces severe soil erosion and environmental problems such as land degradation, river and reservoir sedimentation, and freshwater pollution. A deep understanding of land use and land cover (LULC) change-induced spatiotemporal patterns of soil erosion is of great importance to land use management and sustainable development [2]. widely used at large scales such as at regional and/or global scales [7]. Furthermore, in the study area, it has been verified that it can yield a reasonably reliable soil loss estimation [24].
Therefore, although some work has been done in NEC using different methods, soil erosion characteristics of NEC with LULC changes are scarcely studied. The RUSLE model was used in the present study by integrating long-term RS images to (i) examine the temporal variations of soil erosion at regional, provincial, and county scales following LULC changes from 1980 to 2017, (ii) analyze the impacts of LULC changes on soil erosion during the past decades.

Study Area
NEC (38 • 42 -53 • 36", 115 • 24 -135 • 12") includes three provinces (i.e., Heilongjiang, Jilin, and Liaoning Provinces) and the eastern Inner Mongolia Autonomous Region (IMA) (Figure 1a). It covers an area of 1.24 million km 2 , with elevations ranging from −10 to 2666 m a.s.l. NEC is composed of three mountains, three plains, and three large rivers. The Da Hinggang Mountain, the Xiao Hinggang Mountain, and the Changbai Mountain lie in the northwest, northeast, and southwest, respectively. Three plains including Liao River Plain, Songnen Plain, and Sanjiang Plain are situated between the mountains. They were formed by the piedmont alluvial and flood sedimentation from Heilong River, Songhua River, and Tumen River that join at the east of Sanjiang Plain. There are 2300 rivers and many lakes scattered in the region. The total amount of water resource is around 199 billion m 3 , with 170.4 m 3 on the ground and 68.1 billion m 3 underground. Due to active industrial and agricultural activities as well as increasing population since the 1950s, water pollution and deficiency are severe [25]. The slopes range from zero in the plains to above 60 • in the mountainous regions ( Figure 1b).
The average annual total precipitation in 1970-2018 ranged from around 993 mm in the southeast to less than 242 mm in the northwest (Figure 1c), around 65% of which fell in summer (i.e., June-August). The mean monthly temperature varied greatly, with the highest mean monthly temperature (18 • C) in July, and the lowest temperature (−20 • C) in January.
Affected by the summer monsoon from the southeast and obstruction by the Chanbai Mountain, the vegetation in NEC is distributed zonally, from forest to the east of the Changbai Mountain, to a transition zone between forest and grasslands in the Songnen Plain, and to grassland to the west of the Da Hinggang Mountain (Figure 1a). Coniferous forest is mainly distributed in the mountainous areas, and broad-leaved forests are scattered in the plains. Arable land in the plains covers about one-third of the total area. The main crops are maize, soybean, and rice, occupying around 87% of the total. A single-tillage operation is typically used with a cultivation ploughing depth of around 0.25 m in late September and in early May. From October to April, arable land is left fallow with no vegetation cover due to cold weather.
The main soil types are Luvisols, Phaeozems, Chernozems, and Gleysols, occupying 37.5%, 23.3%, 7.4%, and 6.8% of the total area in NEC (Figure 1d). The parent materials are Quaternary lacustrine and fluvial sand beds or loess sediments [26]. Fuvisols are widely distributed in mountainous areas, while Phaeozems and Chaeozems, which are called black soils by the local people, mainly appear in the Songnen Plain [9]. The textural classes of the topsoil are silt-clay-loam to clay-loam (8-27% sand, 29-66% silt, and 26-40% clay), with soil organic carbon usually above 4%. The hilly regions with long slope in the mountain-plain transition zone are prone to soil erosion. Rill erosion dominates the study area. However, gullies also widely occur [16].
In the study area, agricultural lands occupy around 30% of the provincial area for Heilongjiang, Jilin, and Liaoning Provinces. The major crops are rice, corn, and soybean. In Heilongjiang Province, their area percentages are 0.26, 0.44, and 0.21. In Jilin and Liaoning Provinces, their area percentages are 0.13, 0.70, 0.03 and 0.11, 0.66, and 0.02, respectively (Table 1).

Data Collection
Several data layers are required to run the RUSLE model: a Digital Elevation Model (DEM), land use maps, a rainfall erosivity map, a soil erodibility map, and land use data.
The DEM was prepared from Shuttle Radar Topography Mission (SRTM) with 90 m resolution. It was provided by Geospatial Data Cloud site, Computer Network Information Center, Chinese Academy of Sciences [27].
The Harmonized World Soil Database (HWSD) Version 1.2 with 30 arcsecond resolution was obtained from Food and Agriculture Organization of the United Nations [28]. Soil properties such as percentages of sand, clay, and silt; depth of soil; and organic carbon are available from the HWSD.
The daily precipitation data in 1970-2018 at the 120 meteorological stations were acquired from National Climate Centre of the China Meteorological Administration [29].
The land use datasets with 30 m resolution were acquired from different sources. The datasets in 1980, 1990, and 2000 were obtained from the Resource and Environmental Sciences, Chinese Academy of Sciences (RES-CAS) with an overall accuracy of 92.9% [18], which has a hierarchical classification system of 25 land cover classes. The datasets in 2010 and 2017 were downloaded from the GlobeLand30 [30] and from the research team of Gong Peng in Tsinghua University [31]. The GlobeLand30 product was generated from Landsat TM/STE+ and the Chinese BJ-1 and HJ-1sensors. This product consists of nine land cover types including barren, grass, cultivated, forest, shrub, water, artificial, wetland, and ice. Its overall accuracy is around 80% in 2010 [32]. The land use data by Gong Peng has a similar land use classification system, including cropland, forest, grassland, shrubland, wetland, water, tundra, impervious, barren, and snow/ice, with an overall accuracy of 72.76% [33]. These datasets were reclassified into nine land use types ( Table 2). Table 2. Area

Description of the RUSLE and Data Treatment
The RUSLE was employed to estimate the soil erosion rate. The RUSLE is a linear equation used to estimate soil erosion by water from hillslopes. It is able to reduce a complex system to a quite simple one to predict long-term soil erosion rates over large areas while maintaining a through representation of the main factors that influence the process [7]. It combines five factors, including rainfall erosivity, soil erodibility, slope length and gradient, land use management, and soil conservation practices: where A is the soil erosion rate (Mg ha −1 yr −1 ), R is the rainfall-runoff erosivity factor (MJ mm ha −1 h −1 ), K is the soil erodibility factor (t ha h ha −1 MJ −1 mm −1 ), LS is a combination of slope length L and slope gradient S factor (-), C is the crop/cover management factor (-), and P is the soil conservation factor (-).
The classic method to calculate R factor is to multiply rainfall energy and the maximum 30 min rainfall intensity for rainfall events [34]. However, such datasets are not available for many places in the world. In such instances, alternative methods have been proposed to calculate R factor using readily available precipitation data [35,36]. In the current study, the daily precipitation data in 1970-2018 were used to estimate mean annual R-factor values, using the method proposed by Zhang et al. [36]: where R i is the half-month R-factor (MJ mm ha −1 h −1 yr −1 ), and D j is the erosive rainfall on day j. D j is equal to actual rainfall when it is greater than 12 mm. Otherwise it equal to zero. Here, a and b are empirical parameters. A detailed description was given by Fang and Sun [10]. A co-kriging interpolation was used to derive an R-factor value map in ArcGIS 10.5 software (ESRI, Redlands, CA, USA). In NEC, the mean annual R-factor values ranged from 475 MJ mm ha −1 h −1 yr −1 in western NEC to 5320 MJ mm ha −1 h −1 yr −1 in southeastern NEC ( Figure 2a). Soil texture and soil organic matter of the top 30 cm of soil were available from the HWSD, and the K factor values were calculated by using the EPIC (Environmental Policy Integrated Climate) model method. This method has been successfully used in NEC [10,17]: where SAN is the sand content (%), SIL is the silt content (%), CLA is the clay content (%), SC is the soil organic carbon content (%), and SN = 1 − SAN/100. The K factor values of the study region varied from 0 to 0.045 t ha −1 h MJ −1 mm −1 ha −1 (Figure 2b). The L factor measures the slope length, and the S factor is proportional to the local slope. This study used the 2D approach to calculate the LS factor as proposed by Van Rompaey et al. [37], using the methods by McCool et al. [38] and Desmet and Govers [39]: where L i,j is the slope length factor of the grid cell with coordinates i and j, A i,j is the contribution area at the inlet of a pixel (m 2 ), D is the grid cell size (m); α i,j is the slope aspect direction for a pixel; m is slope length exponent, β is empirical factor (-), and θ is the slope angle. The derived LS factor values ranged from 0.03 to 70.02. Higher values occurred in the mountainous regions and lower values in the plains (Figure 2c). The C and P factors are closely linked to land use types [40]. In order to obtain reasonable C-and P-factor values, a detailed survey of the literatures published in NEC was first carried out [19,24,[41][42][43]. Then, the C-factor values in forest, grassland, shrubland, cultivated land, and residential areas were assigned as 0.004, 0.043, 0.07, 0.228, and 0.03, and those in bare land, wetland and water bodies, and other unused lands were 1, 0, and 0.06, respectively. P-factor values in cultivated land and water body were set to 0.352 and 0, and in other lands they were set to 1 in the study area.
The estimated soil erosion rates were further classified into five erosion classes from low (0T), moderate (1-5T), high (5-10T), and very high (10-20T), to critical (>20T) where T represents soil loss tolerance of 1.41 Mg ha −1 yr −1 in NEC [6]. In order to verify the estimated results, some published results from runoff plots, natural slopes, and catchments were also used in the present study.

LULC Changes in NEC
The land use types in 1980-2017 were shown in Figure 3 and numerically in Table 2. Spatially, forest was distributed in mountainous areas, and arable lands in the three plains ( Figure 1). Among the land use classes, forest was predominant and covered 37% of the whole area, followed by arable land and grasslands covering 30% and 22%, respectively. These three kinds of lands occupied around 90% of the total. Shrubland, wetland, water body, and unused lands each covered 2% of the total area.
During the past decades, great changes in arable and forest lands occurred. Before 2010, arable land percentages increased from 27% in 1980, 29% in 1990, and 30% in 2000 to 35% in 2010, and then decreased to 32% in 2017, fluctuating from 41.34 million km 2 in 1980 to 45.63 million km 2 in 2010. In contrast, forest area percentages decreased from 38% in 1980, 36% in 1990, and 37% in 2000 to 34% in 2010, and then increased to 37% in 2017. The area percentages of grassland decreased too before 2000, ranging from 21% in 1980 and 20% in 1990 to 19% in 2000, and then increased sharply thereafter, occupying around 25% in 2010-2017. Shrubland percentages occupied around 4% before 2000 and less than 0.1% thereafter. The percentages of residential area in general increased from 1% in 1980 to 2% in 2017. On the contrary, the area percentages of wetland and water body decreased from 3% and 3% in 1980 to 0% and 1% in 2017, respectively. Unused lands kept unchanged through time. Since 2010, bare land has almost disappeared. Noticeably, great LULC change occurred in 2000-2010 (i.e., P3 in Table 2), except for those of the wetland and bare land in 2010-2017. Arable land also changed with time on different slopes. In areas with elevation above 200 m a.s.l and slopes above 5 degree, the distribution of arable land showed similar change patterns in different years, i.e., the areas of arable land increased from 1980 to 2010 and decreased in 2017. In contrast, in the low-elevation and gentle-slope regions, arable land area increased continuously (Figure 4c,d).

Changes in Soil Erosion
In 1980-2017, spatial patterns of annual soil erosion rates changed little ( Figure 5). Annual soil erosion rates in NEC ranged from 0 in the plains to over 40 Mg ha −1 yr −1 in the mountainous areas, with an average of 4.22 Mg ha −1 yr −1 (Figure 6a). It was three times the soil loss tolerance (T). Annual soil loss reached around 5.24 million tons from NEC. Figure 6b indicates that 57% of NEC belonged to the low erosion class, followed by moderate erosion with area covering around 30%. The high, very high, and extreme erosion areas covered 6%, 3.5%, and 2.7% of the total, respectively. Land use greatly influences soil erosion. The rates of soil erosion on bare land, shrubland, and arable land were much higher than those on other land use. Soil loss from arable land changed from 40% in 1980 and 46% in 2010 to 34% in 2017 of the total in respective years with an average of 40% (Table 2). The area percentages with different erosion classes differed with time ( Figure 6). The area percentages in the low erosion class decreased from 59% in 1980 to 56% in 2010 and increased to 57% in 2017. Opposite changing patterns occurred for the area percentages with high, very high, and extreme erosion severities. The threshold values occurred in 2010. In contrast, the area percentages with moderate erosion severity increased continuously with time.
The rates of annual soil erosion presented an obvious spatial heterogeneity among the provinces and the IMA (Figure 6). The mean soil erosion rate was the highest in Liaoning Province (i.e., 9.43 Mg ha −1 yr −1 ) in 1980-2017, followed by Jilin Province, and the IMA. The lowest erosion rate occurred in Heilongjiang Province, with an average of 2.46 Mg ha −1 yr −1 (Figure 6c). Considering the area percentages, their annual soil losses accounted for 21%, 17%, 26%, and 35% of the total area in NEC, respectively.  The spatiotemporal characteristics of soil erosion severity were mapped for the 224 counties in NEC (Figure 8). In 1980-2017, the number of the counties with soil erosion rates higher than one T accounted for around 75% of the total. The counties with soil erosion rates higher than 5T and 10T accounted for 24% and 8% of the total, respectively ( Figure 9). For the top 16 counties with erosion rate higher than 10T, 12 and 4 counties were distributed in the southeastern and southwestern NEC, respectively. Almost all the counties with soil erosion rate lower than one T were distributed in the plains (Figure 8).    (Figure 9a). The number percentages of the counties with low, moderate, high, very high, and extreme zones occupied 23%, 53%, 16%, 6%, and 1%, respectively.
In 1980-2017, the counties with low erosion rate ranged from 20% to 25% of the total and remained unchanged after 2000. The counties in moderate soil erosion zone decreased from 54% in 1980 to 49% in 2017, among which the counties with soil erosion rate 1-2T decreased fastest, from 21% in 1980 to 17% in 2017. Noticeably, the numbers of counties with high and very high classes increased from 14% and 5% in 1980 to 18% and 8%, respectively (Figure 9b).

Model Performance
In NEC, soil erosion rates have been identified using runoff plot [14], 137 Cs tracer technique [15,44,45], and model estimation [10,42]. Comparisons of the RUSLE-derived rates of soil erosion in the present study to published results can evaluate the performance of the RUSLE.
At plot and slope scales, the RUSLE-derived rates of soil erosion in the current study were less than the reported results (Table 3). This difference could result from two aspects. One aspect is that the 90 m resolution topographic data can smooth slope gradient [6,46], resulting in lower soil erosion rates. The other aspect is the difference in time spans between the studies in literature and the current study. The rates of soil erosion have also been reported at larger scales, including at a 28.5 ha catchment in Heshan Farmland [15], at a 916 km 2 Shuangyang river catchment [10], and in the Sanjiang Plain [42]. A comparison at a catchment scale is more meaningful since this study was done at a large spatial scale. Table 3 demonstrates that the RUSLE-derived rates of soil erosion in the current study are basically consistent with the published results. Table 3. The published rates of soil erosion in the black soil region (the study area locations of the references were labeled in Figure 1). The average bulk density 1.3 t m −3 was used when the unit mm yr −1 was converted to unit Mg ha −1 yr −1 according to Xie et al. [9]; CGWCR indicates the Complication Group of Water Conservancy Record in Baiquan County; UPC represents up-down cultivated land. The first numbers in the column "Region" represent the counties labeled in Figure 1. * indicates that the locations where the erosion rates were measured were given in Figure 1.

Impact of LULC Change on Soil Erosion
The changes of soil erosion from 1980 to 2017 can be attributed to the LULC with time, due to population pressure and national agricultural policies [5,11,48]. Since the mid-20th century, extensive settlement and land exploitation began in NEC, and a large amount of woodlands were converted into arable land [11]. Because of reduced precipitation and more water use consumption by increased crop land, around 4% forest and 3% wetland were lost, and around 3381 km 2 woodland was converted to arable land in 1990-2000 [12]. As a result, soil loss from arable land accounted for 42% of the total, followed by grassland, forest land, and other land uses (Table 4). NEC is the most important Chinese grain production base, and a series of policies further stimulated arable land reclamation in 2000-2010. For example, the rescission of the tax from agriculture in 2005 and the household responsibility system triggered more cultivated farmland at the expense of woodland (i.e., forest and shrublands), resulting in a 46% soil loss of the total from arable land in 2010 (Table 4). Comparable higher soil loss rates on arable land were also reported in other regions of the world. In the U.K., a total of 1566 individual records indicated that arable land had the highest mean erosion values among different land use types [49]. An average soil erosion rate of 4.04 Mg ha −1 yr −1 was also observed in the agricultural lands in the European territory of Russia [50]. For the part of Europe covered by the CORINE database, soil loss from arable land is 3.6 Mg ha −1 , which is three times of the average rate [51]. Soil loss from the agricultural land in Europe accounts for 93.5% of the total in Europe. A review study by Guo et al. [52] also gave higher soil erosion rates of arable lands in Australia (11.0 Mg ha −1 yr −1 ), Africa (14.9 Mg ha −1 yr −1 ), and USA (18.378 Mg ha −1 yr −1 ). Higher soil loss rate on arable land in the study area and its area can explain that arable land was the major sediment source in NEC. To decrease soil loss, the "Grain to Green Program" (GTGP) was introduced. This program is the largest ecological restoration in the developing world, initiated from 1999 and has made a remarkable contribution to China's vegetation recovery [53]. Similar programs also appeared in other countries. For example, the Great Green Wall of the Sahara and the Sahel (GGW), which covers 21 countries in Africa [54], the First Ten-year Forest Rehabilitation (TYFR) in South Korea, and the Five Million Hectare Reforestation Program (5MHRP) in Vietnam. In NEC, impacted by the GTGP, farmland shelterbelts and sand fixation forests were constructed, and conversion of farmland to forest occurred on steep slopes [11,12,36]. From 2010 to 2017, forest area increased by 39,928 km 2 , with the largest change ratio of 9% during 2010-2017 (Table 4). The changes in land use policies and associated arable land cultivation can explain the increased area percentages with low erosion class in 1980-2010 and the decreased percentages from 2010 to 2017. Impacted by the policies, continuous sprawl of arable land in 1980-2010 and its withdrawal in 2010-2017 occurred in the mountainous areas.

Impact of Topography on Soil Erosion
The spatial distribution of soil erosion rates and their changes with time can be explained by the topography in NEC. First, steep mountainous areas usually have higher LS values [55]. This distribution characteristics has been found in literature [1,56,57]. Second, more rainfall and higher rainfall erosivity occurred in the mountainous areas ( Figure 2a) and lower values occurred in the plains (Figure 3). The spatial pattern of the R value in mountainous regions can also be found in Sudetes Mts., SW Poland [1]. In NEC, soil erosion rates increased with increasing slope gradients (Figure 4a,b). The rate of soil erosion above 300 m a.s.l. was four to five times that in areas below 200 m a.s.l. Due to steep slope and high R values, the rates of soil erosion were high in high-elevation and steep mountainous regions, and low in plains (Figure 4a,b). Therefore, the spatial distributions of soil erosion rates at regional and county scales were quite related to topographic factors (Figures 1 and 6), which is consistent with the studies by Bakker et al. [58] and Latocha et al. [59]. This can explain the spatial distribution of soil erosion rates at regional and county scales in NEC (Figures 4 and 6).

Study Limitations
The soil erosion rate was calculated through multiplying the five factors (i.e., R, K, LS, C, and P). The data accuracies, especially the values of the C and P factors, greatly influence the results. Without consideration, the varied C-and P-factor values with time could result in biased estimation [6]. Furthermore, the datasets with 90 m resolution could influence estimation accuracy at plot or slope scales. However, although finer resolution could improve estimation accuracy, higher computational capacity is required, and it could not change the variation patterns of potential soil erosion.
Except for data resolution, the preclusion of gully erosion from the RUSLE can also induce biased estimation. In NEC, gullies are widely distributed, and over 295,000 (ephemeral) gullies were found in this region [60]. The soil loss from gullies is high. In a reservoir catchment in Baiquan County, the sediment from gullies occupied around 40% of the total [61]. Similar result was also obtained by Huang et al. [62] in Heshan Farmland, NEC. At a catchment scale, models such as WaTEM/SEDEM (Water and Tillage Erosion Model/Sediment Delivery Model) [37] and TeTIS [16] have combined the soil erosion module and the sediment deposition module to simulate soil sediment yield. The combination of such modules with RUSLE could improve the estimation accuracy.
The uncertainty also comes from the applications of LS calculation equations and land use data sources. This kind of limitation has fully been explained by Alewell et al. [63], and thus detailed explanation was not given here. However, it should be stated that the application of the upslope area-based L factor calculation could improve the estimation accuracy. The differences in land use classification systems and data resolution from RES-CAS, GlobeLand30, and Gong Peng research group in Tsinghua University can also lead to some bias of the estimated soil erosion.
Uncertainty also came from the employed model in the present study. The origin version of RUSLE is based on a field plot. With increasing GIS-based computational technique, it can also be applied to regional and global scales, and some studies indicated there are substantial differences between the results from RUSLE and the measurements from microrelief [64,65]. In contrast, Borrelli et al. [7] also pointed out the difference is not substantial. Although the estimated soil erosion rates were potential and approximated ones, they can reflect the variations of soil erosion and its responses to LULC changes in the study region, the focus of the present study.

Conclusions
Assessment of the impact of regional LULC change on soil loss is of paramount importance to appropriately regulate land use. Soil loss characteristics and its changing patterns following LULC changes in 1980-2017 were explored in NEC using the RUSLE method at regional, provincial, and county scales, and some findings were obtained.
In NEC, higher rates of soil erosion occurred in the mountainous regions and lower values in the plains. In the study period, agricultural policies greatly influenced LULC that determine temporal changes in soil erosion rate at different scales. Soil erosion rates increased from 1980 to 2010 and decreased from 2010 to 2017. Spatially, the highest soil erosion rate of 9.43 Mg ha −1 yr −1 occurred in Liaoning Province and the lowest value in Heilongjiang Province. At a county scale, the counties in the southeast and southwest of NEC experienced the highest soil loss rate. In recent years, around 75% of the counties still suffered soil loss with soil erosion rates higher than the tolerance level. Arable land was the major sediment source, appropriate policies should be created to control soil loss through limiting the sprawl of arable lands in unfavorable regions, and land use management should also be done for the bare-, grass-and shrublands because these lands also suffered higher soil erosion rates. Because this study just described soil loss patterns at a large scale based on RUSLE, further study on the efficiencies of soil conservation practices at finer scales should be done in this area.
Author Contributions: H.F. collected information and wrote the paper. Z.F. contributed to statistical analysis and revised the manuscript. All authors have read and agree to the published version of the manuscript.
Funding: This work was financially supported by projects of the National Natural Science Foundation of China (grant number 41977066, 41971358) and Shandong Provincial Department of Education Project (J17KB040).

Conflicts of Interest:
The authors declare no conflict of interest.