Monitoring Cropping Intensity Dynamics across the North China Plain from 1982 to 2018 Using GLASS LAI Products

: China is a large grain producer and consumer. Thus, obtaining information about the cropping intensity (CI) in cultivated land, as well as understanding the intensiﬁed utilization of cultivated land, is important to ensuring an increased grain production and food security for China. This study aims to detect and map the changes in CI over a period of 36 years across China’s core grain-producing area—the North China Plain (NCP)— using remotely sensed leaf area index (LAI) time series data acquired by the Global LAnd Surface Satellite (GLASS) products. We ﬁrst selected 2132 sample points that consisted entirely, or almost entirely, of cultivated cropland from all pixels; the biennial LAI curves for the sample points were then extracted; the Savitzky–Golay ﬁlter and second-order difference algorithm were then applied to reconstruct the biennial LAI curves and obtain the number of peaks in these curves. In addition, the multiple cropping index (MCI) was calculated to represent the CI. Finally, the spatial distribution of the CI of cultivated land on the NCP was mapped from 1982 to 2018 using a geo-statistical kriging approach. Spatially, the results indicate that the CI of cultivated land over the NCP exhibits a distinct spatial pattern that conforms to “high in the south, low in the north”. The single cropping system (SCS) mainly occurred in the higher latitude area ranging from 37.04 ◦ N to 42.54 ◦ N, and the double cropping system (DCS) mainly existed in the lower latitude area between 31.95 ◦ N and 39.97 ◦ N. Temporally, the CI increased over the study period, but there were some large ﬂuctuations in CI from 1982 to 1998 and it maintained relatively stable since 2000. Across the NCP, 68.14% of cultivated land experienced a signiﬁcant increase in CI during the 36-year period, while only 3.87% showed a signiﬁcant decrease. We also found that, between 1982 and 2018, the northern boundary of the area for DCS underwent a signiﬁcant westward expansion and northward movement. Our results show a good degree of consistency with statistical data and previous research and also help to improve the reliability of satellite-based identiﬁcation of CI using low spatial resolution LAI products. The results provide important information that can be used for analyzing and evaluating the rational utilization of cultivated land resources; thus, ensuring food security and realizing agricultural sustainability not only for the NCP, but for China as a whole. These results also highlight the value of satellite remote sensing to the long-term monitoring of cropping intensity at large scales.


Introduction
As a result of climate change, there has been an obvious increase in temperature in recent years [1], and this has had a great impact on the utilization of regional agricul-

Study Area
The NCP mainly consists of the provinces of Hebei, Shandong and Henan as well as the municipalities of Beijing and Tianjin; it also includes parts of northern Jiangsu and Anhui provinces [39]. In this study, our study area was located mainly between 32 • N and 42 • N and 113 • E and 120 • E in Lambert conformal conic projection; thus, excluding the parts of the plain that are in Jiangsu and Anhui provinces ( Figure 1). With flat terrain, deep and fertile soil and a great number of rivers and lakes, the NCP is a typical alluvial plain that has an altitude of less than 50 meters over most of its area [40]. The annual average temperature is 12-15 • C, with an annual accumulated temperature (≥0 • C) of 4100-5400 • C and an annual frost-free period of 190-220 days. The average annual precipitation is 500-1000 mm. These climate conditions meet the requirements of a one-year double cropping system (DCS) [41]. The study area is an important crop production region for China, and its cultivated area accounts for 18% of the total area of cultivated land in China; the area's grain output accounts for 21% of national agricultural production [40]. Dryland farming is the main farming system used on the NCP; the main crops are wheat and maize. Excepting the single cropping system (SCS) over one year in parts of northern region in Beijing and Hebei province and the triple cropping system (TCS) over two years in part of southern Henan province, the DCS with winter wheat and summer maize planting is the main farming system in NCP. This forms a comprehensive agricultural system with well-developed irrigation and features of grain and cotton cultivation.

Remotely Sensed LAI Data
The remotely sensed LAI data used in this study were GLASS LAI data for the period 1982 to 2018, provided by the GLCF (Global Land Cover Facility) (http://glassproduct.bnu.edu.cn/, accessed 3 on May 2020) ( Table 1). The GLASS LAI dataset has a temporal resolution of 8 days and is a 37-year (1981-2018) Leaf Area Index remote sensing product [42][43][44]. The LAI product covering the period 1981-1999 has a spatial resolution of 0.05 • × 0.05 • and uses a geographic projection; the LAI data for 2000 to 2018 uses an Integrated Sinusoidal projection with three spatial resolutions of 0.05 • × 0.05 • , 1 × 1 km and 500 × 500 m [45][46][47]. In this study, the GLASS LAI data with a spatial resolution of 0.05 • × 0.05 • were used to maintain continuity.  1980, 1990, 1995, 2000, 2005, 2010, 2015 and 2018. This database has a scale of 1:100,000 and was provided by the Science Data Center of Resources and Environment, Chinese Academy of Sciences (http://www.resdc.cn, accessed on 3 May 2020). We used these data to extract location and area information related to cultivated land on the NCP. Based on these information, we extracted the data of unchanged cultivated land in the land use dataset to obtain spatio-temporal information about the CI over the NCP; in this way, it could avoid the errors for mapping CI caused by the distribution and area differences of cultivated land in different years.

Statistical Data
Statistical agricultural data in province scale, including the cultivated land area and total sown area in the study region, were obtained from the China Statistical Yearbook (http://www.stats.gov.cn/tjsj/ndsj/, accessed on 11 August 2020) and the Beijing, Tianjin, Hebei, Shandong and Henan Statistical Yearbooks for 1982 to 2018. These data were used to calculate the reference MCI for each subregion as well as the whole study area, and to comparatively analyze and validate the CI extracting results using our proposed method.

Algorithms for Identifying the Cropping Intensity
For crops, long-time series of the LAI can describe the characteristics of annual or interannual growth changes and reflect changes in the life cycles of crops [47,48]. From sowing, emergence, jointing and heading to harvesting, the LAI values of crops also undergo a dynamic process of increase to peak to decrease. Therefore, by constructing LAI time series for crop growth and calculating the MCI from the number of peaks in these curves, the CI of cultivated land can be monitored [47]. The algorithm based on the long-term GLASS LAI series that was used for determining the CI of cultivated land is shown below (Figure 2). Based on the information about unchanged cultivated land on the NCP that was extracted from the land use dataset in 1980,1990,1995,2000,2005,2010, 2015 and 2018, 3500 cultivated land sample points were randomly selected using a gridded spatial sampling strategy (by on grid per 5 km × 5 km; 18,308 grids in total); these, then, formed the preliminary sample used in this study. Next, each sample point was validated and verified using high-resolution imagery (e.g., Google Earth version 7.1.2.2019 imagery) and pixels that consisted of more than 80% cultivated land were then selected, meaning that all of these pixels were either pure or nearly pure cultivated land pixels. In this way, 2132 cultivated land sample points covering the NCP were obtained ( Figure 1B), of which 110 were the validation sample.
The land use on the NCP was highly fragmented and, to a certain extent, this sampling method could not only reduce the errors caused by mixed pixels in low-and mediumresolution data (i.e., GLASS LAI data) when trying to extract the CI, but also improve the accuracy of kriging spatial interpolation by increasing the density of sampling points appropriately (by one point per 250 km 2 on average). In this study, we used the Savitzky-Golay (S-G) filtering method to remove the noisy components in the long LAI time series and to reconstruct the LAI curves for crop growth in the areas of cultivated land [49]. The formula used was: where Y * j (j = 1, . . . , N) is the filter LAI value, Y j (j = 1, . . . , N) is the original LAI value, C i is the coefficient for the ith LAI value of the filter, and N is the number of convoluting integers, which is equal to the size of the smoothing window (2m + 1). The smoothing array consists of 2m + 1 points, where m is the half-width of the smoothing window [50,51]. The C i can be obtained from Steinier et al. [52].
(2) Extracting the number of peaks in the LAI curves We employed a second-order difference algorithm to identify the number of peaks in the reconstructed LAI curves [53,54] and, then, extracted the MCI for cultivated land within the NCP from 1982 to 2018 [55]: where N LAI i is the LAI value of a sampling pixel with i being the time series phase for the GLASS LAI imagery. s (1) , s (2) and s (3) are the detected parameters in Equations (2)-(4), respectively; the meanings of these are explained in Section 3.2.

CI Estimation
(1) Calculating the MCIs of the sample points Based on the relationship between the LAI time series for cultivated land and the MCIs, the MCIs could easily be converted to the number of peaks in the LAI curves [56]. Given the number of peaks for a given sample point, the number of peaks in the LAI curves of the crops (i.e., the MCIs of the cultivated land sample points) could be obtained from the ratio of the number of peaks to the total number of sample points: where F i represents the frequency of wave peaks in the ith year, ∑ s j (i) is the number of peaks in the LAI curves for all of the sample points during the ith year, ∑ p j (i) is the total number of sample points during the ith year and j (equal to 1, 2, 3, . . . ) is the sample point number.
(2) Spatial mapping of CI In this study, the numeric value of MCI was employed to measure the CI. Based on the MCIs of the cultivated land sample points, time series maps showing the spatial distribution of the CI of the cultivated land on the NCP were produced using the geostatistical Kriging interpolation model. According to the characteristics of the sample data, the Ordinary Kriging model was used for the spatial interpolation [57]. The relevant equation was: where Z(s i ) is the measured value at the ith position, λ i is the unknown weight of the measured value at the ith position, s 0 is the predicted position and N is the number of measured values.
Following this, we performed an accuracy analysis by comparing different semivariogram models (i.e., Gaussian, spherical, exponential and linear) to select the most optimal Kriging mode for spatial interpolation [58]. Additionally, the spatial patterns of CI for the cultivated land on the NCP from 1982 to 2018 were mapped using the optimal Kriging mode; i.e., the linear Kriging mode.

Smoothed LAI Time Series Curves
The S-G filtering model was used to smooth and reconstruct the LAI time series, and the results were satisfactory. The reconstructed LAI curves for the sample points on cultivated land reflected the processes of crop growth and development well, and the number of peaks in the curves was clear; thus, allowing the multiple cropping on the NCP to be accurately characterized ( Figure 3). Figure 3A-D shows, respectively, the LAI curves for bare or fallow land, SCS, DCS and triple cropping system. These curves contained 0, 1, 2 and 3 peaks, respectively; curves such as these were used as the basis for extracting the CI of cultivated land from remote sensing data.

Extraction Results of the Peaks of the LAI Curve and the Determination of Valid Correct Peaks
Taking the extraction of the peaks for a sample point from 1982, as an example, we obtained the results for each step of the calculation given by Equations (2)-(4) according to the principles of the second-order difference algorithm. Figure 4A shows the positions and numbers of peaks and troughs in the smoothed LAI time series, whereas Figure 4B-D shows the positions of the peaks and troughs in the curve determined using Equations (2)-(4), respectively. In this way, the complete process used to extract the number of peaks using the second-order difference algorithm can be understood. It can be seen from Figure 4 that peak one and peak two were true peaks, while peak three was a false peak. In order to decide which were the true peaks, the following conditions were applied [59][60][61].

•
The peak shown occurred between day 120 and day 300 of a year; • In the case of a unimodal peak, the LAI value at the peak should not be less than two; • In the case of multimodal peaks, the LAI value at the smaller peak should not be less than 40% of the value of the highest peak.

Comparison of the MCIs in This Study with Statistical Data
As the index to measure CI in this study, the extraction precision of MCI was closely related to the accuracy of the spatio-temporal change analysis for CI. The extraction accuracy of the MCI was determined by comparing the results with the statistical data for the MCI from 1982 to 2018 (the data for some years were missing) ( Figure 5). Our results based on remote sensing LAI data were generally highly consistent with those derived from the statistical data. Although in 1982, 1996 and 2018, the accuracy of the CI for Beijing was less than 70%; for about 85% of the remote sensing extraction results the accuracy was more than 80% (and over 90% for 61% of the results). This showed that the results of this study were fairly consistent with the statistical data, further verifying that the CI results retrieved by remote sensing were reliable. The accuracies of MCIs are calculated by formulas (7) and (8) as follow: where MCI (Extracted) is the MCI in this study; MCI (Statistical) is the MCI calculated based on statistical data; T (sown) is the total sown area of crops; T (real) is the cultivated land area.

Comparison of the MCI in This Study with Other Research Results
Based on the MCI, comparisons were also made between the results of the current study and the monitoring results obtained by Fan and Wu [62] in 2000 and 2002 and Li et al. [26] in 2002, 2008 and 2014 ( Figure 6).
It can be seen from Figure 6 that the results of this study were in good agreement with the earlier researches, especially with regard to the results obtained by Li et al. [26], for which there was an extremely small difference of less than 10% in 2014. The differences of MCI in Tianjin and Beijing were the largest, followed by Shandong Province, and the difference in Henan Province was the smallest. However, for different provinces or years, there were some slight differences between the results of our study and the others. In 2002, Henan had the greatest discrepancy between our results and the results obtained by Li et al. [26] at 28.77%. The difference was less than 5% for Hebei and Tianjin in 2002, Hebei in 2008, as well as Hebei, Henan, Beijing and Tianjin in 2014. This suggests that, overall, there was a high degree of consistency between the results of our study and those of earlier studies; however, due to various factors such as the monitoring methods used and the datasets on which the studies were based, some discrepancies remained. A comparison between the results of our study, the results of previous research and the statistical data showed that our results were closer to the statistical data than the results obtained in previous research. In this study, the rounded numeric value of MCI was employed to measure the cropping system (CS). The value of MCI below 150% and exceeding 50% meant it was SCS in that pixel, while the value 150-250% meant it was DCS and exceeded 50% mean TCS. We also used cross validation to verify the precision of the CS of cultivated land. Firstly, 110 samples were randomly selected from the 2132 cultivated land sample points. Based on the LAI time series data, 2090 group data of cropping systems for every two years from 1982 to 2018 were directly extracted from the 110 validated sample points. In addition, spatial data for the CS used on the NCP were also extracted from the sample points after interpolation and, then, the accuracy of the spatial mapping results was validated by analyzing the two data error matrices ( Table 2). The results showed that the overall accuracy of the spatial mapping in this study was 81.87%; comparing with accuracy for fallow fields and SCS, the precision of the user for the DCS was generally high, while the producer accuracy was relatively lower. The fallow field and SCS were rarely found on the NCP and their user accuracies were relatively low. There was also a high degree of consistency between the extraction results for the SCS and DCS-the corresponding kappa coefficient was 0.65, which indicated that the extraction results for cropping systems on the NCP obtained using interpolation corresponded well to the results obtained directly from LAI data. According to the results shown in Figure 7, in most parts of the NCP, the CI of cultivated land was in the range of 80% to 200%, but the intensity varied widely across the region. In the central and southern parts of the NCP, the CI was relatively high, while the CI was relatively low in the high-latitude and marginal mountainous areas in the north of the NCP. As far as the average CI of cultivated land on the NCP from 1982 to 2018 was  Based on the results for the MCI in the NCP from 1982 to 2018 that were already found, we carried out a further analysis of the spatial distribution patterns of the Cropping System (CS) in the study area. The CS on the NCP was characterized by the SCS and DCS (Figure 9), with the DCS being used over about 76% of the total cultivated area of the plain. The cropping systems in the central and southern parts of the NCP were dominated by the DCS. In the higher latitude regions in the north of the plain, due to less abundant sunshine and water and lower temperatures, the SCS dominated. In the eastern and western parts of the NCP, due to the presence of hills and mountains, the cultivation pattern alternated between the SCS and DCS. The proportions of land where the different cropping systems were used in different provinces and cities were also quite different. From 1982 to 2018, the main cropping pattern used in Shandong and in Henan, in particular, was the DCS. In Henan, the proportion of land where the DCS was used was as high as 97.16%; in Shandong, it was close to 81.12%. The SCS was the most widely used CS in Tianjin and Beijing, accounting for about 65.1% and 58.82% of the total cultivated land area, respectively. In Hebei Province, the DCS and SCS were relatively balanced.

Temporal Patterns in CI
Based on the long-term average of the CI of cultivated land on the NCP (for the period 1982-2018), an interannual trend line for the CI was drawn by using the simple linear regression method ( Figure 10A). Over the 36-year period, the CI showed an overall increase with some large fluctuations from 1982 to 1998 and relative stability since 2000. The lowest value of the CI was 108.66% in 1982 and the highest was 159.83% in 2016. The radar map ( Figure 10B) that was drawn for the five provinces and cities had a fairly wide left-hand side and a narrow right-hand side, which reflected the general increase in the CI. Among these five provinces and cities, Henan Province had the highest CI and, therefore, had the widest 'circle', which reached its maximum value in 1990. The CI values for Shandong and Hebei were also relatively high; Tianjin and Beijing had the lowest values, corresponding to the inner circles on the radar chart.
The spatial distribution pattern of the CI w then analyzed for the periods 1982-1992, 1992-2002, 2002-2012, 2012-2018 and 1982-2018 (Figure 11). It can be seen that the CI from 1982 to 2018 in more than 68.14% of the NCP exhibited an obvious increasing trend, where its value augmented more than 10%, and only 3.87% of the study area experienced a significant fall, where the CI decreased more than 10%. Accordingly, the CI for about 27.99% of the area was in the variation range −10-10% from 1982 to 2018, indicating that the CI in these areas was relatively stable over the 36-year period. The overall change from 2002 to 2012 and from 2012 to 2018 was small: during these periods, the areas with CI changed by −10-10% occupied 61.89% and 59.65% of the study area, respectively. Additionally, the trend was the same for the individual provinces and cities from 2002 to 2012 and 2012 to 2018. There was also an augmentation in the CI during the periods 1982-1992 and 1992-2002, with more than 60.81% and 52.12% of the NCP experiencing an increment.  We also analyzed the proportions of the total land area of the NCP where the different CS were used ( Figure 12) and the related trends to describe the temporal pattern of CI. Except for 1982 and 1984, the DCS accounted for a significantly larger proportion of the cultivated land than the SCS in all years within the period analyzed, with the biggest difference reaching 79.93% in 2008. The results also showed that the proportion of cultivated land on the NCP where the SCS was used was decreasing year by year. In 1982 and 1984, the proportions of land where the SCS was used were as high as 68.82% and 56.98%, respectively. After 1986, the proportion of land where the SCS was used dropped sharply, falling to only 10.03% in 2008. In contrast, the proportion of land where the DCS was used was rising and exceeded 70% in most years within the period analyzed. In 2008, it reached 89.96%, making it the prevailing cropping system on the NCP. This was a comprehensive result of the suitable temperature, sunshine, rainfall conditions, market factors, new technology availability and social structures, etc.

Applicability of GLASS LAI Data
Based on the conclusions of relevant references for the quality of GLASS LAI data [63,64], we set an example by extracting the AVHRR NDVI, MODIS NDVI and GLASS LAI time series of a random sample in 1982 and comparing them with each other. We also compared them for 2012. From the results shown in Figure 13, it can be seen that there were more pseudo peaks in the AVHRR and MODIS NDVI curves than in the GLASS LAI curves, and the outliers in the AVHRR NDVI curve were the most obvious. Meanwhile, the AVHRR and MODIS NDVI data also had a saturation problem under high-density vegetation-this could have increased the error in the corresponding CI and made the research results less reliable. Therefore, GLASS LAI data from 1982 to 2018 were used to obtain the CI of cultivated land covering a long period and over a wide range.

Latitudinal Zonality of CI in Cultivated Land on the NCP
Based on the frequency distribution histogram [65], the frequency of different cropping systems for cultivated land on the NCP at different latitudes was obtained in order to analyze the latitudinal zonality of the cropping systems ( Figure 14). The latitude of the study area ranged from 31.5 • to 42.6 • N, and the SCS was found across the whole area. However, it was mainly used at the higher latitudes ranging from 37.  Based on the known latitudinal limits for the DCS and the other cropping systems used on the NCP, the northern limit of the DCS-, i.e., the boundary between the SCS and DCS in the north of the NCP-was drawn for 1982, 1986, 1990, 1994, 1998, 2002, 2006, 2010, 2014 and 2018 by using the minimum bounding geometry tool in ArcGIS software ( Figure 15).
In this way, the distribution of the DCS by latitude from 1982 to 2018 and the dynamic changes in this distribution could be analyzed. As shown in Figure 15, overall, between 1982 and 2018, the northern boundary of the DCS moved northward. However, this movement was not steady. In 1982In , 1986In and 1990, the area of DCS expanded northward. In 1982, the northern limit of the DCS ran from Shijiazhuang City, Hebei Province, to Baoding City, the southern part of Langfang City and Tianjin City. By 1990, it had moved further north in Hebei Province, crossing Zhangjiakou City, Chengde City and the southern part of Qinhuangdao City. This was the furthest north that this limit reached during the 36-year period. Between 1990 and 1994, the DCS boundary continued to move northward within Qinhuangdao City; however, the rest of the line shifted south to lie across Beijing. Meanwhile, the eastern sections of this northern boundary of the DCS continued to shift southward from 1994 to 2010. In 2010 in particular, due to the strong influence of extreme drought events in the NCP, the northern boundary of the DCS was the closest to its southernmost limit after 1982. After 2010, the line, again, moved rapidly northward and ran from Tianjin to Chengde City in 2018.

Uncertainty Analysis
The results of this study were influenced by multiple sources of uncertainty, which mainly included the following.
(1) The surface topography of the NCP was relatively complex: there were obvious regional differences and a range of crop planting systems were used. In addition, intercropping and interplanting occurred in many regions. This affected crop growth curves, which could have led to confusion in the peak extraction results and, thus, affecting the results of the identification of CI using remote sensing. At the same time, the existence of a large number of mixed pixels in data with a low or medium spatial resolution remained a key problem that restricted the study of CI on the NCP. Therefore, the use of different remote sensing extraction methods corresponding to different regional conditions are important directions for further research. (2) Due to the complex surface conditions on the NCP, grid-based spatial random sampling would produce certain sampling errors, which was another source of uncertainty. Therefore, in subsequent research, we will use adaptive sampling in different regions to improve the sampling accuracy and efficiency.

Conclusions
Based on GLASS LAI data, in this study, S-G filtering and a second-order difference algorithm were used to extract the distribution of CI in cultivated land across the NCP from 1982 to 2018. We explored the temporal and spatial variation characteristics in the CI at different regional levels (for the whole NCP, and for provincial and municipal levels) based on the index of MCI in cultivated land.
By comparing the results of this study with statistical data and the results of previous studies, it was found that there was a high level of consistency between all of the results which showed that this study improved the accuracy of CI mapping of the NCP.
Spatially, based on the results, a spatial pattern of "high in the south and low in the north" was found for the CI and it was observed that there was a continuous increase in CI from 1982 to 2018 for the majority of the region, with the increase in Henan province being most noticeable. The SCS and DCS were the main cropping systems used on the NCP. The SCS was mainly used in the higher latitude area ranging from 37.04 • to 42.54 • N, whereas the DCS was mainly used in the lower latitude area between 31.95 • and 39.97 • N; Henan and Shandong provinces were typical DCS areas.
Temporally, over the 36-year period, the CI showed a rising overall with some large fluctuations from 1982 to 1998 and relative stability since 2000, and 68.14% of the cultivated land across the NCP experienced a significant increase in CI during the thirty-six-year period of the study. Only 3.87% of the study area experienced a significant decrease. Meanwhile, the proportion of cultivated land where the SCS was used was greater in 1982 and 1984, whereas the DCS accounted for a bigger proportion of the land from 1986 to 2018. In addition, between 1982 and 2018, the northern boundary of the DCS on the NCP underwent different degrees of westward expansion and northward movement, reaching Zhangjiakou City in Hebei Province and the northern part of Chengde City at its northernmost extent.
The results of this study provided important information for analyzing and evaluating the rational utilization of cultivated land resources, ensuring food security and realizing ecological protection on the NCP and in China. The results also highlighted the value of satellite remote sensing for the monitoring of crop planting structures at large scales and over a long term.

Data Availability Statement:
Publicly available Global LAnd Surface Satellite Leaf Area Index Products were analyzed in this study. This data can be found here: http://www.glass.umd.edu/ (accessed on 3 May 2020).