**1. Introduction**

Agricultural production is the cornerstone for the survival and development of human beings [1]. Due to the impact of urbanization, development of the market economy, the international food trade and climate change, dramatic changes have taken place in agricultural land systems. The most famous examples are the conversion of croplands into built-up areas in the Yangtze River Delta and Pearl River Delta of China [2,3], the phenomenon of double-cropping transitioning to single-cropping in the middle reaches of the Yangtze River Valley [4], the non-grain use of croplands in China [5], the collapse of soybean planting in Northeast China [6,7], and the abandoned croplands in mountainous area in South China [8]. Furthermore, locust disasters [9], armed conflict [10], COVID-19 [11,12], etc., aggravated the uncertainty of food production. The world has encountered the most serious food

**Citation:** Tao, J.; Zhang, X.; Liu, Y.; Jiang, Q.; Zhou, Y. Estimating Agricultural Cropping Intensity Using a New Temporal Mixture Analysis Method from Time Series MODIS. *Remote Sens.* **2023**, *15*, 4712. https://doi.org/10.3390/rs15194712

Academic Editor: Georgios Mallinis

Received: 9 June 2023 Revised: 10 September 2023 Accepted: 11 September 2023 Published: 26 September 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

crisis in the past 50 years [13]. On the other hand, continuous agricultural intensification in some regions has caused serious ecological and environmental problems, such as excessive consumption of water resources [14], land degradation [15], agricultural non-point source pollution [16] and greenhouse gas emissions [17]. Therefore, agricultural cropping intensity is an important input for evaluating the sustainable development of agriculture.

The existing methods of cropping intensity mapping include cropping frequency (CF) and the multiple cropping index (MCI). The calculation of MCI is mostly based on statistical data, ignoring the spatial heterogeneity within the administrative regions [18]. CF is a nominal measurement that divides crop patterns into single-cropping and doublecropping, lacking quantitative measurement of cropping intensity. Most existing methods for mapping CF are based on counting the number of peaks in coarse resolution remote sensing images such as the Moderate Resolution Imaging Spectroradiometer (MODIS) vegetation index profiles [19–23]. When a coarse resolution image is used to monitor cropping frequency, the mixed pixel problem is particularly serious due to the fragmented landscapes, small patch size of croplands and complicated crop patterns in Central and South China.

There has been some progress recently on the construction of new indicators or new methods for cropping intensity mapping. Time series high-resolution images such as those from Sentinel-2 are used to monitor cropping intensity, alleviating the problem of mixed pixels to some extent [24–26]. However, these methods are highly dependent on data availability and constrained by complicated data pre-processing. The authors used Bayesian network to obtain cropping intensity with interval measurement using time series MODIS data [27]. This method is highly dependent on training samples, and the accuracy of the model output depends on the quality of the samples. In any case, the above existing methods of cropping intensity mapping cannot meet the requirements of precision agriculture, and new methods for fine cropping intensity mapping suitable for Central and South China are still lacking.

Although continuous accumulated high-resolution remote sensing images are helpful for precise cropping intensity mapping, this convenience has only existed since the launch of Landsat 8 in 2013 and the launch of Sentinel-2 in 2015. Medium-resolution data over long time series can bridge the gap and push back observations to the year 2000 or even earlier. Time series remote sensing data contain seasonal variation of vegetation and abundant vegetation phenology information [28]. Intra-annual dense time series images carry information on multiple cropping of crops. However, this information is usually used for crop mapping [29,30], phenology monitoring [31], etc. Vegetation phenology information contained in time series remote sensing data is not yet fully explored for cropping intensity mapping. Dense time-series remote sensing images have enough repeated observation (MODIS MOD13Q1 has 23 observations a year) and make this work possible [32].

Mixed pixels are more common in coarse resolution images and require mixture analysis technology to decompose pixels to the abundance level. There are two mixture analysis techniques: spectral mixture analysis (SMA) and temporal mixture analysis (TMA). SMA methods are generally used to decompose spectral remote sensing data, among which linear spectral mixture analysis (LSMA) is more often used. Based on the principle of SMA, TMA is a method to decompose mixed pixels by replacing spectral data with multi-temporal data.

The selection of endmembers (including their number and types) is another important issue for unmixing. The two-endmember model is suitable for the decomposition of natural land-cover types, while the three-endmember model is suitable for the decomposition of land-cover types under human disturbance [33]. For built-up areas, the most commonly used endmember selection method is the vegetation–impervious surface–soil endmember model proposed by Ridd et al. [34]. For non-built-up areas, vegetation–soil–shadow (or dry vegetation) endmember model is generally used [35]. The endmember selection is usually scene-dependent when the TMA method is used: for example, the vegetation–impervious surface–soil model [36], forests–multiple cropping–single cropping–non-vegetation model [32], and grass–corn–winter wheat–alfalfa model [37]. With proper endmember design, large scale, long temporal spans and fine cropping intensity are possible.

To aim for large-scale, long-temporal-span and fine-resolution agricultural cropping intensity estimation, this study presents a new TMA method for estimating cropping intensity from historical archived time series coarse resolution remote sensing data. The method includes: (1) constructing the feature space, finding the unique endmembers to estimate the abundance level cropping intensity; (2) exploring the feasibility of manually delineating endmembers and the effectiveness of the method on different regions with varied completeness of endmember land-cover types. This work will be of great significance for fine cropping intensity mapping at large scale and long-time series.
