1. Introduction
Increasing demand for land outputs encourages farmers to intensify their use of arable land [1,2,3,4]. Multiple cropping by planting and harvesting crops in the agricultural field multiple times a year is an important form of cultivation in China and other Southeast Asian countries [5,6,7,8]. Cropping intensity refers to the number of cropping cycles per unit area per year, while cropping cycles refers to the number of crops harvested per year, such as single cropping, double cropping, and triple cropping. Variations in cropping intensity not only reflect the changes in agricultural land use practices but also reflect the transformation of agricultural land functions [9,10,11]. Understanding the landscape pattern and changing mechanism in multiple cropping could help provide a scientific reference for the planning of arable land uses [11,12,13,14].
The dynamic change of cropping intensity is influenced by a number of factors such as natural conditions, policy, and social conditions [15,16,17,18,19,20]. Generally, social, economic, and technological factors can affect regional land use changes on a short time scale, while natural factors cause changes in land cover through a relatively long process of succession. Among all the factors, human activities affect the change of cropping intensity directly. In turn, changes in cropping intensity can affect social and economic development [21,22]. Global urbanization has led to a significant abandonment of arable land in many developed countries and regions [23,24,25]. In China, urban and rural employment, food consumption by people, industrial and agricultural relations, and urban–rural relations have experienced great transformations [26,27,28]. The rising cost of labor and the reduction of the residential populations in rural areas have led to large-scale land abandonment. At the same time, mechanized production, agricultural technology promoted by life science advances, economic agriculture, and ecological agriculture transformation are positive for agricultural production [29,30]. By understanding changes in cropping intensity, it is relevant to policy to explore the process of rapid urbanization and the interactions between human activities and land cultivation. However, the lack of continuous and relatively high spatial-resolution cropping intensity products has limited our ability to understand the changes of cropping intensity and the underlying driving forces over regions.
Researchers have made tremendous efforts to identify cropping intensity over large areas using satellite data. There are practical solutions to mapping and monitoring cropping intensity based on a time series of satellite-derived vegetation indexes, and the mapping results of cropping intensity have been validated from regional to continental scales. Ding et al. applied a two-difference algorithm to map multiple cropping from 1982 to 2012 in 17 provinces in northern China using Global Inventory Monitoring and Modeling System normalized difference vegetation index (GIMMS-NDVI) data at an 8 -km resolution and validated the produced maps of the cropping index at the province level with agricultural statistical data [31]. Li et al. developed an iterative moving-window method and produced cropping cycle maps in mainland China during 2007–2012 at a 500 m resolution based on moderate resolution imaging spectroradiometer (MODIS) enhanced vegetation index (EVI) time series and temperature data. The produced maps were validated based on interpreted samples at the pixel level and compared with agricultural statistical data at the provincial level [6]. Li et al. identified cropping intensity at a 30 m resolution using a fusion of Landsat and MODIS data in two scenes and evaluated the maps with results derived from very high resolution images [32]. Estel et al. produced cropping frequency, multi-cropping, fallow, and crop duration maps in Europe from 2000 to 2012 using the MODIS normalized difference vegetation index (NDVI) at a 500 m resolution by classifying the cultivated lands into active or fallow status and then counting the peaks during the growing seasons [33]. Gray et al. mapped cropping intensity across Asia for the time period of 2009–2012 with MODIS network-based application recognition NBAR data by counting the total number of valid peaks. The produced 500 m maps were evaluated with inventory data at district, provincial, and national levels and were also evaluated with field survey data at the pixel level [5]. Liu et al. proposed a shape-matching method and mapped the cropping index in a province in China for a single year. The map was validated at the pixel level with field survey data [34]. Qiu used a wavelet feature-based method to extract cropping intensity from wavelet spectra and derived the cropping intensity maps from 1982 to 2013 using both advanced very high-resolution radiometer (AVHRR) and MODIS vegetation index data. These maps were validated with ground truth data at the pixel level and validated with agricultural census data at the provincial level; They were then used to analyze the influences from socioeconomic factors, topographic factors, and water condition factors at the provincial level by a multiple linear regression model [35]. Although many attempts focused on the development and application of extraction methods at regional scales, none of the cropping intensity datasets have been publicly distributed, and there are still large uncertainties associated with cropping intensity changes over years.
Understanding the mechanism underlying cropping intensity (and its changes) is key to the development of a simulation model for arable land uses. There are Empirical statistical methods and modelling methods needed to explore the driving forces behind land use. Empirical statistical methods, such as multiple linear regression methods and principal component analysis, explore the relationship between land use change types and their driving factors by establishing statistical functions [35,36,37]. The drawback of an empirical statistical method is that the spatial relationship between land use changes and driving factors is not considered. The modelling method aims to construct a complex model that is able to simulate the relationships among the structure, function, and dynamic changes across the entire land use system [38,39,40]. Researchers often use modeling methods to simulate land use changes; Identifying the key driving factors and understanding the mechanisms underlying land use changes are key to successful modeling approaches. The Geo-Detector is a statistical method that was first proposed by Wang et al. [41,42]. A Geo-Detector analyzes the spatial stratified heterogeneity (i.e., the values of an attribute that vary across types and geographic regions) of factors and obtains the importance of determinants on the land use changes. Compared with commonly used empirical statistical methods, the Geo-Detector accounts for the geographic relationships between land use changes and influential factors. Compared with modelling approaches, the Geo-Detector does not require complicated assumptions and is simple and convenient. Geo-Detectors have been used in a number of studies on exploring the mechanisms of socioeconomic factors and natural environmental factors on land use changes and have proven suitable for analyzing complex interactions between socioeconomic driving forces and land use changes over large geographic areas [43,44,45,46].
The extent of Chinese agricultural land covers large areas, and the factors that affect cropping intensity could be distinct in different agricultural zones. To date, there are few studies that have explored the socioeconomic driving forces underlying cropping intensity and their changes over multiple years across large geographic areas. Therefore, the objectives of this study were to quantify China’s cropping intensity from 2001 to 2016 and to evaluate its associated socioeconomic driving forces in different regions across China to derive an extended understanding of the impact of human activities on farmland ecosystems.
2. Study Area and Materials
The agricultural spatial pattern in China is characterized by changes in climate and land resources. Figure 1 shows the nine agricultural zones in China with a backdrop of the Chinese land cover map. Considering the geographical differences in light, heat, water, and soil resources that are related to climate and terrain, the agricultural patterns in China were divided into Northeast China Plain (NE Plain), Yunnan–Guizhou Plateau (YG Plain), Northern arid and semiarid region (NAS Region), Southern China (SC), Sichuan Basin and surrounding regions (SBS Region), Middle-lower Yangtze Plain (MLY Plain), Qinghai Tibet Plateau (QT Plateau), Loess Plateau (L Plateau), and Huang-Huai-Hai Plain (HHH Plain) [47]. The density of population and human activities in these nine regions are different, and a combination of socioeconomic factors and natural geographical backgrounds forms different spatially heterogeneity regions. Given the differences in spatial distribution between socioeconomic factors and agricultural natural zones, this study detected socioeconomic drivers by analyzing the socioeconomic interactions in different regions.
The MODIS surface reflectance 8 day L3 Global 500 m SIN Grid products, based on both the Terra (MOD09A1) and Aqua (MYD09A1) platforms, were download from 2000 to 2017 from the USGS website (
The elevation, aspect, and slope data related to the arable land topography were derived from the NASA Shuttle Radar Topographic Mission (SRTM) digital elevation data at a 90 m spatial resolution (
The agricultural economic census data were collected from the China County Statistical Yearbook 2017 (County and City Volume) (
3. Methods
3.1. Cropping Intensity Quantification
The cropping intensity products from 2001 to 2016 were generated based on our published work (Li et al., 2014) [6] with a refinement of sliding windows and land cover maps. First, the EVI time series were derived from the surface reflectance data from both MODIS Terra and Aqua data, which were then further synthesized based on the minimum blue compositing method [32]. Second, considering the effects of the three fitting methods, including asymmetric Gaussian function, double logistic, and the Savitzky–Golay, the EVI time series data were filtered to obtain a smoothed EVI time series using the SG filter in the software TIMESAT to fill in the missing values and remove noise values and outliers. Then, the peaks and toughs in the smoothed EVI time series were extracted by sliding windows and used to derive the cropping cycles maps at 500-m resolution by counting the peaks and toughs. The cropping cycles included fallow, single cropping, double cropping, and triple cropping, by counting the peaks and toughs. The methodological details for the cropping cycle mapping can be found in Li et al., 2014. Here, the MODIS land surface temperature data were used to eliminate false peaks caused by the sprouting of winter wheat before dormancy in winter. During the night, a land surface temperature for the MODIS products (MOD11A2) with minimum threshold of 5 °C was used to define the real growing season in a year and to deduct the false peak caused by the winter wheat. Moreover, GlobeLand30 was used to refine the cropping intensity maps by defining the class of the cultivated lands. Finally, based on the cropping cycle maps, the cropping intensity was obtained by calculating the ratio of the total sown area of cultivated land to the cultivated land area at a 5 km resolution, while the total sown area was obtained by multiplying the cultivated land area by cropping cycles. When the cropping intensity was between 0 and 1, the cropping cycles were mainly single cropping. When the cropping intensity was between 1 and 2, the cropping cycles are mainly double cropping mixed with single cropping. When the cropping intensity was between 2 and 3, the cropping cycles were triple cropping mixed with double cropping.
3.2. Detecting Socioeconomic Drivers Using the Geo-Detector
The Geo-Detector includes four components: The factor detector, the interaction detector, the risk detector, and the ecological detector. In this study, we mainly used the factor detector to explore the influence of independent variables (i.e., potentially influential factors) on dependent variables (i.e., cropping intensity and its changes). The main idea of the factor detector in the Geo-Detector is to classify geographic units (e.g., pixels or administrative units such as counties, cities, and provinces) based on each individual influential factor using a clustering method and to derive the variance of dependent variables (e.g., cropping intensity or its changes) in different classes of geographic units to obtain the power of determinant (PD) that represent the relative importance of determinants.
The PD in the Geo-Detector is quantified as follows [42]:
where PD denotes the power of determinant; m denotes the number of the stratified classes, donates the overall variance of the dependent variable, donates the variance of the dependent variable in each class, denotes the number of geographic units, and denotes the number of geographic units within each class. The right side of the equation was solved by an F-test based on the assumption that the spatial autocorrelation coefficient is second order spatially stationary. More details about the Geo-Detector algorithm can be found in (Wang et al., 2010) and (Wang et al., 2016) [42,49].The range of the PD value was [0, 1]. The larger the PD value was, the greater the importance of the influential factor on dependent variables. If PD was equal to 0, the dependent variables (e.g., cropping intensity and its changes) are distributed randomly. If PD equals 1, the influential factor completely determined the spatial distribution of the dependent variables.
3.3. Potential Driving Factors and Data
Based on the socioeconomic data, we derived 13 statistical indexes as potential influential factors related to crop yield, agricultural machinery power, farmers’ income, and industrial developments. The 13 factors include grain yield (GY), grain yield per labor (GYlabour), arable land area per capita (ALAcaptita), agricultural machinery power per unit of arable land area (AMPALA), percentage of machinery harvest area in arable land area (PMHAALA), labor force per unit of arable land area (LFALA), terrain elevation (TE), terrain slope (TS), gross product of industrial enterprises above a designated size (GPIE), regional GDP(GDP), the growth rate of the percentage share of the agricultural sector in GDP (GRAG), the growth rate of the percentage share of primary industry sector in GDP (GRPIG), and farmers’ income per capita (FIcapita).
TE was directly extracted from the SRTM data, and TS was calculated from the SRTM data. Both TE and TS were averaged to the county level. The GPIE and GDP were obtained from the agricultural economic census data. The formulas used to derive the other indicators are in Table 1.
In this study, we conducted an analysis on the county level and used the average annual cropping intensity and/or its changes during 2001–2016 as the dependent variables. We classified the counties into three clusters for each individual potential factor as inputs to the Geo-Detector using the K-means clustering method. For example, we selected the GY values of three counties as the seed points and then clustered other GY values according to the shortest Euclidean distance between the GY value to be clustered and the seed GY. The cluster centers were adjusted by the minimum Euclidean distance in each cluster and the maximum Euclidean distance between clusters. After 100 iterations, the GY values were divided into three categories. The counties then were divided into three categories according their GY values’ clustering results.
4. Results
4.1. The Spatiotemporal Patterns of the Cropping Intensity
Figure 2 shows the cropping intensity in China for each individual year during 2001–2016, as derived from MODIS time series data at a 5-km resolution. The spatial distribution of the mapped cropping intensity showed topographic and climatic differentiation patterns. As topography changed from mountain to plain, the cropping intensity increased gradually from west to east. As hydrothermal conditions improved from north to south, the cropping intensity increased gradually from north to south.
The annual mean cropping intensity was summarized in each agriculture zone for 2001–2016 (Figure 3). The SC, MLY Plain, and HHH Plain are the top three regions with the highest cropping intensity. Table 2 illustrates the cropping intensity trends for the nine agricultural zones during 2001–2016, including the annual mean cropping intensity during 2001–2016, trend slope, trend intercept, Mann–Kendall Z-value, Mann–Kendall p-value, and Mann–Kendall significance (0.05). The results indicate that all cropping intensities had a basically constant or slightly increasing trend in the nine regions. SC, L Plateau, and YG Plateau are the top three regions with the largest cropping intensity growth rate.
As shown in Table 2, the main cropping form in NE Plain was single cropping per year, and its annual cropping intensity was 97.8%, on average, from 2001 to 2016. This area is located in a temperate and warm temperate zone, with a short warm summer and a long cold winter. Due to temperature limitations, the hydrothermal conditions are not sufficient for multiple cropping. The fertile agricultural land and the relatively flat terrain are suitable for large-scale machinery planting and harvesting. The cropping intensity in the NE Plain was close to 100% with contiguous cultivated land and did not change significantly.
In NAS regions, single cropping was the main cropping type, with an average annual cropping intensity of 61.3% from 2001 to 2016. Due to water scarcity, this region is not very suitable for agriculture; Its arable land is fragmented, and its cropping intensity is mostly less than 100%.
Many of the agricultural practices in the HHH Plain are double cropping, and the average annual cropping intensity was 135.9% in 2001–2016. The area has a warm temperate monsoon climate with identical seasons. Due to fertile soil and abundant river resources, the HHH Plain is one of the main rainfed agriculture-producing areas. In the north of the Yellow River, there was single cropping every year or triple cropping every two years, while there are mainly double croppings in the south of the Yellow River. Overall, the cropping intensity increased slightly in the HHH Plain.
The L Plateau mainly had single cropping and an average annual cropping intensity of 85.7% in 2001–2016. The cultivated lands are largely rainfed and are located in hilly areas with considerable soil erosion and drought conditions.
The average annual cropping intensity of the QT Plateau was only 41.8%. The area is not suitable for multiple cropping because it has strong radiation, a long photoperiod, a low temperature, a low accumulated temperature, and a large daily temperature range. In addition, this region is the largest plateau in China, and its proportion of suitable agricultural land is small. The cropping intensity did not change significantly in 2001–2016.
The main planting cropping cycles in the MLY Plain were double cropping with little triple cropping. Low terrain, fertile soil, and abundant river resources allow for multiple cropping. As latitude decreases, the agricultural land is small and scattered, resulting in a cropping intensity less than 200%. The average annual cropping intensity was 138.0%, and the cropping intensity increased slightly from 2001 to 2016.
The SBS regions were mainly double cropping with an average annual cropping intensity of 97.3%. A warm and humid climate is helpful for agricultural development, but due to topographical features, the arable lands are fragmental. Cropping intensity growth increased significantly from 2001 to 2016.
Triple cropping mixed with double cropping was dominant in SC, where the average annual cropping intensity was 171.4%. A warm and humid subtropical monsoon climate and sufficient hydrothermal conditions are all prerequisites for multiple growing seasons in one year. This area features high-speed urbanization and industrialization and the fragmentation of cultivated land. The cropping intensity in SC was close to 200% and increased in 2001–2016.
Similar to SC, there was mainly double and triple cropping in the YG Plateau, with an average annual cropping intensity of 128.0%. The cultivated land in this area is mainly distributed in mountain basins, river valleys, and mountainous areas. Although light and water are abundant, and are affected by the topography, the cropping intensity was much lower than 200%. The cropping intensity also increased in 2001–2016.
4.2. Changes in Cropping Intensity
National cropping intensity showed a growing trend in 2001–2016 (Figure 4). More than 81.0% of the arable lands increased their cropping intensity from 2001 to 2016, and approximately 12.2% of cultivated lands increased their cropping intensity by more than 10%. Areas with decreased cropping intensity are mainly located in the NE Plain. The identified changes in cropping intensity are closely related to Chinese agricultural policy in the last 20 years, such as encouraging agriculturally relevant enterprises to sign long-term contracts to rent farmland from local farmers and encouraging agricultural activities to maintain their agro-ecological environments.
Figure 4 also shows the changes in average annual cropping intensity for every four year interval from 2001 to 2016. From 2001 to 2004, 74.4% of the arable lands had increased cropping intensity, 63.3% of the arable lands had more than 10% increased cropping intensity, and 6.5% of the arable lands had more than 10% decreased cropping intensity. Areas with the most obvious growth in cropping intensity were located in the L Plateau, the adjacent belt between the NE Plain and NAS regions, and the adjacent belt between the L Plateau and NAS regions.
From 2005 to 2008, 64.0% of the arable lands had increased cropping intensity, 12.2% of the arable lands had more than 10% increased cropping intensity, and 9.4% of the arable lands had more than 10% decreased cropping intensity. Areas with the most remarkable ascending cropping intensity were in the L Plateau and the fragmental arable lands that are scattered in southern China.
From 2009 to 2012, 66.9% of the arable lands had increased cropping intensity, 16.1% of the arable lands had more than 10% increased cropping intensity, and 7.8% of the arable lands had more than 10% decreased cropping intensity. Areas with increased cropping intensity were in the L Plateau, the adjacent belt between the NE Plain and NAS regions, the southeast of the HHH Plain, and SC.
From 2013 to 2016, 62.0% of the cultivated lands had increased cropping intensity, 12.2% of the arable lands had more than 10% increased cropping intensity, and 8.8% of the arable lands had more than 10% decreased cropping intensity. The visibly increased areas were distributed in the adjacent belt between the HHH Plain and the MLY Plain, SC, the YG Plain, and the SBS regions.
Figure 5 shows the linear trend and its significance to the cropping intensity changes at the pixel scale. Cropping intensity increased in 73.5% of the pixels of the arable land areas in 2001–2016 in China, in which 27.7% of the pixels show significant increases in cropping intensity and 45.8% of the pixels have non-significant increases in cropping intensity. The areas with significantly increased cropping intensities are mainly located in the adjacent belt between the NAS region and the L Plateau, the eastern part of HHH Plain, and SC. The cropping intensity decreased in 26.45% of the pixels of the arable land areas in 2001–2016 in China, in which 3.95% of pixels show a significant reduction in cropping intensity, and 22.5% pixels show a non-significant reduction in cropping intensity. The areas with significantly reduced cropping intensities are mainly located in the western part of the MLY Plain, the western part of the HHH Plain, and a part of Xinjiang.
Table 3 shows the proportions of cropping practices for the arable lands in different agricultural zones in both 2001 From 2016. For the time period from 2001 to 2016, the cropping practices in the nine agricultural divisions did not change much overall. The most obvious areal changes occurred in the NAS regions and the L Plateau. Both single cropping and double cropping largely increased, and fallow lands reduced.
4.3. The Socioeconomic Drivers on Cropping Intensity Distribution
Table 4 shows the PD values from thirteen influential factors on the annual mean cropping intensity as detected by the Geo-Detector at the county level for each individual agricultural zone.
In relatively flat regions, factors that are related to the output efficiency of the arable lands had greater influence than other factors on the spatial distribution of cropping intensity. For example, in the NE Plain, the NAS regions, the HHH Plain, and the L Plateau, the most influential factors were ALAcaptita (0.332), GY (0.193), GY (0.551), and ALAcaptita (0.319), respectively.
By comparison, the topographic factors were the most influential factors in affecting the spatial distribution of cropping intensity in the rugged terrain area. For example, in the QT Plateau, the MLY Plain, the SBS regions, SC, and the YG Plateau, the most influential factors were TS(0.149), TS (0.22), TE (0.468), TE (0.252), and TE (0.307), respectively. Under similar climatic conditions, the topographic factors had greater influence on the distribution of regional cropping intensity.
In essence, both natural conditions and cultivated land outputs, including the efficiency of the cultivated land, input labor, and input machinery, were key to the spatial distribution of cropping intensity.
4.4. The Socioeconomic Drivers on Cropping Intensity Change
Table 5 shows how the influential power of factors on cropping intensity change as detected by the Geo-Detector at the county level for each individual agricultural zone. In the northern region, the factors related to the output efficiency of cultivated lands dominantly affect the spatial distribution of cropping intensity changes.
In the NE Plain, the most relevant factor that affected the changes of cropping intensity was ALAcaptita (0.304). The climatic conditions in the NE Plain only allow for single cropping, and because the agricultural sector plays an important role in the local economy, Therefore the higher ALAcaptita, the more diverse agricultural activities the farmers can carry out, and the more stability the cropping intensity can maintain.
In the NAS regions, TE (0.171), GY (0.126), and GYlabor (0.110) are the top three key factors that influenced cropping intensity changes. Because there are only limited areas suitable for agriculture, and the level of economic development is low in the NAS regions, farmers prefer to abandon lands with poor topography, low GY, or low GYlabour to instead focus on arable lands with high agricultural production efficiency.
In the HHH Plain, the top factors that affected cropping intensity changes were GY (0.409), PMHAALA (0.395), and AMPALA (0.293). The HHH Plain is the largest plain in the north of China and is rich in natural resources such as water, heat, and sunlight. It also includes several major grain producing provinces in northern China. Besides the output factors of the cultivated lands, agricultural investment factors play a decisive role in determining decadal changes in cropping intensity. Agricultural technology development and machinery application help to improve agricultural efficiency, alleviate the dependence of the grain industry on the labor force, and stabilize cropping intensity.
Both the L Plateau and the QT Plateau are not dependent on agricultural development, and the factors affecting cropping intensity variation were diversified. In the L Plateau, the main factors affecting cropping intensity were ALAcaptita (0.196), LFALA (0.188), and TS (0.166). The most influential factors that influenced the cropping intensity changes in the QT Plateau included PMHAALA (0.59), FIcapita (0.427), and GDP (0.275). Farmers tend to determine cropping activities by accounting for the factors that affect both arable land output and farmers’ income.
In the southern region, the factors that were related to regional economic development and agricultural economic benefits were key to influencing cropping intensity changes. The hydrothermal conditions in SC provide ideal natural conditions for multiple cropping, while fast economic development and urbanization promote laborers to leave farms to take jobs in industrial production with higher income. The agricultural sector in SC faces the risk of declination.
The MLY Plain, which has abundant water and heat for crop planting, is situated in the transition zone between double cropping and triple cropping. In the MLY Plain, the top factors that affected cropping intensity changes were GPIE (0.474), FIcapita (0.299), and ALAcaptita (0.211).
Similarly, in SC, a frontier under the China Economic Reform and Opening up policy and the most advanced economic zone in China, the top factors that affected cropping intensity changes were GPIE (0.788), GYlabor (0.643), and AMPALA (0.621). These factors are relevant to industrial and economic developments and the efficiency of agricultural output.
Both the SBS regions and the YG Plateau are located in areas with fine hydrothermal conditions and have natural conditions suitable for double cropping and even triple cropping. The top three influential factors on cropping intensity changes were FIcapita (0.709), GYlabor (0.284), and GY (0.282) for the SBS regions, and FIcapita (0.720), GPIE (0.255), and ALAcaptita (0.246) for the YG Plateau.
Based on the analysis of PD values using the Geo-Detector, the main agricultural areas possessed high PD values, suggesting that the selected factors were decisive in cropping intensity changes. For areas where the agriculture sector is small, their PD values were relatively low, indicating that changes in cropping intensity were relatively small.
5. Discussion
5.1. Influential Factors and the PD Values
Chinese cropping practices vary largely from north to south and from east to west. Even in the same agricultural zone, the cropping form and the main driving factors can vary. This study only selected the main potential factors related to socioeconomic and technological drivers to detect their impacts on cropping intensity and its changes. These factors are related to grain output, machinery application, farmers’ income per capita, and economic development. It would be of interest to further explore these environmental and policy-relevant factors when data are available and obtain a deeper understanding of the driving forces underlying cropping intensity changes.
The PD values derived by the Geo-Detector indicate that influencing factors accounted for 100 × PD% of cropping intensity distribution and cropping intensity change. In general, there are several factors that determine the cropping intensity distribution and cropping intensity changes together in an individual agricultural zone. The ranking of the PD values reflects the relative influence, while the magnitude of the PD values provides a more complete reference. The PD value varies with regions and scientific problems [41,42,49].
5.2. The Influence of the Factor Clustering Method on the Geo-Detector
The clustering method that classifies geographic units into classes might have impacts on the final outputs of the Geo-Detector. Here we compared the Geo-Detector results based on the K-means clustering method and the standard deviation clustering method. We found that the effect of the clustering method adopted in the Geo-Detector did not alter the main conclusions, although the obtained PD values of the influential factors could be different.
Table 6 shows the PD values derived from the Geo-Detector based on the standard deviation clustering method. When comparing the results in Table 5 with those in Table 4, the ranks of the influential factors on cropping intensity changes do not change much in each agricultural zone. According to the PD values, the ranks of the dominant factors did not change much for the six agricultural zones, including the NE Plain, the NAS region, the HHH Plain, the L Plateau, the SBS regions, and SC. For the QT Plateau, the top three influential factors based on the standard deviation clustering method are AMPALA(0.649), LFALA (0.649), and GRPIG (0.242). For the MLY Plain, the top three influential factors on cropping intensity changes were GRAG (0.200), PMHAALA (0.143), and FIcapita (0.135). For the YG Plateau, the top three influential factors on cropping intensity changes were TE (0.505), ALAcaptita (0.327), and GYlabour (0.320). For the agricultural zones of the QT Plateau, the MLY Plain, and the YG Plateau, the rank of the influential factors varies when using different classification methods. There are nonlinear interactions among the top influential factors for these three agricultural zones when using the interaction detector in the Geo-Detector in our preliminary tests. In essence, some influential factors could be correlated in specific geographic regions, and there is a need to further explore the processes and mechanisms for fine-scale zones.
5.3. The Interaction between Human Farming Activities and Rapid Socioeconomic Developments
Multiple cropping is an effective way to maintain and increase crop yield in the process of rapid arable land reduction due to urbanization. The rapid economic development of China during the past few decades has led to negative impacts on agricultural activities [50], largely due to the transformation of rural labor into industrial labor alongside rapid economic development and the encroachment of urban areas on arable lands. This study determined that the cropping intensity in all agricultural regions in China showed an increasing trend, and seven ofnine regions showed significant growth from 2001 to 2016. Especially in the outskirts of big cities, increasing cropping intensity is a beneficial solution to ensure crop yields and guarantee the value of agricultural land as one type of hinterland around urban agglomerations. The SC is located in a demonstration area that has experienced rapid economic growth under the policy of Economic Reform and opening up and increasing cropping intensity could provide more space for economic developments.
The awareness of environmental protection, arising alongside ecological civilization construction precipitated by rapid economic development, promotes cropping intensity. Since 2000, more local governments have paid attention to ecologically sustainable development and promoted a series of policies. The cropping intensity in the NAS Region increased slightly, as it was influenced by these national strategic policies, such as guaranteeing the area of cultivated land and developing agricultural land to reduce desertification land. Benefiting from the construction of the Three-North Protection Project, the ecological environment of the L Plateau and the transition belt between the NAS Region and L Plateau have been significantly improved, and the cropping intensity largely increased from 2001 to 2016.
In addition, agricultural yields and agricultural machinery power inputs have greater impacts on cropping intensity changes in the regions with developed economies than in the regions with undeveloped economies. Economic development may facilitate investment in agricultural production technology research and development, which can improve arable land yields and efficiency. The developments of advanced agricultural technology and agricultural machinery inputs are the intrinsic reasons to maintain cropping intensity in regions with rapid economic development. It is very obvious that the agricultural structure in western China has undergone tremendous changes during 2001–2016. Further, agricultural support policies and agricultural technology development have achieved success in promoting agricultural industry. The driving factors that always have higher PD values in the southern regions are GYlabor, AMPALA, PMHAALA and FIcapita.
6. Conclusions
Monitoring the spatial distribution of multiple cropping and its changes, as well as understanding the underlying driving factors, could support the development of the land use simulation model and decision-making for efficient farmland planning. Due to limitations in the available intensity dataset and comprehensive interactions between complex factors, few studies have focused on the cropping intensity changes under rapid socioeconomic development and the uncertainty of their driving forces. In this study, we analyzed the spatial patterns and changing trends of cropping intensity from 2001 to 2016 and used the Geo-Detector to analyze the socioeconomic driving forces that affect cropping intensity and its changes at the county level in nine major agricultural divisions in China. The results suggest that the annual cropping intensity in nine agricultural zones all increased from 2001 to 2016 under rapid socioeconomic development, and seven of them increased significantly. The cropping cycles that reflect cropping practices in the main regions did not change much, but fallow practices in western China decreased largely from 2001–2016. For example, there was a 23% reduction of fallow in the NAS regions, and a 51% reduction of fallow in the L Plateau.
When using the Geo-Detector to analyze the influential factors that determine the spatial distribution of cropping intensity, it was found that the output efficiency of cultivated lands has considerable impacts on cropping intensity in the plains, while topography is the most influential factor in mountainous areas. In addition, the results suggest that the factors that affect cropping intensity changes are mainly related to the output efficiency of cultivated lands in the northern region, while they are mainly related to regional economic developments in the southern region, but the factors that influence cropping intensity changes are not the same for each individual agricultural region. These results could provide a reference to the driving forces and underlying mechanism behind cropping intensity and its changes in China and also provide profound insight into sustainable policymaking for arable lands. In the plains of northern China, limited by climatic conditions, the government can implement a series of relevant policies to promote cropping intensity, such as improving arable land quality, increasing financial investment in agricultural technology research and development, and increasing financial investment in agricultural machinery. This arid and semiarid region still has room to increase cropping intensity, and the policies for improving the ecological environment by encouraging sustainable agriculture development are worth developing. Further, in the SC, agricultural-promoting policies and industrial production-promoting policies are mutually exclusive, so it will be necessary to develop some policies to achieve a win–win situation for both agricultural and industrial development. Hence, with scientific guidance from influencing factors and driving mechanisms in each agricultural zone, the development of agricultural policies will become more focused.
Author Contributions
L.L. contributed to the literature review, research design, result interpretation, manuscript writing, and revision; Z.A. contributed to the Geo-Dectector in R implementation and data analysis; Y.Z. contributed to research design and manuscript revision; X.L. contributed to data collection and manuscript revision. All authors commented on and approved the final manuscript.
Funding
This research is supported by the Key Lab of Guangdong for Utilization of Remote Sensing and Geographical Information System, Guangzhou Institute of Geography (grant no. 2017B030314138), National Natural Science Foundation Project of China (grant no. 41901345), the Project of the Natural Science Foundation of Guangdong Province, China (grant No. 2018A030313683), and the Project of Science and Technology Planning Project of Guangdong Province, China (grant No. 2019B020208002).
Acknowledgments
We thank anonymous reviewers for their constructive comments.
Conflicts of Interest
The authors declare no conflict of interest.
Figures and Tables
Enlarge this image.Figure 1. The map of nine agricultural zones in China with a backdrop of the Chinese land cover map. Both maps were obtained from the Resource and Environment Data Cloud Platform (http://www.resdc.cn/) supported by the Institute of Geographic Sciences and Natural Resource Research in the Chinese Academy of Sciences.
Enlarge this image.Figure 2. Annual mean cropping intensity maps for each individual year and the sixteen-year mean cropping intensity map in China from 2001 to 2016.
Enlarge this image.Figure 3. Changes in the mean cropping intensity for nine individual agricultural zones in China in 2001–2016.
Enlarge this image.Figure 4. Changes in cropping intensity in China from 2001 to 2016 and at every four year interval. Changes in cropping intensity in China (a) from 2001 to 2016; (b) from 2001 to 2004; (c) from 2005 to 2008; (d) from 2009 to 2012; (e) from 2013 to 2016.
Enlarge this image.Figure 5. The spatial distribution of cropping intensity trends for each pixel within 5 km grids are shown for (a) changing rates in linear trends from 2001 to 2016 and (b) significances in linear trends from 2001 to 2016 (significantly increasing, SigI; increasing but not significant, Non-SigI; decreasing but not significant, Non-SigD; significantly decreasing, SigD).
The equations for the potential driving factors.
| Potential Driving Factors | Equations |
|---|---|
| GY | GY = total grain yield output / total arable land area |
| GY labor | GYlabor = total grain yield output / household registration population |
| ALA captita | ALAcaptita = total arable land area / household registration population |
| AMP ALA | AMPALA = agricultural machinery power / total arable land area |
| PMHA ALA | PMHAALA = total machinery harvest area / total arable land area |
| LF ALA | LFALA = household registration population / total arable land area |
| GR AG | GRAG = increased product in agricultural sector / regional GDP |
| GR PIG | GRPIG = increased product in primary industry / regional GDP |
| FI capita | FIcapita = household balance of farmers’ savings and deposits / household registration population |
Brief summary for the statistics associated with the annual mean cropping intensity in 2001–2016 and its changing trend. Trend slope and intercept denote the linear trend of the annual mean cropping intensity changes in 2001–2016. Z-value, p-value and significance denote the Mann–Kendall tests applied to the trend of annual mean cropping intensity changes in 2001–2016.
| NE Plain | NAS Region | HHH Plain | L Plateau | QT Plain | MLY Plain | SBS Region | SC | YG Plain | |
|---|---|---|---|---|---|---|---|---|---|
| annual mean | 97.8% | 61.3% | 135.9% | 85.7% | 41.8% | 138.0% | 97.3% | 171.4% | 128.0% |
| trend slope | 0.001 | 0.005 | 0.006 | 0.013 | 0.001 | 0.004 | 0.009 | 0.017 | 0.012 |
| trend intercept | −0.696 | −9.725 | −11.192 | −26.019 | −0.653 | −6.795 | −17.119 | −32.596 | −23.651 |
| Z-value | 1.711 | 2.521 | 3.332 | 3.962 | 0.720 | 1.981 | 3.332 | 3.962 | 3.242 |
| p-value | 0.087 | 0.012 | 0.001 | 0.000 | 0.471 | 0.048 | 0.001 | 0.000 | 0.001 |
| significance | FALSE | TRUE | TRUE | TRUE | FALSE | TRUE | TRUE | TRUE | TRUE |
Proportions of cropping practices for the arable lands in different agricultural zones in both 2001 and 2016.
| Agricultural Zones | Single | Double | Triple | Fallow | ||||
|---|---|---|---|---|---|---|---|---|
| 2001 | 2016 | 2001 | 2016 | 2001 | 2016 | 2001 | 2016 | |
| NE Plain | 0.89 | 0.98 | 0.00 | 0.00 | 0.00 | 0.00 | 0.11 | 0.02 |
| NAS Region | 0.57 | 0.77 | 0.01 | 0.04 | 0.00 | 0.00 | 0.42 | 0.19 |
| HHH Plain | 0.59 | 0.47 | 0.36 | 0.51 | 0.00 | 0.00 | 0.05 | 0.02 |
| L Plateau | 0.39 | 0.74 | 0.09 | 0.16 | 0.00 | 0.00 | 0.52 | 0.1 |
| QT Plateau | 0.49 | 0.49 | 0.02 | 0.02 | 0.00 | 0.00 | 0.49 | 0.49 |
| MLY Plain | 0.63 | 0.52 | 0.35 | 0.46 | 0.01 | 0.01 | 0.01 | 0.01 |
| SBS Regions | 0.88 | 0.75 | 0.05 | 0.2 | 0.00 | 0.01 | 0.07 | 0.03 |
| SC | 0.4 | 0.24 | 0.56 | 0.57 | 0.03 | 0.18 | 0.01 | 0.01 |
| YG Plateau | 0.78 | 0.63 | 0.21 | 0.3 | 0.00 | 0.06 | 0.01 | 0.01 |
The power of determinant (PD) values as derived from the Geo-Detector for the influence of 13 factors on cropping intensity. The maximum PD values for each region are underlined.
| PD | NE Plain | NAS Region | HHH Plain | L Plateau | QT Plain | MLY Plain | SBS Region | SC | YG Plain |
|---|---|---|---|---|---|---|---|---|---|
| GY | 0.038 | 0.193 | 0.551 | 0.246 | 0.013 | 0.146 | 0.260 | 0.002 | 0.110 |
| GY labor | 0.093 | 0.108 | 0.276 | 0.012 | 0.008 | 0.220 | 0.391 | 0.101 | 0.074 |
| ALA captita | 0.332 | 0.033 | 0.333 | 0.319 | 0.080 | 0.037 | 0.051 | 0.071 | 0.083 |
| AMP ALA | 0.009 | 0.022 | 0.301 | 0.171 | 0.015 | 0.024 | 0.099 | 0.018 | 0.049 |
| PMHA ALA | 0.044 | 0.009 | 0.547 | 0.096 | 0.007 | 0.180 | 0.047 | 0.025 | 0.232 |
| LF ALA | 0.063 | 0.053 | 0.233 | 0.318 | 0.087 | 0.016 | 0.019 | 0.054 | 0.077 |
| TE | 0.066 | 0.156 | 0.255 | 0.118 | 0.094 | 0.153 | 0.468 | 0.252 | 0.307 |
| TS | 0.123 | 0.070 | 0.208 | 0.133 | 0.149 | 0.222 | 0.096 | 0.219 | 0.102 |
| GPIE | 0.024 | 0.045 | 0.016 | 0.060 | 0.033 | 0.097 | 0.065 | 0.008 | 0.049 |
| GDP | 0.016 | 0.052 | 0.034 | 0.012 | 0.062 | 0.055 | 0.088 | 0.040 | 0.039 |
| GR AG | 0.054 | 0.125 | 0.016 | 0.104 | 0.031 | 0.004 | 0.066 | 0.003 | 0.017 |
| GR PIG | 0.021 | 0.133 | 0.017 | 0.134 | 0.058 | 0.010 | 0.080 | 0.035 | 0.018 |
| FI capita | 0.130 | 0.026 | 0.085 | 0.014 | 0.064 | 0.016 | 0.028 | 0.000 | 0.000 |
The PD values as derived from the Geo-Detector for the influence of thirteen factors on cropping intensity changes. The maximum PD values for each region are underlined.
| PD | NE Plain | NAS Region | HHH Plain | L Plateau | QT Plain | MLY Plain | SBS Region | SC | YG Plain |
|---|---|---|---|---|---|---|---|---|---|
| GY | 0.004 | 0.126 | 0.409 | 0.160 | 0.014 | 0.036 | 0.282 | 0.017 | 0.147 |
| GY labor | 0.074 | 0.110 | 0.204 | 0.027 | 0.138 | 0.064 | 0.284 | 0.643 | 0.210 |
| ALA captita | 0.304 | 0.001 | 0.284 | 0.196 | 0.041 | 0.211 | 0.264 | 0.004 | 0.246 |
| AMP ALA | 0.006 | 0.003 | 0.293 | 0.046 | 0.014 | 0.067 | 0.257 | 0.621 | 0.089 |
| PMHA ALA | 0.039 | 0.013 | 0.395 | 0.025 | 0.590 | 0.038 | 0.001 | 0.008 | 0.120 |
| LF ALA | 0.029 | 0.000 | 0.245 | 0.188 | 0.014 | 0.015 | 0.012 | 0.552 | 0.200 |
| TE | 0.065 | 0.171 | 0.211 | 0.056 | 0.071 | 0.062 | 0.158 | 0.095 | 0.311 |
| TS | 0.099 | 0.004 | 0.162 | 0.166 | 0.251 | 0.006 | 0.051 | 0.062 | 0.075 |
| GPIE | 0.017 | 0.051 | 0.006 | 0.000 | 0.139 | 0.474 | 0.230 | 0.788 | 0.255 |
| GDP | 0.011 | 0.022 | 0.014 | 0.151 | 0.275 | 0.038 | 0.152 | 0.445 | 0.081 |
| GR AG | 0.042 | 0.097 | 0.012 | 0.041 | 0.014 | 0.133 | 0.099 | 0.414 | 0.139 |
| GR PIG | 0.019 | 0.075 | 0.009 | 0.057 | 0.056 | 0.200 | 0.139 | 0.052 | 0.004 |
| FI capita | 0.071 | 0.015 | 0.068 | 0.011 | 0.427 | 0.299 | 0.709 | 0.608 | 0.720 |
The PD values as derived from the Geo-Detector based on a standard deviation classification method for the influence of 13 factors on cropping intensity changes. The maximum PD values for each region are underlined.
| PD | NE Plain | NAS Region | HHH Plain | L Plateau | QT Plain | MLY Plain | SBS Region | SC | YG Plain |
|---|---|---|---|---|---|---|---|---|---|
| GY | 0.007 | 0.134 | 0.463 | 0.203 | 0.131 | 0.057 | 0.167 | 0.070 | 0.222 |
| GY labor | 0.091 | 0.138 | 0.202 | 0.023 | 0.140 | 0.116 | 0.082 | 0.275 | 0.320 |
| ALA captita | 0.148 | 0.046 | 0.324 | 0.226 | 0.133 | 0.062 | 0.125 | 0.206 | 0.327 |
| AMP ALA | 0.024 | 0.042 | 0.266 | 0.108 | 0.649 | 0.120 | 0.218 | 0.202 | 0.148 |
| PMHA ALA | 0.030 | 0.023 | 0.452 | 0.037 | 0.242 | 0.143 | 0.148 | 0.050 | 0.039 |
| LF ALA | 0.155 | 0.006 | 0.277 | 0.254 | 0.649 | 0.026 | 0.171 | 0.360 | 0.068 |
| TE | 0.063 | 0.206 | 0.222 | 0.051 | 0.059 | 0.025 | 0.526 | 0.028 | 0.505 |
| TS | 0.119 | 0.008 | 0.221 | 0.042 | 0.062 | 0.064 | 0.014 | 0.241 | 0.002 |
| GPIE | 0.015 | 0.006 | 0.082 | 0.056 | 0.010 | 0.081 | 0.085 | 0.373 | 0.164 |
| GDP | 0.015 | 0.021 | 0.108 | 0.025 | 0.106 | 0.088 | 0.084 | 0.195 | 0.068 |
| GR AG | 0.079 | 0.089 | 0.022 | 0.069 | 0.160 | 0.200 | 0.221 | 0.095 | 0.065 |
| GR PIG | 0.005 | 0.091 | 0.006 | 0.085 | 0.260 | 0.074 | 0.227 | 0.233 | 0.024 |
| FI capita | 0.100 | 0.007 | 0.126 | 0.007 | 0.213 | 0.135 | 0.569 | 0.124 | 0.296 |
© 2019 by the authors.
Details
; Zhao, Yaolong 3 ; Liu, Xulong 4 1 School of Management, Guangdong University of Technology, Guangzhou 510520, China;
2 Guangdong Provincial Key Laboratory of Urbanization and Geo-Simulation, Sun Yat-sen University, Guangzhou 510275, China
3 School of Geography, Guangdong Research Center of Smart Homeland Engineering, South China Normal University, Guangzhou 510631, China;
4 Key Lab of Guangdong for Utilization of Remote Sensing and Geographical Information System, Guangzhou Institute of Geography, Guangzhou 510070, China;