Biophysical Impact of Sunflower Crop Rotation on Agricultural Fields

Crop rotation is an important determining factor of crop productivity. Sustainable agriculture requires correct rules of crop rotation. Failure to comply with these rules can lead to deterioration of soil biochemical characteristics and land degradation. In Ukraine as well as in many other countries, sunflower monocropping is common practice and the impact of this fact should be studied to find the most precise rules to save the economic potential of land and minimize land degradation factors. This research provides an evaluation of the sunflower monocropping effect on the vegetation indices obtained from MODIS vegetation indices datasets for Ukraine as one of the countries with the biggest sunflower export in Europe. The crop rotation schemes are represented by their area proportions at the village level calculated based on the crop classification maps for 2016 to 2020. This representation gives the possibility to use regression models and f-test feature importance analysis to measure the impact of 3-year and 5-year crop rotation sequences. For these purposes, we use several models: a four-year binary representation model (model A1) and a model with all possible three-year crop rotation scheme representations (model B). These models gave the possibility to evaluate crop rotation schemes based on their biophysical impact on the next sunflower plantings and found that sunflower planting with an interval of three or more years is optimal in terms of the sustainability of


Introduction
Crop rotation is a necessary practice for sustainable food production. In many countries, crop rotation rules are regulated at the state level. This is especially important now, when Europe is moving toward the implementation of Green Deal practices to mitigate anthropogenic influences on climate change. The main task of agricultural management is to preserve the original composition and productivity of soils. The physical and chemical parameters of soils directly affect the condition of plants [1]. At the same time, noncompliance with recognized positive agricultural practices leads to the deterioration of soil fertility and soil quality due to the imbalance of nutrients and the disappearance of beneficial microflora in soils and the biodiversity that supports important processes in soils. Thus, the main goal of this study is to investigate how the ignorance of crop rotation rules for sunflowers affects the biophysical parameters of agricultural lands in Ukraine by analysis of satellite-based vegetation indices and crop classification information at the village council level.

Land Degradation Monitoring with Use of Satellite Data
The main purpose of modern agricultural practices is to achieve high crop productivity with a neutral level of land degradation. The implementation of these practices would allow farmers to use the full potential of agrarian lands without harming the environment and thus to overcome hunger and to ensure food security in the future.
In practice, land productivity indicators can be assessed with remote sensing imagery. In this case, integrated approaches are used to take into account the carbon content in the soil, trends in vegetation indices, and changes in land cover. All these indicators are interrelated-land use changes lead to a significant change in the carbon content of soils. This methodology is used to assess the sustainable development goals (SDGs) indicator 15.3.1, "Proportion of land that is degraded over total land area" [4]. In this case, the maps of changes in land cover are used as a separate sub-indicator that reflects the positive, neutral, or negative human impact on the land surface and as a means to assess changes in carbon content in the soil. The biophysical parameters of land productivity can be obtained from satellite images. Land productivity indicators estimated by remote sensing vegetation indices (VIs) are generally accepted for mapping and assessing land degradation and desertification [5]. They are widely used for countries in Europe [6] as well as in Asia [7] and Africa [8]. Different combinations of vegetation indices, which are studied as time series, are used and reflect the dynamics of vegetation and its biophysical indicators. The most popular indicator is the normalized difference vegetation index (NDVI), which is used for land productivity mapping for indicator 15.3.1 and for supporting the calculation of indicator 2.4.1, "Proportion of agricultural area under productive and sustainable agriculture" in a 2020 study by Kussul et al. [9]. The most common collection of satellite NDVI for land productivity maps is the MODIS data collection (MOD13Q1-coll6). However, today's satellite data processing methods and available satellite missions provide information on land productivity from higher spatial resolution NDVI. For example, the previously mentioned 2020 study by Kussul et al. [9,10] provides an example of using Sentinel-2 and Landsat-8 data to deliver land productivity maps with a spatial resolution of 30 m for Ukraine [9,11]. Such data enable us to analyze land degradation not only at the level of the country or administrative regions but also at the level of fields, thus improving the possibilities for decision making. Another vegetation index, the enhanced vegetation index (EVI), is more sensitive in areas with dense vegetation. Similar to NDVI, it can be used to quantify vegetation greenness. In addition to VIs, biophysical variables could be extracted from satellite data, including the leaf area index (LAI) and the fraction of absorbed photosynthetically active radiation (FAPAR). LAI is extracted from satellite data using complex biophysical models and better conveys the biophysical characteristics of the plant [12]. This index describes the amount of biomass on the earth's surface and its condition, and it makes it possible to qualitatively assess the crop yield and deliver land productivity maps [13]. Another biophysical variable, FAPAR, makes it possible to assess photosynthesis in plants and the absorption of solar energy. This biophysical characteristic directly indicates the primary productivity of photosynthesis. One more interesting vegetation index is the land surface water index (LSWI), which is commonly used for drought monitoring and reflects the total amount of water in vegetation and its soil background. The response of LSWI to rainfall indicates the usability of this index for water state monitoring for agricultural crops, especially in the critical time of the crops' development stages in the early part of the season [14]. In this article, we use all five of these indices to better analyze the biophysical impact of crop rotation on agricultural land drawing on their ability to reflect changes in the land surface from the perspective of various important factors of vegetation.

Crop Classification in Agricultural Monitoring
Most developed countries provide crop mapping based on satellite data on a regular basis or are working on the operationalization of crop mapping technologies. Accurate crop maps are in themselves a valuable product of data processing, and the availability of such products for one area for different years provides great opportunities for the qualitative analysis of land use practices. A good example of a large collection of historical crop maps is the United States' Cropland Data Layer, which has been making crop maps publicly available annually since 1997 [15]. The availability of long-term crop classification maps gives the possibility to get a better understanding of agricultural patterns for different territories and conduct very accurate forecasts. The research by Johnson D. et al. [16] shows that historical crop classification maps can be a valuable source of information for the pre-season or in-season crop classification. So, the pre-season crop classification map for the Corn Belt in the Unites States that was built with use of only previous classification maps can achieve 70% accuracy. In European countries that use the Land Parcel Identification System [17] as part of the Common Agriculture Policy, almost all crop information for agricultural parcels is provided by farmers; even in these cases, there is still a need to deliver such maps to analyze and validate the collected information. For this aim in Europe, there are systems that provide the automatic processing of satellite data and land cover and building of crop type maps, such as Sen-2-agri and Sen-4-CAP [18]. Thus, remote sensing data of the Copernicus program and machine learning land cover and crop type classification models are already having a lot of application in European Union for supporting CAP [19]. In 2021, the first European continental scale crop classification map based on Sentinel-1 data was published by the European Commission team [20]. The next steps will be the operationalization of continental scale crop classification technologies and the creation of crop maps for other years.
Ukraine is not a member of the European Union (EU) yet, and there are no products similar to the Land Parcel Identification System available for Ukraine. Nevertheless, during 2015-2020 we delivered crop type maps for Ukraine based on Sentinel-1 and Sentinel-2 data. For this purpose, we have used modern neural network approaches [21], LSTM and convolutional layers [22], and cloud technologies of high-performance computing on Google Earth Engine and Amazon platforms [23]. In this study, we take advantage of multiyear crop type maps available for Ukraine (2016-2020) at 10 m spatial resolution to identify regions with crop rotation violations and then upscale at village scale to estimate the impact on crop productivity.

Study Area
The study area in this research is Ukraine, which is one of the major producers and exporters of agricultural commodities. In Ukraine, 30 percent of nominal and 11.7 percent of real gross domestic product (GDP) are related to the agricultural sector. Since 1991 (the breakup of the Soviet Union), many alterations in land use have been made because of economic and policy changes as well as from the climate change impact. The past few years were especially dynamic in terms of land use changes and policy making due to (a) the military conflict in Eastern Ukraine and the occupation of 7 percent of the territory by Russian forces, (b) land market opening reform in Ukraine, and (c) the appearance of double-cropping fields due to the climate change warming processes. In terms of land productivity research, Ukraine is especially interesting because this country is a large world exporter of agricultural production. In 2020, Ukraine was listed as one of the five biggest agri-food exporters to the EU, with total exports to the EU of about EUR 1 billion [24]. In addition, Ukraine is the world's largest producer and exporter of sunflower oil, the world's third-largest exporter of corn, the fourth of barley, and the sixth of soybeans. The main crops grown in Ukraine are cereals, sunflowers, corn, soybeans, and sugar beets. Most cereal fields are related to the winter wheat fields, but also barley as a minoritarian crop is occurring. Despite the fact that Ukraine fully ensures its food security, it meets only a third of its potential [25]. The reason for this is the insufficient use of modern agricultural practices and irrigation. The biggest agro-climatic zone in Ukraine is the steppe zone, which occupies 40 percent of the country's territory (240,000 km 2 ). In this zone efficient agriculture is hardly possible without irrigation, but the total irrigated area in Ukraine in 2017 was 5000 km 2 [25].
Ukraine is actively struggling with problems of sustainable agriculture. These problems are compounded by global climate change processes that have significantly affected food production in the world and in Ukraine. According to Ukraine's 2021 Common Country Analysis designed by the United Nations [26], climate change significantly influences the effectiveness of the Ukrainian agricultural sector. The crop growth periods and crop calendars for the majority of crops are affected by the annual average temperature increase by 1.45 Celsius degrees over 50 years, which is double the global increase [27,28]. Particularly sensitive to these changes are winter crops, which have become more likely to sustain loss due to weather conditions. For example, in 2020, droughts in the Odessa region caused a 37 percent loss of winter crops; the loss of winter crops was observed on a smaller scale throughout Ukraine. Farmers had to plow over land, and thus Ukraine suffered great economic and food losses. Although there were no significant changes in cereal planting areas (50.2 percent in 2000 and 55.4 percent in 2012) observed by [29], there is a clear trend toward increasing the land area of industrial crops. As an example, from 2000 to 2014 the total area of sunflowers, soybeans, and rapeseed in Ukraine grew from 8.4 percent to 28.4 percent. In 2016, the area of sunflowers peaked, and it has remained stable for the past several years. The main reason for that stability is the use of sunflower seed as a sort of "insurance policy" against a high loss rate for winter crops [30], helping farmers to salvage income in such conditions. However, this option also leads to an increase in crop rotation violation events [31]. Similar patterns of agricultural practices are common in other Eastern European countries, such as Belarus and the Russian Federation. The planting of industrial crops leads to a decrease in the productivity of land and to its degradation. Therefore, the observed trend to increase the sown area of sunflowers contradicts the principles of sustainable agriculture in Ukraine.
To solve this problem, the government of Ukraine has introduced crop rotation rules that restrict farmers from planting industrial crops, forcing them to follow certain crop rotation schemes. For example, sunflowers can be planted on the same field once per seven years, according to the Resolution of the Cabinet of Ministers of Ukraine, 11 February 2010, Nr 164, "Approval of standards for optimal balance of crop types in crop rotation in different natural and agricultural regions" [32]. However, the introduction of such norms without the creation of a mechanism for detecting violations and for controlling crop rotations has Sustainability 2022, 14, 3965 5 of 23 not led to an improvement in the situation. The lack of a state instrument for crop rotation control still leads to frequent cases of crop rotation violations. In Ukraine, industrial crops such as sunflowers are sometimes planted not just twice in a row, but four or five times in a row. Additionally, farmers compensate for the reduction in soil productivity by increasing the use of fertilizers. Thus, according to Ukrainian State Statistics Service, from 2000 to 2014, the area of crops using fertilizers increased from 22 percent to 84 percent without changes in the use of organic fertilizer (2 percent to 3 percent). Furthermore, the continued industrialization of the agricultural sector leads to the increase in field sizes ( Figure 1). As a result, more than 91% of fields in Ukraine have an area bigger than 5 ha. Even 65% of minoritarian crops that occupy 6.8% of cropland are grown within the fields with an area of 5 or more ha, while most small fields (29% of total area for other classes) are usually owned by local country folk and located in villages close to the village housings. seven years, according to the Resolution of the Cabinet of Ministers of Ukraine, 11 February 2010, Nr 164, "Approval of standards for optimal balance of crop types in crop rotation in different natural and agricultural regions" [32]. However, the introduction of such norms without the creation of a mechanism for detecting violations and for controlling crop rotations has not led to an improvement in the situation. The lack of a state instrument for crop rotation control still leads to frequent cases of crop rotation violations. In Ukraine, industrial crops such as sunflowers are sometimes planted not just twice in a row, but four or five times in a row. Additionally, farmers compensate for the reduction in soil productivity by increasing the use of fertilizers. Thus, according to Ukrainian State Statistics Service, from 2000 to 2014, the area of crops using fertilizers increased from 22 percent to 84 percent without changes in the use of organic fertilizer (2 percent to 3 percent). Furthermore, the continued industrialization of the agricultural sector leads to the increase in field sizes ( Figure 1). As a result, more than 91% of fields in Ukraine have an area bigger than 5 ha. Even 65% of minoritarian crops that occupy 6.8% of cropland are grown within the fields with an area of 5 or more ha, while most small fields (29% of total area for other classes) are usually owned by local country folk and located in villages close to the village housings.

Crop Classification Maps
In this work, crop classification maps produced by Space Research Institute NASU-SSAU (SRI) based on the Sentinel-1 and Sentinel-2 satellite images with 10 m spatial resolution were used to evaluate crop rotations. These maps have been delivered for 2016-2020 using a crop classification methodology based on neural networks and satellite data time series [33]. This approach was used by an SRI team in the World Bank project

Crop Classification Maps
In this work, crop classification maps produced by Space Research Institute NASU-SSAU (SRI) based on the Sentinel-1 and Sentinel-2 satellite images with 10 m spatial resolution were used to evaluate crop rotations. These maps have been delivered for 2016-2020 using a crop classification methodology based on neural networks and satellite data time series [33]. This approach was used by an SRI team in the World Bank project "Supporting Transparent Land Governance in Ukraine", funded by the European Commission [34]. For model training and testing, we used two separate vector layers that contain polygons with different land cover and crop type classes. The ground truth information about crop type classes was collected during ground surveys along the roads throughout Ukraine for each year. Table 1 shows the number of polygons used for training (tr) and testing (tt) for each year. The number of training and testing polygons, manually labeled based on the ground truth points, are balanced. Figure 2 shows the geographical distribution of these polygons on the territory of Ukraine and the Sentinel-1 coverage used for these maps production.
Due to the fact that the data collection requires a lot of resources and time, it is not possible to cover the full territory uniformly. In addition, not all the roads can be used for the ground surveys due to the forest belts and distances to fields.    Table 2 shows the accuracy of classification for main land cover and agricultural classes. The accuracy of these maps for the main crops, such as cereals, maize, sunflowers, and soybeans, is high enough to complete accurate crop rotation analysis.   Table 2 shows the accuracy of classification for main land cover and agricultural classes. The accuracy of these maps for the main crops, such as cereals, maize, sunflowers, and soybeans, is high enough to complete accurate crop rotation analysis.

Vegetation Indices
To analyze crop conditions and the productivity of agricultural lands, it is possible to use five vegetation indices provided by MODIS satellite products with coarse spatial resolution. In particular, for research purposes we used: the normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI) at 250 m resolution and the leaf area index (LAI), fraction of absorbed photosynthetically active radiation (FAPAR), and land surface water index (LSWI) at 500 m resolution. Because the study has been conducted for several years, the numbers of measurements for the vegetation period (from the beginning of April until the end of October) were unified by usage of 16-day composites. The analysis ready data that are available as the 16-day composites were used for the NDVI and EVI estimation. The FAPAR and LSWI composites were obtained using a median filter on vegetation indices images with masked clouds. In this, a cloud-free collection of data with the same number of measurements for each year was produced. After this step, the cumulative vegetation indices (V I cum ) for the whole vegetation period were calculated with the formula: where DoY 16 is the set of days of the target year from April 1 until October 31 with the 16-day step and V I d+ is the vegetation index for the concrete day within the 16

Methodology
The main objective of this study is to investigate the impact of sunflower planting in the different crop rotation schemes during previous years on the biophysical parameters of crops for the current year. This investigation can be done by the analysis of relationships between the areas of various crop rotation schemes on the averaged vegetation indices at the village level. For this purpose, we will use the common technique for the multivariate dependences' analysis-regression analysis [35]. To do this, the dependent and independent variables should be determined first. The dependent variable should reflect the biophysical characteristics of agricultural land for some territory. The independent variables should reflect the impact of each certain crop rotation scheme to the value of the dependent variable. The linear regression model will describe the relations between dependent variable dv and independent variables {iv i } n i=1 : where i is the number of the crop rotation scheme, which takes values from 1 to n (total amount of crop rotation schemes), and a i and b are coefficients of regression. After assessing a i values for each independent variable they can be interpreted as the estimate of the impact of iv i on dv. If the sign of a i is negative, then the impact is negative; otherwise, the impact is positive. To evaluate the significance of the coefficient a i for the dv, we use the f-test of overall significance [36]. This test allows us to assess multiple coefficients simultaneously by checking the null hypothesis-that the fits of the model with and without specific coefficient are equal. For the f-test purpose, two models should be trained for each independent variable successively. If the independent variable is insignificant, the difference between the models' outputs with and without this variable will be small. As a result, the f-test provides the p-values for each regression coefficient-the probabilities dv would get the same or larger effect if the a i is equal to 0 (zero). So, the smaller p-value for a coefficient, the more likely iv i is significant.
First, the space of the features used for the regression analysis S and the dependent variables should be clarified. To do that the separate regression model for each vegetation index NDVI, EVI, LAI, FAPAR, and LSWI was created. The cumulative vegetation index V I cum , calculated by formula (1) and averaged by village boundaries for each village council index (ID) is proposed for use as a dependent variable. The independent variables (or features) should reflect the impact of each considered crop rotation on the vegetation index V I cum . To conduct that, the proportion of the area of each crop rotation and the area of fixed crop in the reference (target) year aggregated at the village council level were used as the features. In this way, the contribution of different crop rotations in the value of the vegetation index was considered. Thus, V I ID j,y is the cumulative value of the vegetation index VI for year y and crop type j, averaged for village ID. The respective vector of features S ID j,y contains the proportions of area of each considered crop rotation to the total area of fields with crop type j in the year y in the village ID. To calculate it, the reference year y and crop type j for the year y were fixed. In this way, n available crop rotations that are sequences of crop types for m years starting with the crop type j in the reference year were constructed. As a result, the vector S ID j,y consists of the following elements, calculated for each crop rotation i.
where S ID i is the area of crop rotation with number i in the village ID and S ID y (j) is the area of crop type j in the village ID and year y. If the crop rotation i will include all possible crop rotations for five years, the number of features, most of which could be present in the training set in only a few examples, will be enormous. This could lead to the overfitting of the model. In particular, analysis of available crop rotation combinations shows that 625 possible variations of crop rotations cover 80 percent of the country's area. In Ukraine, there are only 7834 villages, and none of them include all possible crop rotations. In this situation, any regression model would have the problem of overfitting, and the results will not be adequate. To address the issue, we propose to consider two different features' representation schemes.
The first one is representation of crop rotation by binary values for five years, in which j for the target year (2020) and each pixel for every previous year can have two possible values: j (or 1) or not_j (or 0). Such representation will be named as model A1. In this research, j corresponds to sunflowers (sf) and not_j to not sunflowers (nsf). So, there are 16 possible crop rotations (Table 3). We should mention that such representation is not exhaustive, because it does not take into account all possible combinations of crop rotations. For instance, crop rotation sf-nsf-nsf-nsf-nsf is not a violation of sunflower rotation rules, but in the years when sunflowers have not been planted, the rotation of other crops was not considered. It means that in the years when sunflowers were not planted, crop rotation violations for other crop types still could appear. Thus, this model reflects the impact of sunflower monocropping, but it does not measure the impact of different crop rotations.
To address the issue, taking into account different crop rotations and at the same time avoiding the problem of overfitting, let us consider a smaller number of features.
To reduce the number of crop rotation combinations for the analysis and at the same time increase the number of training points, crop rotation analysis for major crops only for three years with several different reference years will be considered. This representation will be named as model B.
Model B uses 25 crop rotations (25 features) for each value of cumulative VI for sunflowers per each village council. As a result, it will be possible to evaluate 25 possible crop rotations for each reference year and select the best and the worst crop rotations. Table 4

Regression Analysis of Crop Rotation
To analyze the impact of sunflower crop rotation we propose to use a model that represents the cumulative vegetation index for sunflowers, averaged for village ID, as a linear regression function (2) from the independent values S ID j,y,i (3). Calculation of coefficients a i could be considered as a traditional regression problem, which could be solved using an ordinary least squares approach. However, the large number of independent variables in the model leads to a collinearity problem and to potentially overfitting the model. To avoid it, the ridge regression technique was used [37]. The technique uses L2 regularization, which causes regression coefficients to be more balanced and representative in the regression. Thus, the regression model used in our study is this: where α is the L2 regularization penalty. To determine the best value of α, it is possible to use a cross validation technique that fits the regression function with α ∈ 10 k 10 The reliability of the results can be assessed with the use of the f-test of overall significance. In this case, three levels for the p-value for the regression coefficients were defined. The first is * p > 0.05, in which the variable is insignificant and the regression coefficient does not express the impact of independent variables to the dependent variable. The second is ** p > 0.01, in which the regression coefficient expresses the impact of the independent variable to the dependent variable; however, the significance is low. The third is *** p < 0.01, in which the variable is significant and the regression coefficient expresses the impact of the independent variable to the dependent variable.

Model A1 with Binary Crop Rotation Features and Analysis for Five Years
The first model-model A1-evaluates the impact of monocropping on the different vegetation indices. In this case, each crop rotation is described by the sequence of n binary values (we consider n − 1 preceding years), so for each i [1,16] and y [2016, 2020] cl i y = 1 if sunflowers had been grown in the territory in year y, and cl i y = 0 in the other case. Model A1 describes the relation of V I j,2020 at the village level from 16 crop possible monocropping combinations during the previous four years. The main drawback of this approach is the ignorance of the influence of other crop rotations on the vegetation indices.

Model A2 and A3 with Derived Regressors
Binary representation of monocropping schemes allows for the design of a regression model (model A2) of the crop rotation impact prediction for the i-th monocropping scheme (a i ) based on the number of sunflower monocropping (N s f ) for five years and the period of years without sunflowers before the reference year (P s f ): Model A2 uses a i obtained from model A1 as a dependent variable. So, if model A1 provides the evaluation of the crop rotation scheme represented as the sequence of sunflower and not sunflower plantings, model A2 fits to reproduce the result of model A1, but only with the use of statistical characteristics of sunflower appearance in the crop rotation scheme. The main advantage of this model is the small number of regressors, which allows us to avoid overfitting the model and to provide more accurate regression analysis. In addition, based on this model, it is possible to get the evaluation of crop rotation schemes that cannot be estimated from the available time-series of crop maps. Having a statistically significant model (5), we can mathematically estimate the optimal interval between the subsequent planting of sunflowers and can provide predictions (extrapolation) for a longer time period. The main problem in this model is that N s f and P s f are not unique characteristics of concrete crop rotation. For example, model A2 considers crop rotation with N s f = 3 and P s f = 1 as a crop rotation with sunflower planting every two years, so two different crop rotations such as sf-nsf-sf-nsf-sf (10101) and sf-nsf-sf-sf-nsf (10110) can have the same N s f and P s f values. In order to be sure of the reliability of the result, another model-model A3 was fitted. This model represents the relationships between a i and the time intervals between sunflower planting, in the case when all sunflower plantings in the crop rotation scheme have equal intervals. So, the A3 model can be considered a simplification of the A2 model that predicts a i with use of only one variable-the interval between two sunflower plantings in the crop rotation scheme. The fitting strategy is the same as for model A2-dependent values are a i , obtained from model A1, and independent values are P s f , from crop rotation schemes. In this set of points, there are only five schemes that can be used for fitting such a model. We assumed that a i for crop rotations with 0-,1-,2-, 3-, and 4-year planting intervals for 5 years will be same as for 11 years. So, the crop rotation sf-nsf-nsf-nsf-nsf (10000) will be considered as the planting of sunflower once per 5 years in 11 years of observation. This assumption can be done because the impact of crop rotations weakens over the years. So, the plantings made more than 5 years before have a much lower impact than recent plantings do. However, if the p-values of a i for some of the useful crop rotation schemes is too big, these crop rotation schemes cannot be used for the A3 model fitting.

Model B for Three-Year Crop Rotation Analysis out of Five Years
The second regression model-model B-takes into account all possible crop rotations during three years for the VIs of sunflower. In this, it is possible to increase the number of training samples by considering different three-year intervals of observations. So, as an output (dependable variable) of the model we consider cumulative VIs for sunflower for the years 2020, 2019, and 2018: V I j,2020 ∪ V I j,2019 ∪ V I j,2018 = V I j . Model B includes 25 different regressors; each of them corresponds to the sequence of two crop types grown in the preceding two years out of the five major crop types (cereals, sunflowers, maize, soybeans, and other crops). In this case, the example of crop rotation violation is sunflowersunflower-sunflower, and the reverse example is sunflower-soybeans-other crops.

Potential of Sentinel and Landsat Data Usage
MODIS data collection was chosen because of the need for high temporal resolution data that give the possibility to create a uniform collection of vegetation indices time series for the accurate estimation of accumulated vegetation indices for each year of interest. In our experiment, moderate spatial resolution 250 and 500 m is more than enough for crop monitoring at the village council level. However, in the future studies it is possible to improve the accuracy of the experiment by the use of higher spatial resolution data. The growth of available data source numbers as well as the development of new methods of satellite data collection harmonization can give opportunities to improve the described experiment for further years.
There are two more available satellite missions that can be used for such an experiment in the future. The first one is the Landsat mission that launched in 1972. The two most recent satellites of this mission have similar characteristics-Landsat-8 launched in February 2013 and Landsat-9 launched in September 2021. These satellites have two instruments-an Operational Land Imager (OLI) that provides 8-band images with 30 m spatial resolution and 1 panchromatic band with 15 m spatial resolution and a Thermal Infrared Sensor (TIRS) that provides 2 bands of thermal infrared specter with 100 m spatial resolution. The temporal resolution of both satellites is 16 days. The second is the Sentinel mission. The Sentinel-3A satellite was launched in the February 2016. Images obtained from the Ocean and Surface Colour Instrument, installed on the Sentinel-3 satellites, have 21 multispectral bands with 300 m spatial resolution, 1-day temporal resolution, and can be used as an alternative for MODIS data. The Sentinel-2 mission was launched in June 2015, and since March 2017 includes two satellites with multi-spectral instruments that provide 12-band images with a spatial resolution from 10 m to 60 m and temporal resolution of 5 days. Sentinel-2 optical data is an essential source of information in the agricultural monitoring applications that give the possibility to estimate all essential vegetation indices and biophysical characteristics of crops that can be used in the crop state analysis such as NDVI, EVI, LAI [38], and others.
A lot of applications in remote sensing require a dense time-series of measurements with high temporal resolution. Thus, harmonization of multi-satellite high spatial resolution data collections is the cornerstone for the improvement of land monitoring as well as agricultural monitoring systems. Claverie et al.'s study [39] shows the workflow that can be used for the estimation of the harmonized collection of Landsat-8 and Sentinel-2 data. The harmonization process requires the conductance of atmospheric correction, geometric resampling, geographic registration, BRDF normalization, and band pass adjustment. A good example of this harmonization method use is shown in the research by Skakun et al. [40] on yield forecasting in Ukraine. This method has a high potential for the use of 30 m optical data for crop rotation analysis in future research, despite the need for a lot of computational resources for data processing. However, the combined and harmonized archive data of Landsat-8 and Sentinel-2 still have a lower temporal resolution, if we take into account cloud conditions, in comparison with MODIS or Sentinel-3. Kirovohradska oblast, located in central Ukraine, on average for the winter crop vegetation period (from March until the end of June) has 6 cloud-free observations for 2016, 8 for 2017, and 11 for 2018 [41]. The absence of uniformity in the data coverage at the regional or country level causes complications in the use of the historical harmonized collection of Landsat-8 and Sentinel-2 data in the multi-year crop rotation analysis experiment. Such uniformity for these years can be explained not only by atmospheric conditions, but also by the fact that before the launch of Sentinel-2B satellite in March 2017, the temporal resolution of Sentinel-2 data was 10 days. This is why, in our experiment on the historical data from 2016 to 2020, we are using MODIS data. However, in future experiments from 2018 with 3 satellites (Landsat-8, Sentinel-2A and Sentinel-2B) or from 2022 with 4 satellites (by adding Landsat-9), the temporal resolution of harmonized high spatial resolution image collections can improve the results obtained by the methods described in the Methodology section.
As an additional data source in the crop rotation analysis, we can use the Sentinel-1 mission that was launched in the April 2014. The spatial resolution of Sentinel-1 is 20 m for the Ground Range Detected (GRD) data used for the land monitoring applications. Usage of Synthetic Aperture Radar data usually require the usage of complex data processing workflows that include calibration, geometric correction, terrain correction, and resampling to 10 m spatial resolution. Due to the presence of noise that seriously influence the quality of data, an additional filtering step with a refined Lee algorithm is required [42] to acquire qualitative data. As a result, it is possible to obtain SAR data with VV and VH polarization, 10 m spatial resolution, and 6-day temporal resolution that can be used for crop phenology estimation [42] and crop classification [35]. In addition, Filgueiras et al.'s [43] study showed high dependences between Sentinel-1 SAR indices and NDVI. This study shows that it is possible to model Sentinel-2 NDVI synthetic data with the usage of regression functions based on the VV and VH characteristics of Sentinel-1, expanding in this way the time-series of Sentinel-2 vegetation indices. If we consider that Sentinel-1 is an active sensor that is not vulnerable to the clouds, these synthetic data can be used for missing value recovery or as an alternative to absent Sentinel-2 images due to the atmospheric conditions.

Results for Model A1 with Binary Crop Rotation Features for Five Years
In this sub-section, we provide an analysis of 4 years of sunflower monocropping impact on the biophysical characteristics with the binary representation of crop rotation schemes. To evaluate this impact, we used a regression model A1 with 16 regressors (pre-dictors). Figure 4 (distribution of areas of each crop rotation scheme) demonstrates one of the challenges of such a study-a significant imbalance of crop rotation area representation. An increase in the time period for the analysis can lead to an increase in imbalance; in addition, very important crop rotation schemes, such as sunflower planting five, four, or even three years in a row cover only 1.1, 0.8, and 1.7 percent of the total sunflower area, respectively. It is a big problem because the fields with the highest frequencies of sunflower monocropping are the best sites for such analysis.

Results for Model A1 with Binary Crop Rotation Features for Five Years
In this sub-section, we provide an analysis of 4 years of sunflower monocropping impact on the biophysical characteristics with the binary representation of crop rotation schemes. To evaluate this impact, we used a regression model A1 with 16 regressors (pre dictors). Figure 4 (distribution of areas of each crop rotation scheme) demonstrates one o the challenges of such a study-a significant imbalance of crop rotation area representa tion. An increase in the time period for the analysis can lead to an increase in imbalance in addition, very important crop rotation schemes, such as sunflower planting five, four or even three years in a row cover only 1.1, 0.8, and 1.7 percent of the total sunflower area respectively. It is a big problem because the fields with the highest frequencies of sun flower monocropping are the best sites for such analysis.  (Table 3), and the Y axis the areas they occupy in the percentages.
After fitting the regression function, the regression coefficients for each monocrop ping scheme and each vegetation index were accessed (Table 5). To make it easier to ana lyze, these coefficients were numerated on the basis of the binary representation of the schemes-0 for sunflower and 1 for no sunflower. Such representation contains infor mation about the number of years of growing sunflowers within five years and the  (Table 3), and the Y-axis the areas they occupy in the percentages.
After fitting the regression function, the regression coefficients for each monocropping scheme and each vegetation index were accessed (Table 5). To make it easier to analyze, these coefficients were numerated on the basis of the binary representation of the schemes-0 for sunflower and 1 for no sunflower. Such representation contains information about the number of years of growing sunflowers within five years and the position (year) of the closest sunflower planting to the reference year (2020 or 2019). These characteristics can be mathematically obtained with use of logarithm function and subtractions with 5 (length of sequence) and 17 (group order). If the number is N, then the first value can be extracted from the N as 5 − [log 2 (N)] and the second as [5 − log 2 (17 − N)], where [ ] -is the integer part of the expression. In general, the lowest values of regression coefficients correspond to those schemes with monocropping of sunflowers. The biggest values correspond to those schemes in which during the previous two years sunflowers had not been planted. The worst three-crop rotations-those that have the lowest values for most of the biophysical crop parameters-are sf-sf-nsf-nsf-nsf (11000), sf-sf-sf-nsf-nsf (11100), and sf-sf-nsf-sf-nsf (11010). The best crop rotations from the prospective of all VIs and biophysical parameters are those with no more than one planting of sunflowers in five years. In general, it is possible to draw the following conclusion from these results: the best crop rotation schemes are those in which sunflowers are as rare as possible. Most schemes with violations of sunflower crop rotation in the past three years have negative coefficients with big absolute values. In that case, Table 5 shows that the impact depends on the year of violation and the number of repeated sunflower plantings.   The proposed binary representation gives the possibility to estimate the number of years of sunflower planting in each scheme and the number of years without sunflower planting before the reference year. The first value has a negative moderate correlation with a i equal to −0.74 and an R 2 score equal to 0.55. The second value has a strong correlation of 0.83 with a i and an R 2 score equal to 0.69. These results uncover the fact that the number of years without monocropping is a more important indicator than the number of years of sunflower plantings. In addition, we can conclude that the impact of the effect of crop rotation violations weakens over time and that a sunflower crop rotation violation four years ago has a much lower effect than a violation in recent years (one or two years ago). Figure 6 shows the dependences between the crop rotation scheme's areas and regression coefficients. It shows a moderate correlation (0.45) between the crop rotation scheme's area and regression coefficients. So, for the long-term crop rotation strategies, land users prefer to use crop rotation schemes with a neutral or positive effect on the vegetation indices, rather than schemes that lead to the degradation and decrease in the vegetation indices.

Results for A2 and A3 Models with Derived Regressors
The results from section 5.1 give the possibility to conduct modeling of other crop rotation schemes' impact on the land fertility that cannot be analyzed without long-term measurements and statistical significance of their areas. Using coefficients of model A1 (Table 5), it is possible to estimate coefficients of model A2 described by Equation (5) with the traditional least squares approach. After the model training we obtained the relationship (6) for coefficients for such crop rotations: This model has a 0.71 R 2 score and a very low mean squared error equal to 0.0052. This means that these two characteristics (number of sunflower monocroppings for five years and the period without sunflowers before reference year ) can be used to predict the effect of crop rotation on vegetation indexes.
After this, we used model A2 as a simplification of model A1, described by relationship (6), based on the and . For policy making, it is important to estimate the optimal interval between subsequent plantings of sunflowers. If we consider the time period of observations as 11 years, the period without sunflowers before reference year takes values from 0 to 10 and the number of sunflower monocroppings for five years = � 10 +1 � So, having and for each value in the interval from 0 to 10, it is possible to model the influence of monocropping schemes as a polynomial approximation of the output of model A2 (Figure 7a). Now, it is observable that with an increasing interval, the positive effect of the indicators increases, but the slope of the curve decreases. These observations are consistent with the results obtained from model A3 (Figure 7b) that trained based on the crop rotation schemes that reflect the planting of sunflower from once per 2 years to once per 5 years.

Results for A2 and A3 Models with Derived Regressors
The results from Section 5.1 give the possibility to conduct modeling of other crop rotation schemes' impact on the land fertility that cannot be analyzed without long-term measurements and statistical significance of their areas. Using coefficients a i of model A1 (Table 5), it is possible to estimate coefficients of model A2 described by Equation (5) with the traditional least squares approach. After the model training we obtained the relationship (6) for a i coefficients for such crop rotations: This model has a 0.71 R 2 score and a very low mean squared error equal to 0.0052. This means that these two characteristics (number of sunflower monocroppings for five years N s f and the period without sunflowers before reference year P s f ) can be used to predict the effect of crop rotation on vegetation indexes.
After this, we used model A2 as a simplification of model A1, described by relationship (6), based on the N s f and P s f . For policy making, it is important to estimate the optimal interval between subsequent plantings of sunflowers. If we consider the time period of observations as 11 years, the period without sunflowers before reference year P s f takes values from 0 to 10 and the number of sunflower monocroppings for five years N s f = 10 P s f +1 . So, having N s f and P s f for each value in the interval from 0 to 10, it is possible to model the influence of monocropping schemes as a polynomial approximation of the output of model A2 (Figure 7a). Now, it is observable that with an increasing interval, the positive effect of the indicators increases, but the slope of the curve decreases. These observations are consistent with the results obtained from model A3 (Figure 7b) that trained based on the crop rotation schemes that reflect the planting of sunflower from once per 2 years to once per 5 years. confidence interval. Therefore, for sustainable food production, the restriction can be reduced to planting sunflowers once every three years.

Results for Model B for Three Years of Crop Rotation out of Five Years
In this sub-section, we investigated the relationships between vegetation indices for sunflowers and crop rotation schemes over three years using the B model. To conduct this analysis, we calculated cumulative VIs for 2018, 2019, and 2020, which were selected as reference years for the dependent variables according to (1). The regressors for Equation (4) have been estimated by the same methodology as for model A, but for three years and all five possible crop types (cereals, maize, sunflowers, soybeans, and other crops) used in the crop rotation schemes. Data sets for all three reference years have been joined into one data set for training the model. The distribution of three years of crop rotation schemes ( Figure 8) is more preferable for our model than the distribution of model A because the crop rotations with sunflower plantings two or three times per three years (indices 3, 11-15, 18, and 23) are presented with a large number of samples. According to the crop rotation rules in Ukraine, sunflowers can be planted once per seven years. If the highest positive effect is observed for the 10-year interval, a 6-year interval is 60 percent of its value. At the same time, the negative effect on the vegetation characteristics disappears after the interval equals three, taking into account the confidence interval. Therefore, for sustainable food production, the restriction can be reduced to planting sunflowers once every three years.

Results for Model B for Three Years of Crop Rotation out of Five Years
In this sub-section, we investigated the relationships between vegetation indices for sunflowers and crop rotation schemes over three years using the B model. To conduct this analysis, we calculated cumulative VIs for 2018, 2019, and 2020, which were selected as reference years for the dependent variables according to (1). The regressors for Equation (4) have been estimated by the same methodology as for model A, but for three years and all five possible crop types (cereals, maize, sunflowers, soybeans, and other crops) used in the crop rotation schemes. Data sets for all three reference years have been joined into one data set for training the model. The distribution of three years of crop rotation schemes ( Figure 8) is more preferable for our model than the distribution of model A because the crop rotations with sunflower plantings two or three times per three years (indices 3, 11-15, 18, and 23) are presented with a large number of samples.  (Table 4), and the Y-axis the areas they occupy in the percentage. Table 6 represents the impact of crop rotations on the vegetation indices of sunflowers. The analysis shows that for sunflowers, the worst crop rotations are sunflower-sunflower-cereals, sunflower-cereals-cereals, sunflower-other-sunflower, and sunflowersunflower-sunflower. The best crop rotations are sunflower-soybeans-other and sunflower-other-other. For most of the vegetation indices, negative crop rotations include sunflowers planted at least once per two years or cereals crop rotation violations. This model is not sensitive to the violation of maize and soybeans. To mitigate the negative effect on sunflower production, farmers would be better to rotate it with maize or soybeans and, at the same time, plant sunflowers no more than once per three years.   (Table 4), and the Y-axis the areas they occupy in the percentage. Table 6 represents the impact of crop rotations on the vegetation indices of sunflowers. The analysis shows that for sunflowers, the worst crop rotations are sunflower-sunflowercereals, sunflower-cereals-cereals, sunflower-other-sunflower, and sunflower-sunflowersunflower. The best crop rotations are sunflower-soybeans-other and sunflower-other-other. For most of the vegetation indices, negative crop rotations include sunflowers planted at least once per two years or cereals crop rotation violations. This model is not sensitive to the violation of maize and soybeans. To mitigate the negative effect on sunflower production, farmers would be better to rotate it with maize or soybeans and, at the same time, plant sunflowers no more than once per three years. Figure 9 shows the dependences between the crop rotation scheme's areas and regression coefficients. The cumulative areas for three years, used for the experiment, reflect the tendency of crop rotation scheme occurrence from 2016 to 2020. Now, the correlation between the crop rotation scheme areas and regression coefficients is moderately negative (−0.46). So, for the short-term crop rotation strategy, land users prefer to use crop rotation schemes with a negative effect on the vegetation indices. The crop rotation schemes that occupy the biggest part of the land are sunflower-cereals-sunflower, sunflower-cerealsother, sunflower-cereals-cereals, and sunflower-maize-sunflower. So, we can see that crop rotations with sunflower planting rule violations are very popular. It can be explained by the economic benefits that can be obtained from the planting of sunflower or other industrial crops. Additionally, it is important to explain why, in long-term planning, the land users have a tendency to use more sustainable crop rotation schemes, while in the short-term, unsustainable are more preferable. Due to the use of the accumulated three-year intervals areas statistics, this tendency does not mean that most fields have negative crop rotation schemes every three years. It means that most fields have negative three-year crop rotation schemes at least once per 5 years. The negative effect of some three-year crop rotation schemes can be negated by further crop rotations on the fields.  Note: *** p<0.01. Figure 9 shows the dependences between the crop rotation scheme's areas and regression coefficients. The cumulative areas for three years, used for the experiment, reflect the tendency of crop rotation scheme occurrence from 2016 to 2020. Now, the correlation between the crop rotation scheme areas and regression coefficients is moderately negative (-0.46). So, for the short-term crop rotation strategy, land users prefer to use crop rotation schemes with a negative effect on the vegetation indices. The crop rotation schemes that occupy the biggest part of the land are sunflower-cereals-sunflower, sunflower-cerealsother, sunflower-cereals-cereals, and sunflower-maize-sunflower. So, we can see that crop rotations with sunflower planting rule violations are very popular. It can be explained by the economic benefits that can be obtained from the planting of sunflower or other industrial crops. Additionally, it is important to explain why, in long-term planning, the land users have a tendency to use more sustainable crop rotation schemes, while in the shortterm, unsustainable are more preferable. Due to the use of the accumulated three-year intervals areas statistics, this tendency does not mean that most fields have negative crop rotation schemes every three years. It means that most fields have negative three-year crop rotation schemes at least once per 5 years. The negative effect of some three-year crop rotation schemes can be negated by further crop rotations on the fields.  6. Discussion

Crop Classification Maps Time Series
The analysis of the impact of crop rotations on biophysical parameters of plants became possible because of the use of crop classification maps for 2016-2020. One of the important questions that arises when performing this analysis is whether five years of data are enough for the analysis. On the one hand, five years is not a very long period of time.
However, on the other hand, the modernization of production and the introduction of modern technologies in agriculture in Ukraine are carried out at a very fast pace. Therefore, the information extracted from a very long time series of measurements may be biased by climate change and the modernization of agricultural practices. Climate change strongly influences the biophysical parameters of plants growing in Ukraine and the nomenclature of crops and species. Taking into account these factors, five years is the optimal time for such modeling.
Crop rotation scheme representation was chosen as proportions of areas of each crop rotation at the village council level because it is the best way to represent crop rotation information in the linkable way for biophysical characteristics. The use of pixel or field level representation does not give the possibility to compare different crop rotation schemes because for each unit of analysis only one type of it will be presented. At the same time, the averaging of vegetation indices at village council level give the possibility to estimate which crop rotation, in the proportion, gives a higher increase or loss of vegetation indices.

Crop Classification Accuracy
To perform this study, we use crop classification maps delivered on the basis of time series of optical and radar satellite data from Copernicus Sentinel missions. As a classification algorithm, modern approaches for crop classification based on neural networks were used, which provide rather high accuracy of classification for different crops. Neural networks have been trained on the in situ data collected each year during ground surveys along the roads in Ukraine. The obtained crop classification maps have high accuracies with the F1-score from 81 percent up to 99 percent for majoritarian crop classes. Therefore, it can be concluded that the crop classification maps used in the analysis provide sufficient accuracy for this analysis.

Crop Rotation Valuation Drawbacks
The main problem with the proposed methodology is the influence of weather conditions and agricultural practices on the value of vegetation indices. These factors lead to the possibility that the values of vegetation indices for crop rotations with a low positive or negative effect may be biased. Thus, the sign of the regression coefficients does not show in general whether they have a positive effect on the vegetation, but they show that with such crop rotation under these conditions, the vegetation index decreases. However, it is still possible to determine which crop rotations are better and which are worse. In this case, the worst crop rotations would have the lowest values of regression coefficients, while the best crop rotations would have the highest values of regression coefficients.
Analyzing the results for concrete conditions in Ukraine, we can determine for which specific crop rotations there is an improvement or deterioration of biophysical indicators of vegetation in the climatic conditions of the country and in local agricultural practices.

Conclusions
In this study, we analyzed the impact of crop rotation on the biophysical parameters of sunflowers at the level of village councils using satellite data. The research method is regression analysis, which allows us to assess the impact of different indicators on the vegetation indices. As independent variables, we used the ratios of different crop rotation scheme areas for the villages.
The five-year monocropping model with the binary representation of sunflower planting confirms the negative effect of frequent planting of sunflowers. It shows the strong relationship between the number of years of sunflower plantings in the crop rotation and a decline in the vegetation indices. Additionally, it demonstrates that this effect is more significant for small intervals between sunflower plantings. This analysis also showed that planting sunflowers once per seven years according to the crop rotation rules in Ukraine is rational because of the high positive influence on the vegetation indices. However, starting from the sunflower planting interval equal to three years, the negative effect from sun-flower planting is not observed. Thus, planting sunflowers once every four years would be infrequent enough to avoid negative consequences and land degradation. Analysis of crop rotations of all major crops over three years shows that different crop rotations with and without violations have different effects on sunflower productivity. All crop rotations with previously planted sunflowers had the highest negative effect on the vegetation indices. The worst crop rotation with the highest negative effect for most crop types and biophysical characteristics is crop rotation violations for sunflowers. In addition, we found a tendency to use more sustainable crop rotation schemes for long-term crop rotation strategy and common violations of rules and use of unsustainable crop rotations for the short-term strategies. The results can be used for the implementation and modernization of agricultural norms in Ukraine and Europe. Funding: The authors acknowledge the funding received by the National Research Foundation of Ukraine from the state budget 2020/01.0273 "Intelligent models and methods for determining land degradation indicators based on satellite data" (NRFU Competition "Science for human security and society").

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.