Biophysical Impact of Sunflower Crop Rotation on Agricultural Fields

Abstract

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 soil fertility.

1. 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.

1.1. Crop Rotation Violation Impact on the Agricultural Productivity

Crop rotation violation leads to a significant change in the carbon content of soils. Currently, studies of the impact of agricultural management and of changes in plant species on agricultural lands is focused on soil quality. Interesting research in this area was made by Zuber et al. in 2018 [1]—the authors evaluate the effect of long-term crop rotation and tillage on the quantity of carbon and nitrogen stored in the soil organic matter. The study was conducted for the state of Illinois in the United States on the basis of in situ measurements on the test sites. It shows that some crop rotation schemes, such as corn-soybean-wheat, can increase the soil organic carbon (SOC) stocks in the soil, whereas crop rotation schemes such as soybean-soybean-soybean reduce SOC and total nitrogen content. Another research project was done for the Anhui province, China [2]. In the article, the authors simulated the effect of crop rotation using Fourier transformations based on the crop change patterns and vegetation indices. The result was validated using in situ measurements. The validation shows the high accuracy of the estimation of organic carbon in the topsoil based on a combination of vegetation variables with crop rotation types. This study showed the importance of crop rotation information and the effectiveness of its use to assess the change and state of organic carbon. A similar conclusion on the strong impact of crop rotation on carbon content can be drawn from a study [3] for Canada. In the article, the authors consider the important impact of crop rotation, the tillage system, and precipitation on the SOC stocks by using 30 years of observations. Looking at these results, the authors conclude that the plantation of pulse crops after wheat could be better for the SOC stocks than plantation of similar wheat crops.

1.2. 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.

1.3. 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.

2. 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²). In this zone efficient agriculture is hardly possible without irrigation, but the total irrigated area in Ukraine in 2017 was 5000 km² [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 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.

3. Data

3.1. 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.

Figure 3 shows these maps for major agriculture classes. The full time-series of Sentinel-1 images (15–16 images per year) with VV and VH polarization and cloud-free composites of Sentinel-2 images (4–5 composites per year) with red, green, blue, and near infrared (NIR) bands were used for each year. 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.

3.2. 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:

(1)$V{I}_{cum}={\displaystyle \sum }_{d\in Do{Y}_{16}}median{\left(V{I}_{d+i}\right)}_{r=1}^{16}$

where $Do{Y}_{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-day period. For instance, $Do{Y}_{16}$ for 2021 can be calculated as ${\left\{91+d\ast 16\right\}}_{d=0}^{13}$ (1 April 2021, is the 91st DoY). This information was extracted from MODIS data collections MOD13Q1 V6 for NDVI and EVI, MOD15A2H V6 for LAI and FAPAR and MOD09A1 V6 surface reflectance data for LSWI using Google Earth Engine platform. So, $V{I}_{cum}$ for each pixel is the sum of 16-day medians of correspondent VIs over the vegetation period of the target year.

4. 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 ${\left\{i{v}_i\right\}}_{i=1}^n$:

(2)$dv={\displaystyle \sum }_{i=1}^{n}i{v}_i\ast a_i+b$

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 $i{v}_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 $i{v}_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}_{j,y}^{ID}$ 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_{j,y}^{ID}$ 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_{j,y}^{ID}$ consists of the following elements, calculated for each crop rotation i.

(3)$S_{j,y,i}^{ID}=\frac{S_i^{ID}}{S_y^{ID}\left(j\right)}$

where $S_i^{ID}$ is the area of crop rotation with number i in the village ID and $S_y^{ID}\left(j\right)$ 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 illustrates the investigated crop rotation schemes in model B. This model combined the data for the years 2020-2019-2018, 2019-2018-2017, and 2018-2017-2016 in one training data set and used it for the regression function fitting. This approach, on the one hand, increases the number of training points and, on the other hand, reduces the effect of errors in the classification of extreme weather conditions for a specific year.

4.1. 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_{j,y,i}^{ID}$ (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:

(4)$dv={\displaystyle \sum }_{i=1}^{n}i{v}_i\ast a_i+b+\alpha \ast {\displaystyle \sum }_{i=1}^n{a}_i$

where $\alpha$ is the L2 regularization penalty. To determine the best value of $\alpha$, it is possible to use a cross validation technique that fits the regression function with $\alpha \in {\left\{{10}^k\right\}}_{k=-10}^{10}$ and selects the value $\alpha$ based on the best coefficient of determination (R² score). The $a_i$ coefficients defined through the model fitting determine the influence of crop rotation i in the biophysical characteristics of crops. The biggest values of $a_i$ correspond to the best crop rotations with positive impact on VI, while the lowest values of $a_i$ correspond to the worst crop rotations. As a result, it is possible to combine all crop rotations into three groups: those (a) with a negative effect on the vegetation, (b) with a positive effect on the vegetation, and (c) with a low effect on the vegetation. The usage of different vegetation indices allows us to evaluate this impact for different biophysical parameters of vegetation, such as NDVI, EVI, LAI, FAPAR, and LSWI.

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.

4.2. 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ϵ\left[1,16\right]$ and $yϵ\left[2016,2020\right]$ $c{l}_y^i=1$ if sunflowers had been grown in the territory in year y, and $c{l}_y^i=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.

4.3. 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_{sf})$ for five years and the period of years without sunflowers before the reference year ($P_{sf})$:

(5)$a_i=k1\ast N_{sf}+k2\ast P_{sf}+k3$

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_{sf}$ and $P_{sf}$ are not unique characteristics of concrete crop rotation. For example, model A2 considers crop rotation with $N_{sf}$ = 3 and $P_{sf}$ = 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_{sf}$ and $P_{sf}$ 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_{sf}$, 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.

4.4. 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}\cup V{I}_{j,2019}\cup 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 sunflower-sunflower-sunflower, and the reverse example is sunflower-soybeans-other crops.

4.5. 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.

5. Results

5.1. 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 (predictors). 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.

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-\left[lo{g}_2\left(N\right)\right]$ and the second as $\left[5-lo{g}_2\left(17-N\right)\right],$ 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.

Figure 5 shows the relationship between the position of crop rotation in the binary-sorted table and the values of regression coefficients. This figure does not include crop rotation schemes with numbers 1, 2, 3, 5, and 7, because they have not passed the f-test of overall significance. The trend of the values of regression coefficients for different vegetation indices is positive. It is possible to see that all schemes starting with sf-sf have the smallest values of the coefficients and are mostly negative. The Pearson correlation coefficient [44] between the position and average $a_i$ for five vegetation indices is 0.81, while the R² score [45] is 0.67. Those results mean that the regression coefficients are strongly correlated with the position of the scheme in the sorted list of binary representations.

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² score equal to 0.55. The second value has a strong correlation of 0.83 with $a_i$ and an R² 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.

5.2. 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:

(6)$a_i=-0.053\ast N_{sf}+0.071\ast P_{sf}+0.06$

This model has a 0.71 R² 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_{sf}$ and the period without sunflowers before reference year $P_{sf})$ 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_{sf}$ and $P_{sf}$. 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_{sf}$ takes values from 0 to 10 and the number of sunflower monocroppings for five years $N_{sf}=\left[\frac{10}{P_{sf}+1}\right]$.

So, having $N_{sf}$ and $P_{sf}$ 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.

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.

5.3. 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 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 sunflower-sunflower-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.

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-cereals-other, 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.

6. Discussion

6.1. 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.

6.2. 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.

6.3. 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.

7. 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 sunflower 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.

Author Contributions

Conceptualization N.K., K.D., L.S., M.L., D.A.A. and O.N.; methodology N.K., K.D., L.S., M.L. and D.A.A.; validation, L.S., M.L. and D.A.A.; investigation, L.S. and M.L.; data curation, N.K., K.D., L.S., M.L. and D.A.A.; writing—original draft preparation, N.K., L.S., M.L., K.D. and D.A.A.; writing—review and editing N.K., L.S., K.D. and D.A.A.; visualization, L.S.; supervision, N.K., K.D., D.A.A. and O.N. All authors have read and agreed to the published version of the manuscript.

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.

Data Availability Statement

This study is used publicly available datasets for the vegetation indices analysis: NDVI and EVI from the MOD13Q1 V6 collection (https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MOD13Q1 accessed on 25 January 2022), LAI and FAPAR from MOD15A2H V6 collection (https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MOD15A2H accessed on 25 January 2022) and MOD09A1 V6 surface reflectance data collection (https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MOD09A1 accessed on 25 January 2022)for LSWI calculation.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Major crop classes’ field sizes in Ukraine 2020.

Figure 2. Geographical distribution of training (tr) and testing (tt) polygons, collected during ground survey along the roads throughout Ukraine for each year and the coverage of territory by Sentinel-1 data.

Figure 3. Crop classification maps for major agricultural classes for Ukraine, 2016–2020.

Figure 4. The distribution of sunflower monocropping scheme areas in Ukraine for five years, used for the model A1 training. Note: The X-axis lists the classes of crop rotations (Table 3), and the Y-axis the areas they occupy in the percentages.

Figure 5. The dependences between the binary representation of monocropping scheme and regression coefficients for different vegetation indices.

Figure 6. The dependences between crop rotation scheme’s areas and regression coefficients for model A1.

Figure 7. The dependences between influence ([Forumla omitted. See PDF.]) of crop rotations and the interval between the sunflower planting in the corresponding crop rotations. (a) is the model A2 curve and (b) is the model A3 curve. Note: The x axis shows the intervals for the sunflower planting in the crop rotations. The y axis is the [Forumla omitted. See PDF.] value for the corresponding crop rotation. Points represent the training data for the model fitting; the curve is the fitted regression model based on these points. Shaded zones demonstrate the confidence intervals for prediction.

Figure 8. The distribution of the areas for different sunflower crop rotation schemes in Ukraine for three years, used for the training of model B. Note: The X-axis lists the classes of crop rotations (Table 4), and the Y-axis the areas they occupy in the percentage.

Figure 9. The dependences between crop rotation scheme’s areas and regression coefficients for the model B.

Table 1

The number of polygons, used for the model training (tr) and testing (tt) crop classification maps for 2016–2020.

Class	2016		2017		2018		2019		2020
	tr	tt	tr	tt	tr	tt	tr	tt	tr	tt
Artificial	68	67	68	67	164	164	222	223	244	237
Bare land	57	57	54	54	81	80	135	115	138	122
Cereals	1108	1106	586	586	1146	1146	1797	1776	1911	1896
Forest	338	338	342	341	666	665	223	224	297	291
Grassland	399	399	272	272	651	650	888	866	965	957
Maize	514	513	344	344	486	486	576	566	855	850
Peas	16	16	25	24	38	38	37	34	61	55
Soybeans	273	272	191	190	325	324	373	367	389	382
Sugar beet	30	29	20	19	69	69	36	29	36	30
Sunflower	823	822	596	595	876	875	690	681	794	790
Water	88	88	96	96	240	240	351	327	350	340
Wetland	24	24	35	35	66	66	107	114	133	119
Winter rapeseed	94	94	107	106	196	195	288	278	273	265
Total	3832	3825	2736	2729	5004	4998	5723	5600	6446	6334

Table 2

Accuracies of crop classification maps on the independent test sets for 2016–2020 (UA—User Accuracy, PA—Producer Accuracy, F1—F1 score, OA—Overall Accuracy).

Class	Classification Accuracies (%)
	2016			2017			2018			2019			2020
	UA	PA	F1	UA	PA	F1	UA	PA	F1	UA	PA	F1	UA	PA	F1
Artificial	73.1	77.2	75.0	69.0	71.7	70.3	66.6	95.0	78.3	66.4	96.4	78.7	71.5	88.5	79.1
Bare land	58.2	63.5	60.0	55.9	66.6	60.8	90.7	67.5	77.4	61.0	67.6	64.1	84.1	60.6	70.4
Cereals	92.0	88.1	90.0	97.3	96.5	96.9	99.0	98.7	98.8	99.1	98.5	98.8	94.5	92.2	93.3
Forest	96.9	99.3	98.1	99.2	98.2	98.7	99.7	99.8	99.8	83.4	99.2	90.6	85.5	99.8	92.1
Grassland	72.9	89.1	80.2	76.2	94.2	84.3	93.5	95.9	94.7	73.6	96.7	83.6	80.5	95.8	87.5
Maize	93.4	92.6	93.0	90.3	94.1	92.1	97.7	96.2	96.9	98.4	97.6	98.0	95.4	96.8	96.1
Peas	59.0	88.8	70.9	96.8	87.2	91.8	91.4	94.5	92.9	88.4	64.1	74.3	80.7	91.6	85.8
Soybeans	81.6	83.4	82.5	92.5	73.0	81.6	90.5	94.2	92.3	93.1	93.6	93.3	87.8	83.9	85.8
Sugar beet	90.9	96.6	93.6	99.8	97.5	98.6	96.4	94.6	95.5	96.6	99.3	97.9	93.3	78.7	85.4
Sunflower	94.3	94.3	94.3	94.2	97.1	95.6	98.0	98.1	98.0	98.7	98.8	98.7	96.8	97.0	96.9
Water	80.2	98.9	88.6	99.7	98.7	99.2	99.9	99.9	99.0	99.0	98.2	99.1	99.0	99.0	99.0
Wetland	77.1	77.6	77.3	83.8	49.8	62.5	85.1	64.4	73.3	88.5	81.6	84.9	86.2	86.4	86.3
Winter rapeseed	91.7	76.8	83.6	95.8	98.3	97.0	99.7	98.5	99.1	98.1	99.1	98.1	96.7	94.5	95.6
OA	88.3			91.0			97.3			96.2			91.6

Table 3

All crop rotations considered for model A1 (sf—sunflower, nsf—not sunflower).

Index	Binary Representation	2020	2019	2018	2017	2016
1	11111	sf	sf	sf	sf	sf
2	11110	sf	sf	sf	sf	nsf
3	11101	sf	sf	sf	nsf	sf
4	11100	sf	sf	sf	nsf	nsf
5	11011	sf	sf	nsf	sf	sf
6	11010	sf	sf	nsf	sf	nsf
7	11001	sf	sf	nsf	nsf	sf
8	11000	sf	sf	nsf	nsf	nsf
9	10111	sf	nsf	sf	sf	sf
10	10110	sf	nsf	sf	sf	nsf
11	10101	sf	nsf	sf	nsf	sf
12	10100	sf	nsf	sf	nsf	nsf
13	10011	sf	nsf	nsf	sf	sf
14	10010	sf	nsf	nsf	sf	nsf
15	10001	sf	nsf	nsf	nsf	sf
16	10000	sf	nsf	nsf	nsf	nsf

Table 4

All crop rotations considered for model B.

Index	2020\2019\2018	2019\2018\2017	2018\2017\2016
1	Sunflower	Cereals	Cereals
2	Sunflower	Cereals	Maize
3	Sunflower	Cereals	Sunflower
4	Sunflower	Cereals	Soybeans
5	Sunflower	Cereals	Other
6	Sunflower	Maize	Cereals
7	Sunflower	Maize	Maize
8	Sunflower	Maize	Sunflower
9	Sunflower	Maize	Soybeans
10	Sunflower	Maize	Other
11	Sunflower	Sunflower	Cereals
12	Sunflower	Sunflower	Maize
13	Sunflower	Sunflower	Sunflower
14	Sunflower	Sunflower	Soybeans
15	Sunflower	Sunflower	Other
16	Sunflower	Soybeans	Cereals
17	Sunflower	Soybeans	Maize
18	Sunflower	Soybeans	Sunflower
19	Sunflower	Soybeans	Soybeans
20	Sunflower	Soybeans	Other
21	Sunflower	Other	Cereals
22	Sunflower	Other	Maize
23	Sunflower	Other	Sunflower
24	Sunflower	Other	Soybeans
25	Sunflower	Other	Other

Table 5

The effect of sunflower monocropping on the vegetation indices from model A1.

Binary Representation	Crop Rotation	NDVI	EVI	FAPAR	LAI	LSWI
11111	sf-sf-sf-sf-sf	0.041 *	−0.032 *	−0.09 *	−0.047 *	−0.068 *
11110	sf-sf-sf-sf-nsf	−0.026 *	−0.158 *	−0.176 *	−0.065 *	−0.132 *
11101	sf-sf-sf-nsf-sf	−0.025 **	−0.097 *	−0.1 *	−0.048 *	−0.08 *
11100	sf-sf-sf-nsf-nsf	−0.232 ***	−0.194 ***	−0.121 ***	−0.046 ***	−0.046 ***
11011	sf-sf-nsf-sf-sf	−0.244 *	−0.192 *	−0.193 *	−0.105 *	−0.193 *
11010	sf-sf-nsf-sf-nsf	−0.195 ***	−0.09 ***	−0.06 ***	−0.115 ***	0.046 ***
11001	sf-sf-nsf-nsf-sf	−0.04 **	−0.143 *	−0.152 **	−0.087 *	−0.066 **
11000	sf-sf-nsf-nsf-nsf	−0.265 ***	−0.245 ***	−0.213 ***	−0.07 ***	−0.177 ***
10111	sf-nsf-sf-sf-sf	0.168 ***	0.035 **	0.126 **	0.009 **	0.007 ***
10110	sf-nsf-sf-sf-nsf	0.163 ***	0.221 ***	0.098 ***	0.022 ***	0.077 ***
10101	sf-nsf-sf-nsf-sf	−0.035 ***	0.006 ***	−0.018 ***	−0.028 ***	−0.042 ***
10100	sf-nsf-sf-nsf-nsf	0.144 ***	0.227 ***	0.206 ***	0.126 ***	0.19 ***
10011	sf-nsf-nsf-sf-sf	0.169 ***	0.037 ***	0.156 ***	0.023 ***	0.082 ***
10010	sf-nsf-nsf-sf-nsf	0.055 ***	0.104 ***	0.089 ***	0.101 ***	0.014 ***
10001	sf-nsf-nsf-nsf-sf	0.051 ***	0.177 ***	0.13 ***	0.094 ***	0.124 ***
10000	sf-nsf-nsf-nsf-nsf	0.27 ***	0.343 ***	0.319 ***	0.238 ***	0.263 ***

Note: * p > 0.05, ** p > 0.01, *** p < 0.01.

Table 6

The crop rotations’ effect on sunflower vegetation indices.

Crop Rotation	NDVI	EVI	LAI	FAPAR	LSWI
Sunflower-Cereals-Cereals	−0.254 ***	−0.249 ***	−0.145 ***	−0.21 ***	−0.187 ***
Sunflower-Cereals-Maize	−0.039 ***	−0.031 ***	−0.028 ***	−0.085 ***	−0.098 ***
Sunflower-Cereals-Sunflower	−0.194 ***	−0.22 ***	−0.166 ***	−0.208 ***	−0.133 ***
Sunflower-Cereals-Soybeans	0.16 ***	0.24 ***	0.124 ***	0.112 ***	0.057 ***
Sunflower-Cereals-Other	0.003 ***	0.075 ***	−0.009 ***	−0.009 ***	0.056 ***
Sunflower-Maize-Cereals	0.012 ***	0.056 ***	0.048 ***	0.01 ***	0.074 ***
Sunflower-Maize-Maize	−0.016 ***	0.064 ***	0.053 ***	−0.027 ***	−0.035 ***
Sunflower-Maize-Sunflower	−0.108 ***	−0.059 ***	−0.016 ***	−0.067 ***	−0.074 ***
Sunflower-Maize-Soybeans	0.071 ***	0.13 ***	0.062 ***	0.009 ***	−0.012 ***
Sunflower-Maize-Other	0.044 ***	0.055 ***	0.069 ***	0.047 ***	−0.008 ***
Sunflower-Sunflower-Cereals	−0.319 ***	−0.431 ***	−0.289 ***	−0.348 ***	−0.371 ***
Sunflower-Sunflower-Maize	0.059 ***	0.061 ***	−0.084 ***	−0.083 ***	−0.028 ***
Sunflower-Sunflower-Sunflower	−0.134 ***	−0.313 ***	−0.213 ***	−0.281 ***	−0.336 ***
Sunflower-Sunflower-Soybeans	−0.061 ***	−0.021 ***	0.025 ***	0.029 ***	0.073 ***
Sunflower-Sunflower-Other	0.005 ***	−0.039 ***	0.007 ***	0.047 ***	0.004 ***
Sunflower-Soybeans-Cereals	−0.012 ***	0.048 ***	0.027 ***	0.011 ***	−0.033 ***
Sunflower-Soybeans-Maize	0.047 ***	0.098 ***	0.199 ***	0.171 ***	0.147 ***
Sunflower-Soybeans-Sunflower	0.052 ***	0.07 ***	−0.014 ***	0.089 ***	0.115 ***
Sunflower-Soybeans-Soybeans	0.03 ***	0.088 ***	0.025 ***	0.02 ***	0.052 ***
Sunflower-Soybeans-Other	0.15 ***	0.137 ***	0.115 ***	0.206 ***	0.118 ***
Sunflower-Other-Cereals	−0.034 ***	0.026 ***	−0.009 ***	−0.027 ***	−0.063 ***
Sunflower-Other-Maize	−0.008 ***	0.03 ***	0.002 ***	0.065 ***	0.063 ***
Sunflower-Other-Sunflower	−0.219 ***	−0.224 ***	−0.157 ***	−0.173 ***	−0.133 ***
Sunflower-Other-Soybeans	0.065 ***	0.087 ***	0.049 ***	0.072 ***	0.018 ***
Sunflower-Other-Other	0.131 ***	0.112 ***	0.158 ***	0.195 ***	0.105 ***
R2 score	0.31	0.39	0.44	0.4	0.21

Note: *** p < 0.01.

References

1. Zuber, S.M.; Behnke, G.D.; Nafziger, E.D.; Villamil, M.B. Carbon and nitrogen content of soil organic matter and microbial biomass under long-term crop rotation and tillage in Illinois, USA. Agriculture; 2018; 8, 37. [DOI: https://dx.doi.org/10.3390/agriculture8030037]

2. Yang, L.; Song, M.; Zhu, A.-X.; Qin, C.; Zhou, C.; Qi, F.; Li, X.; Chen, Z.; Gao, B. Predicting soil organic carbon content in croplands using crop rotation and Fourier transform decomposed variables. Geoderma; 2019; 340, pp. 289-302. [DOI: https://dx.doi.org/10.1016/j.geoderma.2019.01.015]

3. Maillard, É.; McConkey, B.G.; Luce, M.S.; Angers, D.A.; Fan, J. Crop rotation, tillage system, and precipitation regime effects on soil carbon stocks over 1 to 30 years in Saskatchewan, Canada. Soil Tillage Res.; 2018; 177, pp. 97-104. [DOI: https://dx.doi.org/10.1016/j.still.2017.12.001]

4. McCallum, I.; Montzka, C.; Bayat, B.; Kollet, S.; Kolotii, A.; Kussul, N.; Lavreniuk, M.; Lehmann, A.; Maso, J.; Mazzetti, P. et al. Developing food, water and energy nexus workflows. Int. J. Digit. Earth; 2020; 13, pp. 299-308. [DOI: https://dx.doi.org/10.1080/17538947.2019.1626921]

5. Bai, Z.G.; Dent, D.L.; Olsson, L.; Schaepman, M.E. Proxy global assessment of land degradation. Soil Use Manag.; 2008; 24, pp. 223-234. [DOI: https://dx.doi.org/10.1111/j.1475-2743.2008.00169.x]

6. Kussul, N.; Kolotii, A.; Shelestov, A.; Yailymov, B.; Lavreniuk, M. Land degradation estimation from global and national satellite based datasets within UN program. Proceedings of the 2017 9th IEEE International Conference on Intelligent Data Acquisition and Advanced Computing Systems: Technology and Applications (IDAACS); Bucharest, Romania, 21–23 September 2017; Volume 1, pp. 383-386.

7. Wang, T.; Giuliani, G.; Lehmann, A.; Jiang, Y.; Shao, X.; Li, L.; Zhao, H. Supporting SDG 15, Life on Land: Identifying the Main Drivers of Land Degradation in Honghe Prefecture, China, between 2005 and 2015. ISPRS Int. J. Geo-Inf.; 2020; 9, 710. [DOI: https://dx.doi.org/10.3390/ijgi9120710]

8. Nhemachena, C.; Matchaya, G.; Nhemachena, C.R.; Karuaihe, S.; Muchara, B.; Nhlengethwa, S. Measuring baseline agriculture-related sustainable development goals index for Southern Africa. Sustainability; 2018; 10, 849. [DOI: https://dx.doi.org/10.3390/su10030849]

9. Kussul, N.; Lavreniuk, M.; Kolotii, A.; Skakun, S.; Rakoid, O.; Shumilo, L. A workflow for Sustainable Development Goals indicators assessment based on high-resolution satellite data. Int. J. Digit. Earth; 2020; 13, pp. 309-321. [DOI: https://dx.doi.org/10.1080/17538947.2019.1610807]

10. Dubovyk, O. The role of Remote Sensing in land degradation assessments: Opportunities and challenges. Eur. J. Remote Sens.; 2017; 50, pp. 601-613. [DOI: https://dx.doi.org/10.1080/22797254.2017.1378926]

11. Sims, N.C.; England, J.R.; Newnham, G.J.; Alexander, S.; Green, C.; Minelli, S.; Held, A. Developing good practice guidance for estimating land degradation in the context of the United Nations Sustainable Development Goals. Environ. Sci. Policy; 2019; 92, pp. 349-355. [DOI: https://dx.doi.org/10.1016/j.envsci.2018.10.014]

12. Hosseini, M.; McNairn, H.; Mitchell, S.; Robertson, L.D.; Davidson, A.; Ahmadian, N.; Bhattacharya, A.; Borg, E.; Conrad, C.; Dabrowska-Zielinska, K. et al. A comparison between support vector machine and water cloud model for estimating crop leaf area index. Remote Sens.; 2021; 13, 1348. [DOI: https://dx.doi.org/10.3390/rs13071348]

13. Kussul, N.; Lavreniuk, M.; Shumilo, L.; Kolotii, A. Nexus Approach for Calculating SDG Indicator 2.4. 1 Using Remote Sensing and Biophysical Modeling. Proceedings of the IGARSS 2019-2019 IEEE International Geoscience and Remote Sensing Symposium; Yokohama, Japan, 28 July–2 August 2019; pp. 6425-6428.

14. Chandrasekar, K.; Sesha Sai, M.V.R.; Roy, P.S.; Dwevedi, R.S. Land Surface Water Index (LSWI) response to rainfall and NDVI using the MODIS Vegetation Index product. Int. J. Remote Sens.; 2010; 31, pp. 3987-4005. [DOI: https://dx.doi.org/10.1080/01431160802575653]

15. Han, W.; Yang, Z.; Di, L.; Mueller, R. CropScape: A Web service based application for exploring and disseminating US conterminous geospatial cropland data products for decision support. Comput. Electron. Agric.; 2012; 84, pp. 111-123. [DOI: https://dx.doi.org/10.1016/j.compag.2012.03.005]

16. Johnson, D.M.; Mueller, R. Pre-and within-season crop type classification trained with archival land cover information. Remote Sens. Environ.; 2021; 264, 112576. [DOI: https://dx.doi.org/10.1016/j.rse.2021.112576]

17. Inan, H.I.; Sagris, V.; Devos, W.; Milenov, P.; Van Oosterom, P.; Zevenbergen, J. Data model for the collaboration between land administration systems and agricultural land parcel identification systems. J. Environ. Manag.; 2010; 91, pp. 2440-2454. [DOI: https://dx.doi.org/10.1016/j.jenvman.2010.06.030]

18. Defourny, P.; Bontemps, S.; Bellemans, N.; Cara, C.; Dedieu, G.; Guzzonato, E.; Hagolle, O.; Inglada, J.; Nicola, L.; Rabaute, T. et al. Near real-time agriculture monitoring at national scale at parcel resolution: Performance assessment of the Sen2-Agri automated system in various cropping systems around the world. Remote Sens. Environ.; 2019; 221, pp. 551-568. [DOI: https://dx.doi.org/10.1016/j.rse.2018.11.007]

19. Sarvia, F.; Xausa, E.; Petris, S.D.; Cantamessa, G.; Borgogno-Mondino, E. A Possible Role of Copernicus Sentinel-2 Data to Support Common Agricultural Policy Controls in Agriculture. Agronomy; 2021; 11, 110. [DOI: https://dx.doi.org/10.3390/agronomy11010110]

20. d’Andrimont, R.; Verhegghen, A.; Lemoine, G.; Kempeneers, P.; Meroni, M.; van der Velde, M. From parcel to continental scale--A first European crop type map based on Sentinel-1 and LUCAS Copernicus in-situ observations. Remote Sens. Environ.; 2021; 266, 112708. [DOI: https://dx.doi.org/10.1016/j.rse.2021.112708]

21. Kussul, N.; Lavreniuk, M.; Skakun, S.; Shelestov, A. Deep learning classification of land cover and crop types using remote sensing data. IEEE Geosci. Remote Sens. Lett.; 2017; 14, pp. 778-782. [DOI: https://dx.doi.org/10.1109/LGRS.2017.2681128]

22. Kussul, N.; Lavreniuk, M.; Shumilo, L. Deep Recurrent Neural Network for Crop Classification Task Based on Sentinel-1 and Sentinel-2 Imagery. Proceedings of the IGARSS 2020–2020 IEEE International Geoscience and Remote Sensing Symposium; Waikoloa, HI, USA, 26 September–2 October 2020; pp. 6914-6917.

23. Shelestov, A.; Lavreniuk, M.; Kussul, N.; Novikov, A.; Skakun, S. Exploring Google Earth Engine platform for big data processing: Classification of multi-temporal satellite imagery for crop mapping. Front. Earth Sci.; 2017; 5, 17. [DOI: https://dx.doi.org/10.3389/feart.2017.00017]

24. Monitoring EU Agri-Food Trade: Developments January–May 2020. Available online: https://ec.europa.eu/info/sites/default/files/food-farming-fisheries/trade/documents/monitoring-agri-food-trade_may2020_en.pdf (accessed on 15 February 2022).

25. Irrigation and Drainage Strategy in Ukraine until 2030. Available online: http://www.fao.org/faolex/results/details/en/c/LEX-FAOC190984/ (accessed on 15 February 2022).

26. Ukraine Common Country Analysis. Available online: https://ukraine.un.org/sites/default/files/2021-05/UN%20CCA%20Ukraine_April%202021%20%281%29.pdf (accessed on 15 February 2022).

27. Kogan, F.; Adamenko, T.; Guo, W. Global and regional drought dynamics in the climate warming era. Remote Sens. Lett.; 2013; 4, pp. 364-372. [DOI: https://dx.doi.org/10.1080/2150704X.2012.736033]

28. Kogan, F.; Kussul, N.; Adamenko, T.; Skakun, S.; Kravchenko, O.; Kryvobok, O.; Shelestov, A.; Kolotii, A.; Kussul, O.; Lavrenyuk, A. Winter wheat yield forecasting in Ukraine based on Earth observation, meteorological data and biophysical models. Int. J. Appl. Earth Obs. Geoinf.; 2013; 23, pp. 192-203. [DOI: https://dx.doi.org/10.1016/j.jag.2013.01.002]

29. Fileccia, T.; Guadagni, M.; Hovhera, V.; Bernoux, M. Ukraine: Soil Fertility to Strengthen Climate Resilience; Food and Agriculture Organization of the United Nations: Rome, Italy, Available online: http://www.fao.org/3/a-i3905e.pdf (accessed on 15 February 2022).

30. Sobolev, D. Oilseeds and Products Annual. USDA FAS: Kyiv, Ukraine, Available online: https://apps.fas.usda.gov/newgainapi/api/Report/DownloadReportByFileName?fileName=Oilseeds%20and%20Products%20Annual_Kyiv_Ukraine_04-15-2020 (accessed on 15 February 2022).

31. Kussul, N.; Lavreniuk, M.; Shelestov, A.; Skakun, S. Crop inventory at regional scale in Ukraine: Developing in season and end of season crop maps with multi-temporal optical and SAR satellite imagery. Eur. J. Remote Sens.; 2018; 51, pp. 627-636. [DOI: https://dx.doi.org/10.1080/22797254.2018.1454265]

32. Resolution of the Cabinet of Ministers of Ukraine, 11 February 2010, Nr 164. Available online: https://zakon.rada.gov.ua/laws/show/164-2010-%D0%BF#Text (accessed on 15 February 2022).

33. Shelestov, A.; Lavreniuk, M.; Vasiliev, V.; Shumilo, L.; Kolotii, A.; Yailymov, B.; Kussul, H.; Yailymova, H. Cloud approach to automated crop classification using Sentinel-1 imagery. IEEE Trans. Big Data; 2019; 6, pp. 572-582. [DOI: https://dx.doi.org/10.1109/TBDATA.2019.2940237]

34. Kussul, N.; Shelestov, A.; Yailymova, H.; Yailymov, B.; Lavreniuk, M.; Ilyashenko, M. Satellite Agricultural Monitoring in Ukraine at Country Level: World Bank Project. Proceedings of the IGARSS 2020-2020 IEEE International Geoscience and Remote Sensing Symposium; Waikoloa, HI, USA, 26 September–2 October 2020; pp. 1050-1053.

35. Cox, D.R.; Wermuth, N. Multivariate Dependencies: Models, Analysis and Interpretation. Chapman and Hall/CRC: London, UK, 2014.

36. Sureiman, O.; Mangera, C.M. F-test of overall significance in regression analysis simplified. J. Pract. Cardiovasc. Sci.; 2020; 6, 116. [DOI: https://dx.doi.org/10.4103/jpcs.jpcs_18_20]

37. Marquardt, D.W.; Snee, R.D. Ridge regression in practice. Am. Stat.; 1975; 29, pp. 3-20.

38. Pasqualotto, N.; Delegido, J.; Van Wittenberghe, S.; Rinaldi, M.; Moreno, J. Multi-crop green LAI estimation with a new simple Sentinel-2 LAI Index (SeLI). Sensors; 2019; 19, 904. [DOI: https://dx.doi.org/10.3390/s19040904] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30795571]

39. Claverie, M.; Ju, J.; Masek, J.G.; Dungan, J.L.; Vermote, E.F.; Roger, J.-C.; Skakun, S.V.; Justice, C. The Harmonized Landsat and Sentinel-2 surface reflectance data set. Remote Sens. Environ.; 2018; 219, pp. 145-161. [DOI: https://dx.doi.org/10.1016/j.rse.2018.09.002]

40. Skakun, S.; Vermote, E.; Franch, B.; Roger, J.C.; Kussul, N.; Ju, J.; Masek, J. Winter wheat yield assessment from Landsat 8 and Sentinel-2 data: Incorporating surface reflectance, through phenological fitting, into regression yield models. Remote Sens.; 2019; 11, 1768. [DOI: https://dx.doi.org/10.3390/rs11151768]

41. Rubel, O.; Lukin, V.; Rubel, A.; Egiazarian, K. Prediction of Lee filter performance for Sentinel-1 SAR images. Electron. Imaging; 2020; 9, 371. [DOI: https://dx.doi.org/10.2352/ISSN.2470-1173.2020.9.IQSP-371]

42. Veloso, A.; Mermoz, S.; Bouvet, A.; Le Toan, T.; Planells, M.; Dejoux, J.F.; Ceschia, E. Understanding the temporal behavior of crops using Sentinel-1 and Sentinel-2-like data for agricultural applications. Remote Sens. Environ.; 2017; 199, pp. 415-426. [DOI: https://dx.doi.org/10.1016/j.rse.2017.07.015]

43. Filgueiras, R.; Mantovani, E.C.; Althoff, D.; Fernandes Filho, E.I.; Cunha, F.F.D. Crop NDVI monitoring based on sentinel 1. Remote Sens.; 2019; 11, 1441. [DOI: https://dx.doi.org/10.3390/rs11121441]

44. Benesty, J.; Chen, J.; Huang, Y.; Cohen, I. Pearson correlation coefficient. Noise Reduction in Speech Processing; Springer: Berlin/Heidelberg, Germany, 2009; pp. 1-4.

45. Helland, I.S. On the interpretation and use of R2 in regression analysis. Biometrics; 1987; 43, pp. 61-69. [DOI: https://dx.doi.org/10.2307/2531949]