[go: up one dir, main page]
More Web Proxy on the site http://driver.im/
Next Article in Journal
Forecasting GRACE Data over the African Watersheds Using Artificial Neural Networks
Previous Article in Journal
Accuracy of the TanDEM-X Digital Elevation Model for Coastal Geomorphological Studies in Patagonia (South Argentina)
You seem to have javascript disabled. Please note that many of the page functionalities won't work as expected without javascript enabled.
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Remote Sensing
Volume 11
Issue 15
10.3390/rs11151768
remotesensing-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Winter Wheat Yield Assessment from Landsat 8 and Sentinel-2 Data: Incorporating Surface Reflectance, Through Phenological Fitting, into Regression Yield Models

by
Sergii Skakun
1,2,3,*,
Eric Vermote
3,
Belen Franch
1,3,
Jean-Claude Roger
1,3,
Nataliia Kussul
4,
Junchang Ju
5,6 and
Jeffrey Masek
6
1
Department of Geographical Sciences, University of Maryland, College Park, MD 20742, USA
2
College of Information Studies (iSchool), University of Maryland, College Park, MD 20742, USA
3
NASA Goddard Space Flight Center Code 619, 8800 Greenbelt Road, Greenbelt, MD 20771, USA
4
Space Research Institute NAS Ukraine & SSA Ukraine, 03680 Kyiv, Ukraine
5
Earth System Science Interdisciplinary Center (ESSIC), University of Maryland, College Park, MD 20742, USA
6
NASA Goddard Space Flight Center Code 618, 8800 Greenbelt Road, Greenbelt, MD 20771, USA
*
Author to whom correspondence should be addressed.
Remote Sens. 2019, 11(15), 1768; https://doi.org/10.3390/rs11151768
Submission received: 9 July 2019 / Revised: 24 July 2019 / Accepted: 25 July 2019 / Published: 27 July 2019
(This article belongs to the Section Remote Sensing in Agriculture and Vegetation)
Browse Figures
Graphical abstract
"> Figure 1
<p>Location of Kirovohradska oblast in Ukraine and its division into districts. Shown also is the location of Sentinel-2 tiles covering the area. The background image shows a winter crop mask derived from Landsat 8 and Sentinel-2 for 2018.</p> "> Figure 2
<p>General workflow for winter wheat mapping and yield assessment.</p> "> Figure 3
<p>(<b>a</b>) Average number of cloud-free observations for winter crop pixels depending on satellite data usage. The number of pixels was taken from March until the end of June, which was the period of winter crop growth. (<b>b</b>) Spatial distribution of the number of cloud-free observations from Landsat 8 and Sentinel-2 from March through the end of June in 2018.</p> "> Figure 4
<p>Peak and accumulation approaches for winter wheat yield assessment. Shown is the difference vegetation index (DVI), derived from Landsat 8 and Sentinel-2 data over a winter crop pixel, and fitted against accumulated growing degree days (AGDD) using a quadratic relationship.</p> "> Figure 5
<p>DVI (<b>a</b>) and NIR-red (<b>b</b>) dynamics for the same district over three years with variations in winter wheat yields.</p> "> Figure 6
<p>Comparison between winter crop areas from official statistics and satellite at the district scale in Kirovohradska oblast in 2016 to 2018.</p> "> Figure 7
<p>Regression crop yield models based on DVI and AUC for 2016 (<b>a</b>), 2017 (<b>b</b>), and 2018 (<b>c</b>).</p> "> Figure 8
<p>Performance of the AUC-DVI empirical winter wheat yield model on calibration data depending on the combined use of Landsat 8/OLI and Sentinel-2/MSI and a single usage. The coefficient of determination <span class="html-italic">R</span><sup>2</sup> (<b>a</b>) and relative RMSE (<b>b</b>) are shown.</p> "> Figure 9
<p>Estimates vs. reference winter wheat yields within the temporal cross-validation for the best AUC, coefficients–NIR + red model.</p> ">
Versions Notes

Abstract

:
A combination of Landsat 8 and Sentinel-2 offers a high frequency of observations (3–5 days) at moderate spatial resolution (10–30 m), which is essential for crop yield studies. Existing methods traditionally apply vegetation indices (VIs) that incorporate surface reflectances (SRs) in two or more spectral bands into a single variable, and rarely address the incorporation of SRs into empirical regression models of crop yield. In this work, we address these issues by normalizing satellite data (both VIs and SRs) derived from NASA’s Harmonized Landsat Sentinel-2 (HLS) product, through a phenological fitting. We apply a quadratic function to fit VIs or SRs against accumulated growing degree days (AGDDs), which affects the rate of crop development. The derived phenological metrics for VIs and SRs, namely peak, area under curve (AUC), and fitting coefficients from a quadratic function, were used to build empirical regression winter wheat models at a regional scale in Ukraine for three years, 2016–2018. The best results were achieved for the model with near infrared (NIR) and red spectral bands and derived AUC, constant, linear, and quadratic coefficients of the quadratic model. The best model yielded a root mean square error (RMSE) of 0.201 t/ha (5.4%) and coefficient of determination R2 = 0.73 on cross-validation.

Graphical Abstract">

Graphical Abstract

1. Introduction

The combination of data acquired by Landsat 8 and Sentinel-2 remote sensing satellites can provide a high temporal resolution (3–5 days) [1], which is critical for various applications requiring dense satellite data time series. Previously, such (or better) high temporal resolution was available mainly for Earth observation sensors, which acquire daily data over Earth’s surface, but at a coarser spatial resolution (>250 m) [2,3]. The latter, for example, includes space-borne remote sensing sensors, such as Moderate Resolution Imaging Spectroradiometer (MODIS), Visible Infrared Imaging Radiometer Suite (VIIRS), Advanced Very High Resolution Radiometer (AVHRR), and SPOT-VEGETATION. Taking into account an increased frequency of observations at moderate spatial resolution (<30 m), the assumption is that methods for generating products from coarse spatial resolution sensors can be ported to moderate (Landsat 8/OLI, Sentinel-2/MSI) or high (Planet/PlanetScope) spatial resolution sensors. However, the practice shows that such a transition is not always straightforward due to larger data gaps because of clouds and uneven coverage, sensor characteristics, and increased spatial resolution (at least at the order of 100, when going from 250 to 30 m).
Consider, for example, a crop yield assessment/forecasting application [4]. The hypothesis is that satellite-based features, such as vegetation indices (VIs) or biophysical parameters derived at a single date or accumulated over some time period, can be correlated to crop yields [5,6]. Since the reference data on crop yields are mainly available at the regional scale, the corresponding empirical models are built by averaging satellite-based features over those regions and correlating these derived variables with crop yields [7,8,9]. It is assumed that there is a homogeneity within the region in terms of crops grown and agricultural practices applied and, therefore, the averaging should be performed for satellite data acquired at the same (or approximately the same) stages of crop growth, meaning that the data are normalized. This is usually the case for coarse spatial resolution remote sensing sensors, which enable a higher likelihood of obtaining a high temporal frequency of cloud-free data over the Earth’s surface [10,11]. This is also evidenced by multiple successful applications of coarse spatial resolution satellite data to crop yield assessment and forecasting [5,7,8,9,12,13,14,15].
However, this is not the case for moderate spatial resolution satellite data (<30 m). Irregular spatial coverage, when the area in question is covered by several “stripes” sensed at different times, and lower revisit cycles lead to discrepancies in dates of cloud-free observations. For example, Azzari et al. (2017) [16] encountered this problem, when they used Landsat data for mapping the yields of corn and soybean in the USA. The direct application of their model to available Landsat observations resulted in obvious artefacts, which were due not to spatial variability of crop yields but rather to irregularities in satellite observations at a moderate spatial resolution. To address this problem, they used MODIS observations to further rectify (normalize) Landsat data, so to avoid the effect of irregular Landsat observations on the resulting product. MODIS data are frequently used to assist with the processing of Landsat-like sensors to generate near-daily observations at a 30 m resolution [17]. These studies showed promising results and have already been applied to assessing yields at the field scale [18,19]. However, the effect of combining various spatial resolutions, i.e., >250 m MODIS data with 30 m Landsat data, usually is not quantified in terms of crop yield model performance. For example, Gao et al. (2019) [20] indicated that the error of Landsat-MODIS fusion can be up to 0.083 for the normalized difference vegetation index (NDVI), 0.026 for the red band, and 0.052 for the near infrared (NIR) band. Moreover, previous studies showed saturation of the commonly applied NDVI at a 30 m spatial resolution for winter wheat at high yields (>4 t/ha) or for leaf-area index (LAI) > 2 [21,22], something that was not usually observed for MODIS due to mixed pixels. This means that the behavior of commonly applied VIs at coarse and moderate spatial resolutions might be different.
Recent studies took advantage of Landsat and Sentinel-2 data to address crop yield assessment at a moderate spatial resolution. For example, Lambert et al. (2018) [23] used Sentinel-2 data and a peak LAI approach to estimate the yields of cotton, maize, millet, and sorghum in Mali. The R2 ranged from 0.48 to 0.80 for different crops on training data. However, they did not address the issues of proving that they are actually capturing the peak. Also, this study was based on a single year, therefore the model’s robustness in time was not explored. Lai et al. (2018) [24] used time-integrated Landsat NDVI for wheat yield prediction in Australia. They used an asymmetric bell-shaped growth model to fit NDVI against time. Their approach yielded a root mean square error (RMSE) of 0.79 t/ha on cross-validation. However, they did not go beyond NDVI and did not explore surface reflectance values.
A frequently used approach to smooth satellite data is through phenological fitting. Multiple fitting functions have been proposed, tested, and validated [24,25,26,27,28,29,30,31,32]. The current general consensus is that no single model would fit all requirements and selection of the smoothing method and/or model is predicated by the application, and prior knowledge of the land–surface dynamics [25]. Also, almost all existing remote sensing-based studies for crop yields are focused on using vegetation indices and do not exploit spectral information.
Though a combination of Landsat 8 and Sentinel-2 offers a high frequency of observations, discrepancies in available cloud-free data, however, will still exist. With more and more applications transitioning to a higher spatial resolution, it is important to explore and analyze the effect of those discrepancies on the quality of the resulting products. In the case of crop yield assessment, one of the important questions is as follows: How do discrepancies in the dates of satellite observations influence the performance of empirical crop yield models at the regional scale? The present paper aims to address this question by documenting the results of a study building empirical models for winter wheat yield assessment with Landsat 8 and Sentinel-2 data in Ukraine. Not performing satellite data smoothing may lead to poorer performance of the empirical crop yield model, which would be attributed not to the lack of correlation with satellite-derived variables but rather to observation irregularities. We will not limit satellite data fitting to vegetation indices only, and will also explore surface reflectance values in multiple spectral bands for the assessment of winter wheat yields. Finally, we will evaluate the use of multi-source data, namely Landsat 8 and Sentinel-2, for crop yield assessment. In particular, we assess quantitively how the performance of an empirical model of crop yield depends on the combined use of Landsat 8 and Sentinel-2 data, compared to a single-sensor usage.

2. Materials

2.1. Study Area and Reference Data

The study was performed for Kirohohradska oblast in Ukraine for 2016 to 2018. An “oblast” is a high-level administrative division of the country, and each oblast is further divided into districts. There are 24 oblasts in Ukraine and Autonomous Republic of Crimea. Kirovhradska oblast is located in the central part of Ukraine and consists of 21 districts with a geographical area ranging from 65 to 165 thousand ha and cropland area ranging from 27 to 112 thousand ha (Figure 1).
The reason for selecting this region is that it is a top 10 wheat producer in Ukraine and because of the availability of reference crop yield and harvested area data at the district scale for 2016 to 2018. Winter wheat is one of the major crops in Kirovhradska oblast, accounting for 20% of the production of all crops in the region. Winter wheat is mainly rain-fed in the region and usually planted in September to October. After dormancy during the winter, it re-starts its growth in early spring, reaching maturity by the end of June. The harvest of winter crops is typically undertaken in July.
Reference data on the crop yield and harvested area at the district level were collected from the Department of Agro-Industry Development of Kirovohrad State Administration (http://apk.kr-admin.gov.ua). The data were made available online as the harvest progressed and were based on farm surveys of all large agricultural enterprises (that account for more than 90% of all winter crops produced in the region) and samples of household farms, the same way that official statistics are collected [33]. The final estimates for winter crop yields and harvested areas were available at the end of November and were used as a reference in this study. The reference data went through a quality-control procedure, which resulted in two districts in 2017 and five districts in 2018, being excluded from the analysis, because the reported production was >+/−10% than the area multiped by yields. The average winter wheat yield among the districts was 4.0 ± 0.3 t/ha in 2016, 3.5 ± 0.4 t/ha in 2017, and 3.6 ± 0.4 t/ha in 2018.

2.2. Landsat-8/OLI and Sentinel-2A /MSI Datasets

Remote sensing images acquired by the Operational Land Imager (OLI) instrument aboard the Landsat-8 satellite and by the Multi-Spectral Instrument (MSI) aboard the Sentinel-2A/B satellites were used in the study. Landsat 8/OLI captures images of the Earth’s surface in nine spectral bands at a 30 m spatial resolution (15 m for panchromatic band) [34], while Sentinel-2A/B/MSI captures images of the Earth’s surface in 13 spectral bands at a 10, 20, and 60 m spatial resolution [35]. We used NASA’s Harmonized Landsat Sentinel-2 (HLS) product [36], which provides a Level-2 Nadir BRDF (Bidirectional Reflectance Distribution Function)-Adjusted surface Reflectance (NBAR) at a 30 m spatial resolution. HLS combines, in a single data set, observations of the land surface from the Landsat 8’s OLI [37] and Sentinel-2’s MSI [38]. MSI data are also spectrally adjusted to OLI spectral bands, so to generate a single harmonized dataset. We used HLS version 1.4. A detailed description of HLS product generation is described in [37].
The HLS product is organized using the Sentinel-2 110 km × 110 km tiles. The following tiles cover the study region (Figure 1): 35UQQ, 35UQP, 36UUV, 36UUU, 36UVV, 36UVU, 36UWV, and 36UWU. Overall, 3565 scenes acquired by Landsat 8 and Sentinel-2A/B and covering eight tiles were used in the study. The scene covered the period of March to July within 2016 to 2018 to capture the main development stages of wheat. It should be noted that the availability of Sentinel-2 data varied through 2016 to 2018. In 2016 and 2017, only Sentinel-2A data were available, and only in 2017 did Sentinel-2A start operational acquisitions on a 10-day revisit cycle. In 2018, both Sentinel-2A/B were available with a combined five-day revisit cycle.
HLS data were used to extract surface reflectances (SRs) and to calculate corresponding vegetation indices (VIs). We used spectral bands common to Landsat 8 and Sentinel-2A/B, namely: Green (wavelength ~0.560 µm), red (~0.660 µm), near infrared (NIR) (~0.865 µm, we used band 8A from Sentinel-2), and two short-wave infrared (SWIR) (~1.5 and ~2.2 µm). This is an advantage of using analysis ready data (ARD), such as HLS, since data have been harmonized and as if they are acquired by a single sensor. We used a quality assessment (QA) layer provided in HLS to disregard unreliable pixels identified either as cloud, adjacent to cloud, shadow, snow/ice, or water.
Red and NIR bands were used to calculate the following commonly used VIs:
  • Normalized difference vegetation index (NDVI) [39]:
    N D V I = ρ N I R − ρ r e d ρ N I R + ρ r e d .
  • Enhanced vegetation index (EVI) [40]:
    E V I 2 = 2.5 ρ N I R − ρ r e d ρ N I R + 2.4 ρ r e d + 1 .
  • Difference vegetation index (DVI) [41]:
    D V I = ρ N I R − ρ r e d ,
    where ρ N I R and ρ r e d are SR values in NIR and red spectral bands, respectively.

2.3. Meteorological Data

We used air temperature derived from NASA’s Modern-Era Retrospective analysis for Research and Applications (MERRA2) reanalysis data set [42] to compute the growing degree days (GDDs) for winter wheat. Santamaria-Artigas et al. (2019) [43] showed that MERRA2 is one of the best data sources for computing GDD for winter wheat in terms of error (<2K compared to in situ measurements), spatial resolution, and low delay in data availability. The data are provided on a regular grid that has 576 points in the longitudinal direction and 361 points in the latitudinal direction, corresponding to a resolution of 0.625° × 0.5°. We used the time-averaged two-dimensional data collection (M2SDNXSLV) to extract daily averaged two-meter air temperature (T2MMEAN). The data sets were extracted from the netCDF format, transformed to the geo-referenced raster, subset for study areas, and linearly interpolated to the 30-m spatial resolution.

3. Methods

3.1. General Overview

Figure 2 shows main components and dataflows of the study with generation of winter crop masks and winter wheat yield estimation.
The crop-specific masks were generated on a yearly basis, so the signal for building empirical yield models was extracted for winter crops only. Satellite data, namely SRs and VIs, were first smoothed before incorporating them into the crop yield model. This was performed through phenological fitting against meteorological data, namely accumulated GDD (AGDD), since temperature is a primary factor affecting the rate of crop development [44]. The following subsections describe in detail the methodologies used.

3.2. Winter Crop Type Mapping

In order to provide in-season winter crop maps (2016–2018), we implemented an automatic approach, which was initially developed for the coarse spatial resolution MODIS sensor [45] and further extended and prototyped for Landsat 8/OLI and Sentinel-2/MSI at a 30 m spatial resolution [22]. The approach uses a phenological metric during the green-up stage of winter crop development, namely the maximum NDVI, which is used as a main feature to discriminate winter crops from spring and summer crops. A cropland mask [46] is an essential data input for the method, to constrain the analysis to cultivated areas only, and to reduce potential confusion with natural vegetation and forested areas. To perform discrimination in an automatic way, a Gaussian mixture model (GMM) [47] is applied with the number of GMM components set to 2. Setting this number is based on the assumption that there will be two modes in the distribution of maximum NDVI (phenological metric), which would correspond to winter and non-winter crops, respectively. This method was run in an automatic unsupervised way during the vegetation season for 2016 to 2018 without the need to provide ground truth. However, this method requires ancillary data inputs and a priori knowledge about the crop calendar. We refer readers to [22,45] for details.

3.3. Winter Wheat Yield Assessment

Crop yield assessment was done through the development of an empirical regression model between satellite-derived features (regressors) and reference crop yields at the district scale (Figure 1). Satellite-derived features were averaged over the crop mask (Section 3.2) for each district and correlated with the available reference winter wheat yields (Section 2.1). Before averaging satellite-based features over pixels identified as winter crops, one has to ensure that satellite data are phenologically normalized, so averaging is performed for the same stages of crop growth development. Though a combination of Landsat 8 and Sentinel-2 offers a high frequency of observations, discrepancies in available cloud-free data, however, will still exist (Figure 3).
We used a quadratic model between satellite data (VI or SR) and AGDD [26,48,49]:
ysat = a0 + a1*AGDD + a2*AGDD2,
where ysat is the VI (NDVI, EVI2, or DVI) or SR (green, red, NIR, SWIR1, or SWIR2); and a0, a1, and a2 are the constant, linear, and quadratic coefficients, respectively. GDD was defined as the average daily maximum (Tmax) and minimum temperatures (Tmin) minus a base temperature (Tbase):
G D D = T m a x + T m i n 2 − T b a s e ,
where if [(Tmax + Tmin)/2 < Tbase], then GDD = 0 [50], and with Tbase = 0. AGDD was calculated by summing GDDs (Equation (5)).
Having derived per-pixel relationships between VI/SR and AGDD (Equation (4)), it is possible to derive various phenological metrics to correlate with yields or incorporate coefficients (a0, a1, a2) directly into the yield model. The generalized relationship between crop yields and satellite-derived features can be expressed as follows:
y y i e l d l = ω 0 + ∑ j = 1 n r ω j r j l ,
where y y i e l d l is the crop yield for district l; n r is the number of regressors used; ω j | j = 0 n r are regression coefficients, and r j l are satellite-derived features averaged for district l over the crop mask:
r j l = 1 | A l | ∑ p ∈ A l r j , p l ,
where A l is the set of winter crop pixels; r j , p l is the satellite-derived feature for the winter crop pixel at location p; and | · | is the cardinality number, i.e., the number of elements in the set.
Multiple approaches were considered in the study that varied in the regressors used (VI, SR, or combination of SRs) and the number of regressors, as well as phenological metrics. For VIs, two phenological metrics were derived from the quadratics relationship: Peak and area under curve (AUC) (Figure 4).
The peak approach assumes that yield is correlated to the peak biomass during a crop growth season [7,8,9,23,37,51,52], which would result in the following regression (we removed l notations for simplicity):
y y i e l d = ω 0 + ω 1 r 1 ,   where   r 1   is   p V I = a 0 − a 1 2 4 a 2 ,
where p V I is the peak of the VI during the vegetation season.
The accumulation, or integrated, approach assumes that yield is proportional to the accumulated biomass over the crop growth season [18,53,54,55,56]:
y y i e l d = ω 0 + ω 1 r 1 ,   r 1 ≡ s V I = ∫ A G D D m i n A G D D p e a k y s a t , V I d x = ∫ A G D D m i n A G D D p e a k ( a 0 + a 1 x + a 2 x 2 ) d x ,
where A G D D m i n and A G D D p e a k define ranges, for which accumulation is performed. A G D D m i n was determined based on the minimum VI values, namely 0.05 for DVI and 0.15 for NDVI and EVI2. These thresholds for VIs were determined empirically to make sure that accumulation was performed for vegetated stages and not bare land. A G D D p e a k is the AGDD value, for which the VI’s peak is achieved and is equal to − a 1 2 a 2 (Equation (4)).
The p-value was also calculated to test the null hypothesis that a coefficient for this regressor is equal to zero, i.e., the regressor has no effect on the dependent variable.
It should be noted that almost all existing approaches dealing with satellite data and crop yields utilize VIs (e.g., Equations (1)–(3)), which incorporate two or more spectral bands into a single variable. With fitted SRs, one can explore trajectories in a multi-dimensional space of SRs for different spectral bands. Figure 5 shows examples of 1-D and 2-D trajectories for DVI and NIR-red, respectively, for different winter wheat yield values.
Therefore, the following question arises: What is the added value that can be achieved, when incorporating surface reflectances into empirical regression crop yield models? However, adding more regressors in the model (Equation (6)), while having a limited number of reference data, might lead to model overfitting—the model’s performance will improve on training (calibration) data, while failing to generalize, i.e., to provide a reasonable performance on independent (previously unseen) data. In order to prevent overfitting of the model (Equation (6)) with multiple regressors, we used a ridge regression [57] with the L2 (or Tikhonov’s) regularization [58,59]:
W ^ = argmin W ( ‖ Y − W R ‖ 2 + α   ‖ W 2 ‖ ) ,
where Y, R, and W are vectors of yields, regressors, and weights of the regression, respectively, and α is the regularization coefficient. The regularization term, α   ‖ W 2 ‖ , prevents coefficients of the regression increasing during the minimization process, and therefore prevents overfitting, and might improve the generalization ability of the model. Such a constraint can either force the coefficients, W, to be “sparse” or reflect prior knowledge, namely correlation between regressors. Table 1 summarizes the regression models used in the study. The optimal regularization parameter, α, is determined during the training (calibration) process by dividing calibration data into training and testing datasets, where α is selected in such a way that it maximizes a model’s performance for testing data. Note that validation data were not used for selecting the optimal regularization parameter.

3.4. Implementation and Performance Evaluation

VIs (NDVI, EVI2, and DVI—Equations (1)–(3)) and SRs (green, red, NIR, SWIR1, SWIR2) were derived from HLS over the study area. Regression models of yields as a function of satellite-derived features were evaluated using a standard set of metrics, such as coefficient of determination (R2), root mean square error (RMSE) between estimate and reference, and relative RMSE (RRMSE):
R 2 = 1 − ∑ i = 1 n ( y e s t − y r e f ) 2 ∑ i = 1 n ( y ¯ r e f − y r e f ) 2 ,
R M S E = 1 n ∑ i = 1 n ( y e s t − y r e f ) 2 ,
R R M S E = 100 % × R M S E y ¯ r e f ,
where y e s t and y r e f are estimate and reference values, respectively; y ¯ r e f = 1 n ∑ i = 1 n y r e f is the average of the reference values, and n is the number of samples.
The developed crop yield models were evaluated with two versions of leave-one-out cross-validations (CVs). Within spatial CV, we iteratively removed all samples for a particular district, which was reserved for testing purposes, and the model was calibrated on the rest of the districts using all three years (2016–2018). Within temporal CV, the model was trained on two years, while the third year was reserved for testing. Spatial CV allowed us to explore the spatial robustness of the empirical models, while temporal CV allowed us to explore the temporal robustness.

4. Results

4.1. Winter Crop Type Mapping

Figure 6 shows results of the validation for 2016 to 2018, when comparing the areas of winter crops from official statistics and those derived from satellite data, namely Landsat 8 and Sentinel-2, using the GMM approach.
Overall, there is a good correspondence between satellite-derived and reference areas, with an average R2 of 0.91 and a bias of +4.3 thousand ha. This bias might be attributed to the confusion between winter crops and early spring cereal crops [60].

4.2. Fitting VI and SR

Table 2 provides a summary of fitting metrics for various VI and SR values derived from HLS for 2016 to 2018.
All VIs (DVI, EVI2, and NDVI) showed a very good correspondence to a quadratic model (Equation (4)), with an average R2 > 0.9. Among SRs, the best performance was for NIR, with an average R2 = 0.88. Other SRs, except green, showed a variable performance, with R2 ranging from 0.7 to 0.9, with the green spectral band yielding an average R2 of 0.6. Nevertheless, all VIs and SRs (except green) exhibited good performance with a relative error of fitting less than 10% (for the green band, the error was 10.2%). It should also be noted that the fitting results were stable across multiple years.

4.3. Comparison Between VI-Based Yield Models

A comparison of VI-based empirical models for different years is presented in Table 3.
The biggest impact of smoothing satellite data for peak-VI models was observed in 2016, where fewest cloud-free data were available compared to 2017 and 2018 (Figure 3a). The coefficient of determination, R2, was 0.179, 0.056, and 0.088 for the peak-DVI, EVI2, and NDVI models, respectively, without phenological fitting, when the peak was calculated directly from available HLS datasets; those R2 values increased to 0.332, 0.282, and 0.485, when calculating the peak p V I (Equation (8)) from phenological fitting. The best performance for 2016 was, however, for the AUC-DVI approach, with an R2 = 0.588. This approach gave the best performance among all VI-based approaches considered, providing a consistent performance of R2 = 0.589 and R2 = 0.608 for 2017 and 2018, respectively (Figure 7).
Performance of EVI2 and NDVI varied within 2016 to 2018, yielding a good performance for one year and failing to produce good results for another year. Interestingly, all models consistently improved from 2016 to 2018 thanks, most probably, to the availability of more cloud-free data (Figure 3a). On average, the number of cloud-free observations for winter crops within the March to June period for 2018 was 1.3 and 1.8 times larger compared to 2017 and 2016, respectively. This shows the importance of having dense and high-frequency satellite observations at a moderate spatial resolution (10–30 m) for crop yield studies.

4.4. Impact of the Combined Use of Landsat 8 and Sentinel-2 Data

We used the AUC-DVI crop yield model as a benchmark to evaluate the combined use of Landsat 8/OLI and Sentinel-2/MSI (both A and B satellites) data as compared to a single sensor usage. Performance metrics, namely R2 and RRMSE, are shown in Figure 8. As with the results from the previous subsection, the biggest impact was observed in 2016. With a single sensor usage (OLI or MSI), it was impossible to derive satellite-based features for some districts in 2016, because of not having enough cloud-free observations. In 2016 (March–June period), there were on average 2.6 ± 1.4 cloud-free observations from Landsat 8/OLI and 3.6 ± 1.6 cloud-free observations from Sentinel-2/MSI. With a single sensor usage in 2016, there was a poor correlation observed between crop yields and satellite-derived features, even with phenological fitting.
In 2017 and 2018, the combined use of Landsat 8/OLI and Sentinel-2/MSI outperformed a single sensor usage, though less impact was observed in 2018, when both Sentinel-2A and -2B provided a nominal five-day revisit cycle. Though Landsat 8/OLI has a 16-day revisit cycle (8 days for overlapping scenes) compared to 5 days from both Sentinel-2’s, it provided a valuable contribution, accounting for, on average, 30% to 40% of total cloud-free observations in the March to June period in the study area. These results highlight the importance of multi-source satellite data usage for crop yield studies.

4.5. Comparison Between VI- and SR-Based Yield Models

Results of the cross-validation (CV) of different models incorporating VIs, SRs, and derived AUC (Equation (9)) and coefficients a0, a1, and a2 of the parabola (Equation (4)) are presented in Table 4.
Among the single regressor-based models, the best performance was achieved using an NIR spectral band, with an RMSE = 0.236 t/ha (6.3%) and R2 = 0.63 for temporal CV. This model outperformed VI- and SR-based models. Among VIs, the best performance was achieved using DVI, with an RMSE = 0.257 t/ha (6.9%) and R2 = 0.6.
Having two AUC-SR-based regressors did not improve on the model compared to the AUC–NIR based model. The best results were achieved for the AUC–NIR + green model (RMSE = 0.244 t/ha, 6.5%, R2 = 0.62), which did not improve on the single regressor AUC–NIR; however, it did improve the AUC–VI-based models. This shows the importance of adding separate SR values, rather than using VIs with pre-calculated weights for spectral bands. Surprisingly, the AUC–NIR outperformed the AUC–NIR + red model and other combinations of spectral bands, meaning that AUC-derived features for green, red, and SWIR1 spectral bands did not add much information in terms of winter wheat yields. SWIR2 data were not included in the analysis, because of considerable correlation between SWIR1 and SWIR2.
The addition to the AUC of the coefficients of the quadratic relationship (Equation (4)) as regressors, which define the full dynamics in the multi-dimensional spectral space (Figure 5), improved estimations of winter wheat yields. In particular, the AUC, coefficients–NIR + red model showed the best overall results, with RMSE = 0.201 t/ha (5.4%) and R2 = 0.73 (Figure 9). It shows that not only is the AUC important but also the full trajectories itself, defined by coefficients a0, a1, a2. Since this model has eight regressors (AUC, a0, a1, a2 for NIR and red), it easily overfits for a small sample size, such as in this study, and, therefore, the role of regularization (Equation (10)) becomes extremely important. The derived optimal regularization coefficient was α = 10−5. Without regularization, the AUC, coefficients–NIR + red model yields an RMSE = 0.296 t/ha (7.9%) and R2 = 0.49 performance, which is at least 30% worse than with regularization. Adding all spectral bands (NIR, red, green, SWIR1) to the AUC, coefficients-based model did not improve the AUC, coefficients–NIR + red model; however, it showed the second best result among all the models considered in the study along with AUC, coefficients–DVI model. It is also worth noting that performance metrics within spatial CV were consistently better than metrics for temporal CV, meaning a better robustness of the models in the spatial context, rather than temporal.

5. Discussion

Remote sensing has always been a valuable source of information for crop yield. Previous success has been achieved chiefly by using coarse spatial resolution satellite data, thanks to its high temporal resolution (daily), which is needed to capture the behavior of crop growth. Vegetation indices that incorporate two or more spectral bands with pre-defined weights have usually been applied to correlate with crop yields. With recent advancements in technologies, high temporal resolution satellite data have been available at moderate (10–30 m, combination of free Landsat 8 and Sentinel-2) and high (~3 m, commercial Planet) spatial resolutions. However, most existing studies continue to exploit vegetation indices, disregarding mainly surface reflectance values (SRs) in the available spectral bands. This study aimed to address this question and explored the use of SRs, namely green (~0.560 µm), red (~0.660 µm), NIR (~0.865 µm), and SWIR1 (~1.5 µm), when estimating winter wheat yields. In order to account for discrepancies in satellite data acquisitions, we applied a quadratic relationship (Equation (4)) to SRs of individual spectral bands against accumulated GDD (Equation (5)). This relationship was previously applied to VIs only and has never been applied to SRs. Our results showed that this quadratic relationship can be effectively applied to SRs as well, with a relative error of fitting ranging from 5% to 12%, and R2 from 0.59 to 0.91. The best results of phenological fitting for SR values were observed for NIR (average R2 = 0.88), while the green spectral band showed the worst performance of R2 = 0.6. While these results show the adequacy of a quadratic function, we would like to emphasize the potential limitations of this fitting model. The use of a parabola implies that the upward phase of crop growth is mirrored in the senescence, which is the case for our study area. However, in the regions, where water is more limiting to crop growth than sunlight, such an assumption might be incorrect. Therefore, when fitting the parabola to the satellite-derived VIs or SRs, one should always check the fitting performance, and in case of poor fitting, other fitting models should be explored to better capture crop phenology.
By having those quadratic relationships in place for SRs, multiple phenological metrics can be derived. A similar approach with multiple phenological metrics has been successfully applied to cropland and grassland mapping in the USA for multiple years using Landsat and MODIS data, though only EVI was used for fitting [50]. Therefore, phenological metrics, such as area under curve (AUC) (Figure 4) and coefficients of the quadratic relationship (Equation (4)), can be used to estimate crop yields. Since adding multiple regressors into a regression model might lead to model overfitting (especially in the case of a small sample size), we proposed the use of ridge regression [53], which incorporates an L2 regularization [58,59]. Our best winter wheat yield model was a combination of AUC and constant, linear, and quadratic coefficients of the quadratic model for NIR and red spectral bands, with an RMSE = 0.201 t/ha (5.4%) and R2 = 0.73 (Figure 9). This model outperformed VI-based models. Adding other spectral bands, namely green and SWIR1, did not improve the model performance for winter wheat yields. These results confirm previous findings for winter wheat yields studies, where NIR and red play the most important role in capturing the biophysical parameters of winter wheat [12,61]. For other crops, for example, corn and soybean, a similar approach can be applied to identify spectral bands that maximize the correlation with crop yields. However, these should be investigated further.
We also investigated— in a numerical way—the efficiency of the combined use of Landsat 8/OLI and Sentinel-2/MSI for winter wheat yield estimation. All three years featured in the study (2016–2018) had different scenarios of satellite data acquisitions, which benefited such an analysis. The most substantial impact was observed in 2016, when only one of the two Sentinel-2 satellites was available, however, was not fully operational (revisit cycle >10 days). Using a single sensor only in 2016, either Landsat 8/OLI or Sentinel-2A/MSI, did not give any reasonable results for the crop yield model, since each satellite contributed with 2.6 ± 1.4 and 3.6 ± 1.6 cloud-free observations in the March to June period, respectively, which was not enough to capture the growth dynamics of winter wheat. The single-sensor-only models provided R2 = 0.01 and R2 = 0.34 for Landsat 8 and Sentinel-2A, respectively. A combination of satellites increased the number of observations to 6.2 ± 2.2, and the performance of the winter wheat yield model improved to R2 = 0.59. For 2017 to 2018, when data from Sentinel-2B started to be available and both Sentinel-2A/B became operational with a combined five-day revisit cycle, that impact was not so dramatic; nevertheless, the combination of Landsat 8 and Sentinel-2 improved the performance of the winter wheat yield model by 18% and 15% in 2017 and 2018, respectively, compared to Sentinel-2A/B only. Though Landsat 8 has a nominal 16-day revisit cycle (8-day for overlapping scenes), which is by far less than two satellites Sentinel-2A/B, it consistently added on average 3.2 ± 1.6 cloud-free observations (during the March to June period), which constituted on average 39% of the total observations. These results further highlight the importance of multi-source data usage in these types of studies.

6. Conclusions

In this paper, we focused on the estimation of winter wheat yield at a regional scale in Ukraine, using the combined Landsat 8 and Sentinel-2 data at a 30-m spatial resolution. We emphasized the importance of normalizing (smoothing) satellite data in crop yield studies through, for example, phenological fitting. Not doing such normalization may result in a considerable decrease in crop yield model performance. Unlike most existing approaches utilizing vegetation indices only, we derived per-pixel phenological metrics directly from all spectral bands (green, red, NIR, and SWIR), when fitting them against AGDD. The derived phenological metrics were correlated against crop yields using a ridge regression (with L2 regularization) in order to avoid model overfitting. The best results were achieved for the model with NIR and red spectral bands and derived area under curve (AUC), constant, linear, and quadratic coefficients of the quadratic model. The best model yielded RMSE = 0.201 t/ha (5.4%) and R2 = 0.73 on cross-validation.
We also analyzed the importance of multi-source data usage for crop yield estimation. For the area considered in this study, both Landsat 8/OLI and Sentinel-2A/B/MSI data were valuable, and the combined use yielded better performance than a single sensor usage for all three years (2016–2018). These results reiterate the need for high-frequency observations at moderate and higher spatial resolutions for crop yield studies.
Further works should be directed towards an investigation of other crops, such as corn and soybeans, as well as identification of efficient ways of incorporating synthetic aperture radar (SAR) observations to improve the frequency of data in areas with high cloud cover.

Author Contributions

Conceptualization and methodology, S.S., E.V. and B.F.; data curation, S.S., J.J. and J.M.; validation, S.S. and N.K.; writing—original draft preparation, S.S.; writing—review and editing, E.V., B.F., J.-C.R. and N.K.; supervision, E.V. and J.-C.R.; funding acquisition, S.S., B.F. and J.-C.R.

Funding

This research was funded by the NASA Land-Cover/Land-Use Change (LCLUC) Program, grant number 80NSSC18K0336, and Cooperative Agreement NASA Harvest, grant number 80NSSC18M0039.

Acknowledgments

We would like to thank NASA for providing support in generating the Harmonized Landsat Sentinel-2 (HLS) product. This study used modified Copernicus Sentinel data (2016–2018).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Li, J.; Roy, D.P. A Global Analysis of Sentinel-2A, Sentinel-2B and Landsat-8 Data Revisit Intervals and Implications for Terrestrial Monitoring. Remote Sens. 2017, 9, 902. [Google Scholar] [Green Version]
  2. Justice, C.O.; Vermote, E.; Townshend, J.R.; DeFries, R.; Roy, D.P.; Hall, D.K.; Salomonson, V.V.; Privette, J.L.; Riggs, G.; Strahler, A.; et al. The Moderate Resolution Imaging Spectroradiometer (MODIS): Land remote sensing for global change research. IEEE Trans. Geosci. Remote Sens. 1998, 36, 1228–1249. [Google Scholar] [CrossRef]
  3. Justice, C.O.; Román, M.O.; Csiszar, I.; Vermote, E.F.; Wolfe, R.E.; Hook, S.J.; Friedl, M.; Wang, Z.; Schaaf, C.B.; Miura, T.; et al. Land and cryosphere products from Suomi NPP VIIRS: Overview and status. J. Geophys. Res. Atmos. 2013, 118, 9753–9765. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  4. Becker-Reshef, I.; Justice, C.; Sullivan, M.; Vermote, E.; Tucker, C.; Anyamba, A.; Small, J.; Pak, E.; Masuoka, E.; Schmaltz, J.; et al. Monitoring global croplands with coarse resolution earth observations: The Global Agriculture Monitoring (GLAM) project. Remote Sens. 2010, 2, 1589–1609. [Google Scholar] [CrossRef]
  5. Johnson, D.M. A comprehensive assessment of the correlations between field crop yields and commonly used MODIS products. Int. J. Appl. Earth Obs. Geoinf. 2016, 52, 65–81. [Google Scholar] [CrossRef] [Green Version]
  6. Meroni, M.; Fasbender, D.; Balaghi, R.; Dali, M.; Haffani, M.; Haythem, I.; Hooker, J.; Lahlou, M.; Lopez-Lozano, R.; Mahyou, H.; et al. Evaluating NDVI data continuity between SPOT-VEGETATION and PROBA-V missions for operational yield forecasting in North African countries. IEEE Trans. Geosci. Remote Sens. 2016, 54, 795–804. [Google Scholar] [CrossRef]
  7. Becker-Reshef, I.; Vermote, E.; Lindeman, M.; Justice, C. A generalized regression-based model for forecasting winter wheat yields in Kansas and Ukraine using MODIS data. Remote Sens. Environ. 2010, 114, 1312–1323. [Google Scholar] [CrossRef]
  8. Franch, B.; Vermote, E.F.; Becker-Reshef, I.; Claverie, M.; Huang, J.; Zhang, J.; Justice, C.; Sobrino, J.A. Improving the timeliness of winter wheat production forecast in the United States of America, Ukraine and China using MODIS data and NCAR Growing Degree Day information. Remote Sens. Environ. 2015, 161, 131–148. [Google Scholar] [CrossRef]
  9. 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, 192–203. [Google Scholar] [CrossRef]
  10. Whitcraft, A.K.; Becker-Reshef, I.; Killough, B.D.; Justice, C.O. Meeting earth observation requirements for global agricultural monitoring: An evaluation of the revisit capabilities of current and planned moderate resolution optical earth observing missions. Remote Sens. 2015, 7, 1482–1503. [Google Scholar] [CrossRef]
  11. Whitcraft, A.K.; Vermote, E.F.; Becker-Reshef, I.; Justice, C.O. Cloud cover throughout the agricultural growing season: Impacts on passive optical earth observations. Remote Sens. Environ. 2015, 156, 438–447. [Google Scholar] [CrossRef]
  12. Franch, B.; Vermote, E.F.; Skakun, S.; Roger, J.C.; Becker-Reshef, I.; Murphy, E.; Justice, C. Remote sensing based yield monitoring: Application to winter wheat in United States and Ukraine. Int. J. Appl. Earth Obs. Geoinf. 2019, 76, 112–127. [Google Scholar] [CrossRef]
  13. Kolotii, A.; Kussul, N.; Shelestov, A.; Skakun, S.; Yailymov, B.; Basarab, R.; Lavreniuk, M.; Oliinyk, T.; Ostapenko, V. Comparison of biophysical and satellite predictors for wheat yield forecasting in Ukraine. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 2015, 40, 39–44. [Google Scholar] [CrossRef]
  14. López-Lozano, R.; Duveiller, G.; Seguini, L.; Meroni, M.; García-Condado, S.; Hooker, J.; Leo, O.; Baruth, B. Towards regional grain yield forecasting with 1km-resolution EO biophysical products: Strengths and limitations at pan-European level. Agric. For. Meteorol. 2015, 206, 12–32. [Google Scholar] [CrossRef]
  15. Mkhabela, M.S.; Bullock, P.; Raj, S.; Wang, S.; Yang, Y. Crop yield forecasting on the Canadian Prairies using MODIS NDVI data. Agric. For. Meteorol. 2011, 151, 385–393. [Google Scholar] [CrossRef]
  16. Azzari, G.; Jain, M.; Lobell, D.B. Towards fine resolution global maps of crop yields: Testing multiple methods and satellites in three countries. Remote Sens. Environ. 2017, 202, 129–141. [Google Scholar] [CrossRef]
  17. Gao, F.; Masek, J.; Schwaller, M.; Hall, F. On the blending of the Landsat and MODIS surface reflectance: Predicting daily Landsat surface reflectance. IEEE Trans. Geosci. Remote Sens. 2006, 44, 2207–2218. [Google Scholar]
  18. Gao, F.; Anderson, M.C.; Zhang, X.; Yang, Z.; Alfieri, J.G.; Kustas, W.P.; Mueller, R.; Johnson, D.M.; Prueger, J.H. Toward mapping crop progress at field scales through fusion of Landsat and MODIS imagery. Remote Sens. Environ. 2017, 188, 9–25. [Google Scholar] [CrossRef] [Green Version]
  19. Yang, Y.; Anderson, M.C.; Gao, F.; Wardlow, B.; Hain, C.R.; Otkin, J.A.; Alfieri, J.; Yang, Y.; Sun, L.; Dulaney, W. Field-scale mapping of evaporative stress indicators of crop yield: An application over Mead, NE, USA. Remote Sens. Environ. 2018, 210, 387–402. [Google Scholar] [CrossRef]
  20. Gao, F.; Anderson, M.; Daughtry, C.; Johnson, D. Assessing the Variability of Corn and Soybean Yields in Central Iowa Using High Spatiotemporal Resolution Multi-Satellite Imagery. Remote Sens. 2018, 10, 1489. [Google Scholar] [CrossRef]
  21. Herrmann, I.; Pimstein, A.; Karnieli, A.; Cohen, Y.; Alchanatis, V.; Bonfil, D.J. LAI assessment of wheat and potato crops by VENμS and Sentinel-2 bands. Remote Sens. Environ. 2011, 115, 2141–2151. [Google Scholar] [CrossRef]
  22. Skakun, S.; Vermote, E.; Roger, J.C.; Franch, B. Combined Use of Landsat-8 and Sentinel-2A Images for Winter Crop Mapping and Winter Wheat Yield Assessment at Regional Scale. AIMS Geosci. 2017, 3, 163–186. [Google Scholar] [CrossRef] [PubMed]
  23. Lambert, M.J.; Traoré, P.C.S.; Blaes, X.; Baret, P.; Defourny, P. Estimating smallholder crops production at village level from Sentinel-2 time series in Mali’s cotton belt. Remote Sens. Environ. 2018, 216, 647–657. [Google Scholar] [CrossRef]
  24. Lai, Y.R.; Pringle, M.J.; Kopittke, P.M.; Menzies, N.W.; Orton, T.G.; Dang, Y.P. An empirical model for prediction of wheat yield, using time-integrated Landsat NDVI. Int. J. Appl. Earth Obs. Geoinf. 2018, 72, 99–108. [Google Scholar] [CrossRef]
  25. Atkinson, P.M.; Jeganathan, C.; Dash, J.; Atzberger, C. Inter-comparison of four models for smoothing satellite sensor time-series data to estimate vegetation phenology. Remote Sens. Environ. 2012, 123, 400–417. [Google Scholar] [CrossRef]
  26. De Beurs, K.M.; Henebry, G.M. Land surface phenology, climatic variation, and institutional change: Analyzing agricultural land cover change in Kazakhstan. Remote Sens. Environ. 2004, 89, 497–509. [Google Scholar] [CrossRef]
  27. Fisher, J.I.; Mustard, J.F.; Vadeboncoeur, M.A. Green leaf phenology at Landsat resolution: Scaling from the field to the satellite. Remote Sens. Environ. 2006, 100, 265–279. [Google Scholar] [CrossRef]
  28. Ganguly, S.; Friedl, M.A.; Tan, B.; Zhang, X.; Verma, M. Land surface phenology from MODIS: Characterization of the Collection 5 global land cover dynamics product. Remote Sens. Environ. 2010, 114, 1805–1816. [Google Scholar] [CrossRef] [Green Version]
  29. Jönsson, P.; Cai, Z.; Melaas, E.; Friedl, M.A.; Eklundh, L. A Method for Robust Estimation of Vegetation Seasonality from Landsat and Sentinel-2 Time Series Data. Remote Sens. 2018, 10, 635. [Google Scholar] [CrossRef]
  30. Roy, D.P.; Yan, L. Robust Landsat-based crop time series modelling. Remote Sens. Environ. 2019. [Google Scholar] [CrossRef]
  31. Zhang, X.; Friedl, M.A.; Schaaf, C.B.; Strahler, A.H.; Hodges, J.C.; Gao, F.; Reed, B.C.; Huete, A. Monitoring vegetation phenology using MODIS. Remote Sens. Environ. 2003, 84, 471–475. [Google Scholar] [CrossRef]
  32. Zhang, X.; Friedl, M.A.; Schaaf, C.B. Global vegetation phenology from Moderate Resolution Imaging Spectroradiometer (MODIS): Evaluation of global patterns and comparison with in situ measurements. J. Geophys. Res. Biogeosci. 2006, 111, G04017. [Google Scholar] [CrossRef]
  33. Gallego, F.J.; Kussul, N.; Skakun, S.; Kravchenko, O.; Shelestov, A.; Kussul, O. Efficiency assessment of using satellite data for crop area estimation in Ukraine. Int. J. Appl. Earth Obs. Geoinf. 2014, 29, 22–30. [Google Scholar] [CrossRef]
  34. Roy, D.P.; Wulder, M.A.; Loveland, T.R.; Woodcock, C.E.; Allen, R.G.; Anderson, M.C.; Helder, D.; Irons, J.R.; Johnson, D.M.; Kennedy, R.; et al. Landsat-8: Science and product vision for terrestrial global change research. Remote Sens. Environ. 2014, 145, 154–172. [Google Scholar] [CrossRef] [Green Version]
  35. Drusch, M.; Del Bello, U.; Carlier, S.; Colin, O.; Fernandez, V.; Gascon, F.; Hoersch, B.; Isola, C.; Laberinti, P.; Martimort, P.; et al. Sentinel-2: ESA’s optical high-resolution mission for GMES operational services. Remote Sens. Environ. 2012, 120, 25–36. [Google Scholar] [CrossRef]
  36. Claverie, M.; Ju, J.; Masek, J.G.; Dungan, J.L.; Vermote, E.F.; Roger, J.C.; Skakun, S.; Justice, C. The Harmonized Landsat and Sentinel-2 surface reflectance data set. Remote Sens. Environ. 2018, 219, 145–161. [Google Scholar] [CrossRef]
  37. Vermote, E.; Justice, C.; Claverie, M.; Franch, B. Preliminary analysis of the performance of the Landsat 8/OLI land surface reflectance product. Remote Sens. Environ. 2016, 185, 46–56. [Google Scholar] [CrossRef]
  38. Doxani, G.; Vermote, E.; Roger, J.C.; Gascon, F.; Adriaensen, S.; Frantz, D.; Hagolle, O.; Hollstein, A.; Kirches, G.; Li, F.; et al. Atmospheric correction inter-comparison exercise. Remote Sens. 2018, 10, 352. [Google Scholar] [CrossRef]
  39. Tucker, C.J. Red and photographic infrared linear combinations for monitoring vegetation. Remote Sens. Environ. 1979, 8, 127–150. [Google Scholar] [CrossRef] [Green Version]
  40. Jiang, Z.; Huete, A.R.; Didan, K.; Miura, T. Development of a two-band enhanced vegetation index without a blue band. Remote Sens Environ. 2008, 112, 3833–3845. [Google Scholar] [CrossRef]
  41. Richardson, A.J.; Wiegand, C.L. Distinguishing vegetation from soil background information. Photogramm. Eng. Remote Sens. 1977, 43, 1541–1552. [Google Scholar]
  42. Molod, A.; Takacs, L.; Suarez, M.; Bacmeister, J. Development of the GEOS-5 atmospheric general circulation model: Evolution from MERRA to MERRA2. Geosci. Model Dev. 2015, 8, 1339–1356. [Google Scholar] [CrossRef]
  43. Santamaria-Artigas, A.E.; Franch, B.; Guillevic, P.; Roger, J.-C.; Vermote, E.F.; Skakun, S. Evaluation of Near-Surface Air Temperature from Reanalysis Over the United States and Ukraine: Application to Winter Wheat Yield Forecasting. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2019. [Google Scholar] [CrossRef]
  44. Hatfield, J.L.; Prueger, J.H. Temperature extremes: Effect on plant growth and development. Weather Clim. Extrem. 2015, 10, 4–10. [Google Scholar] [CrossRef] [Green Version]
  45. Skakun, S.; Franch, B.; Vermote, E.; Roger, J.-C.; Becker-Reshef, I.; Justice, C.; Kussul, N. Early season large-area winter crop mapping using MODIS NDVI data, growing degree days information and a Gaussian mixture model. Remote Sens. Environ. 2017, 195, 244–258. [Google Scholar] [CrossRef]
  46. Lavreniuk, M.; Kussul, N.; Skakun, S.; Shelestov, A.; Yailymov, B. Regional retrospective high resolution land cover for Ukraine: Methodology and results. In Proceedings of the 2015 IEEE International Geoscience and Remote Sensing Symposium (IGARSS), Milan, Italy, 26–31 July 2015; IEEE: Piscataway, NJ, USA, 2015; pp. 3965–3968. [Google Scholar]
  47. Bishop, C.M. Pattern Recognition and Machine Learning; Springer: New York, NY, USA, 2006. [Google Scholar]
  48. Alemu, W.G.; Henebry, G.M. Land surface phenologies and seasonalities using cool earthlight in mid-latitude croplands. Environ. Res. Lett. 2013, 8, 045002. [Google Scholar] [CrossRef]
  49. McMaster, G.S.; Wilhelm, W.W. Growing degree-days: One equation, two interpretations. Agric. For. Meteorol. 1997, 87, 291–300. [Google Scholar] [CrossRef]
  50. Nguyen, L.H.; Joshi, D.R.; Clay, D.E.; Henebry, G.M. Characterizing land cover/land use from multiple years of Landsat and MODIS time series: A novel approach using land surface phenology modeling and random forest classifier. Remote Sens. Environ. 2018. [Google Scholar] [CrossRef]
  51. Sun, L.; Gao, F.; Anderson, M.; Kustas, W.; Alsina, M.; Sanchez, L.; Sams, B.; McKee, L.; Dulaney, W.; White, W.; et al. Daily mapping of 30 m LAI and NDVI for Grape yield prediction in California Vineyards. Remote Sens. 2017, 9, 317. [Google Scholar] [CrossRef]
  52. Manjunath, K.R.; Potdar, M.B.; Purohit, N.L. Large area operational wheat yield model development and validation based on spectral and meteorological data. Int. J. Remote Sens. 2002, 23, 3023–3038. [Google Scholar] [CrossRef]
  53. He, M.; Kimball, J.; Maneta, M.; Maxwell, B.; Moreno, A.; Beguería, S.; Wu, X. Regional crop gross primary productivity and yield estimation using fused Landsat-MODIS data. Remote Sens. 2018, 10, 372. [Google Scholar] [CrossRef]
  54. Monteith, J.L. Solar radiation and productivity in tropical ecosystems. J. Appl. Ecol. 1972, 9, 747–766. [Google Scholar] [CrossRef]
  55. Prince, S.D.; Goward, S.N. Global primary production: A remote sensing approach. J. Biogeogr. 1995, 815–835. [Google Scholar] [CrossRef]
  56. Tucker, C.J.; Holben, B.N.; Elgin, J.H., Jr.; McMurtrey, J.E., III. Relationship of spectral data to grain yield variation. Photogramm. Eng. Remote Sens. 1980, 46, 657–666. [Google Scholar]
  57. Hoerl, A.E.; Kennard, R.W. Ridge regression: Biased estimation for nonorthogonal problems. Technometrics 1970, 12, 55–67. [Google Scholar] [CrossRef]
  58. Phillips, D.L. A technique for the numerical solution of certain integral equations of the first kind. J. ACM (JACM) 1962, 9, 84–97. [Google Scholar] [CrossRef]
  59. Tikhonov, A.N.; Leonov, A.S.; Yagola, A.G. Nonlinear Ill-Posed Problems; Chapman & Hall: London, UK, 1998. [Google Scholar]
  60. Kussul, N.; Skakun, S.; Shelestov, A.; Lavreniuk, M.; Yailymov, B.; Kussul, O. Regional scale crop mapping using multi-temporal satellite imagery. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 2015, 40, 45–52. [Google Scholar] [CrossRef]
  61. Shelestov, A.; Kolotii, A.; Skakun, S.; Baruth, B.; Lozano, R.L.; Yailymov, B. Biophysical parameters mapping within the SPOT-5 Take 5 initiative. Eur. J. Remote Sens. 2017, 50, 300–309. [Google Scholar] [CrossRef]
Remotesensing 11 01768 g001 550
Figure 1. Location of Kirovohradska oblast in Ukraine and its division into districts. Shown also is the location of Sentinel-2 tiles covering the area. The background image shows a winter crop mask derived from Landsat 8 and Sentinel-2 for 2018.
Figure 1. Location of Kirovohradska oblast in Ukraine and its division into districts. Shown also is the location of Sentinel-2 tiles covering the area. The background image shows a winter crop mask derived from Landsat 8 and Sentinel-2 for 2018.
Remotesensing 11 01768 g001
Remotesensing 11 01768 g002 550
Figure 2. General workflow for winter wheat mapping and yield assessment.
Figure 2. General workflow for winter wheat mapping and yield assessment.
Remotesensing 11 01768 g002
Remotesensing 11 01768 g003 550
Figure 3. (a) Average number of cloud-free observations for winter crop pixels depending on satellite data usage. The number of pixels was taken from March until the end of June, which was the period of winter crop growth. (b) Spatial distribution of the number of cloud-free observations from Landsat 8 and Sentinel-2 from March through the end of June in 2018.
Figure 3. (a) Average number of cloud-free observations for winter crop pixels depending on satellite data usage. The number of pixels was taken from March until the end of June, which was the period of winter crop growth. (b) Spatial distribution of the number of cloud-free observations from Landsat 8 and Sentinel-2 from March through the end of June in 2018.
Remotesensing 11 01768 g003
Remotesensing 11 01768 g004 550
Figure 4. Peak and accumulation approaches for winter wheat yield assessment. Shown is the difference vegetation index (DVI), derived from Landsat 8 and Sentinel-2 data over a winter crop pixel, and fitted against accumulated growing degree days (AGDD) using a quadratic relationship.
Figure 4. Peak and accumulation approaches for winter wheat yield assessment. Shown is the difference vegetation index (DVI), derived from Landsat 8 and Sentinel-2 data over a winter crop pixel, and fitted against accumulated growing degree days (AGDD) using a quadratic relationship.
Remotesensing 11 01768 g004
Remotesensing 11 01768 g005 550
Figure 5. DVI (a) and NIR-red (b) dynamics for the same district over three years with variations in winter wheat yields.
Figure 5. DVI (a) and NIR-red (b) dynamics for the same district over three years with variations in winter wheat yields.
Remotesensing 11 01768 g005
Remotesensing 11 01768 g006 550
Figure 6. Comparison between winter crop areas from official statistics and satellite at the district scale in Kirovohradska oblast in 2016 to 2018.
Figure 6. Comparison between winter crop areas from official statistics and satellite at the district scale in Kirovohradska oblast in 2016 to 2018.
Remotesensing 11 01768 g006
Remotesensing 11 01768 g007 550
Figure 7. Regression crop yield models based on DVI and AUC for 2016 (a), 2017 (b), and 2018 (c).
Figure 7. Regression crop yield models based on DVI and AUC for 2016 (a), 2017 (b), and 2018 (c).
Remotesensing 11 01768 g007
Remotesensing 11 01768 g008 550
Figure 8. Performance of the AUC-DVI empirical winter wheat yield model on calibration data depending on the combined use of Landsat 8/OLI and Sentinel-2/MSI and a single usage. The coefficient of determination R2 (a) and relative RMSE (b) are shown.
Figure 8. Performance of the AUC-DVI empirical winter wheat yield model on calibration data depending on the combined use of Landsat 8/OLI and Sentinel-2/MSI and a single usage. The coefficient of determination R2 (a) and relative RMSE (b) are shown.
Remotesensing 11 01768 g008
Remotesensing 11 01768 g009 550
Figure 9. Estimates vs. reference winter wheat yields within the temporal cross-validation for the best AUC, coefficients–NIR + red model.
Figure 9. Estimates vs. reference winter wheat yields within the temporal cross-validation for the best AUC, coefficients–NIR + red model.
Remotesensing 11 01768 g009
Table
Table 1. Summary of the types of regression models used in the study.
Table 1. Summary of the types of regression models used in the study.
ParametersMetric Number   of   Regressors ,   n r Examples
Single VIPeak1 y y i e l d = ω 0 + ω 1 p D V I
Single VI or SRAUC1 y y i e l d = ω 0 + ω 1 s D V I
y y i e l d = ω 0 + ω 1 s N I R
Multiple SRsAUC2–4 y y i e l d = ω 0 + ω 1 s N I R + ω 2 s r e d
Multiple SRsa0, a1, a2, AUC4–16 y y i e l d = ω 0 + ω 1 a 0 , N I R        + ω 2 a 1 , N I R + ω 3 a 2 , N I R + ω 4 s N I R        + ω 5 a 0 , r e d        + ω 6 a 1 , r e d + ω 7 a 2 , r e d + ω 8 s r e d
Table
Table 2. Summary (average and standard deviation) of metrics (RMSE and R2), when fitting satellite-derived VI and SR values against AGDD using a quadratic model for 2016 to 2018. The percentages show a fraction of winter crop pixels with >6 valid HLS observations, for which metrics were calculated.
Table 2. Summary (average and standard deviation) of metrics (RMSE and R2), when fitting satellite-derived VI and SR values against AGDD using a quadratic model for 2016 to 2018. The percentages show a fraction of winter crop pixels with >6 valid HLS observations, for which metrics were calculated.
2016 (47%, Total 8,207,432 Winter Crop Pixels)
ValueRMSEAverage ysatR2
av.std.av.std.av.std.
DVI0.0160.0120.1840.0560.940.09
EVI20.0250.0200.3930.0830.940.10
NDVI0.0330.0170.6030.1070.940.09
Green0.0050.0020.0530.0070.610.25
Red0.0060.0030.0480.0110.810.16
NIR0.0170.0110.2320.0500.900.12
SWIR10.0150.0070.1810.0260.760.20
SWIR20.0150.0080.1300.0310.850.14
2017 (70%, Total 9,137,931 Winter Crop Pixels)
ValueRMSEAverage ysatR2
av.std.av.std.av.std.
DVI0.0140.0090.1680.0510.930.10
EVI20.0220.0150.3500.0760.930.11
NDVI0.0270.0150.5780.1060.940.11
Green0.0040.0020.0530.0060.600.26
Red0.0040.0020.0500.0110.830.19
NIR0.0140.0090.2190.0450.910.11
SWIR10.0110.0060.1860.0270.780.21
SWIR20.0100.0070.1380.0320.880.15
2018 (85%, Total 10,334,332 Winter Crop Pixels)
ValueRMSEAverage ysatR2
av.std.av. std.av.std.
DVI0.0150.0100.1810.0590.880.13
EVI20.0250.0160.3490.0890.880.13
NDVI0.0300.0170.6160.1190.900.12
Green0.0050.0020.0520.0070.590.27
Red0.0050.0030.0500.0120.770.21
NIR0.0150.0090.2300.0520.830.17
SWIR10.0130.0070.1800.0300.720.24
SWIR20.0120.0080.1290.0360.790.20
Table
Table 3. Winter wheat yield model performance (R2, RMSE, RRMSE) on calibration data depending on VI and peak/AUC feature for 2016 to 2018. Best models for each year are marked with italics.
Table 3. Winter wheat yield model performance (R2, RMSE, RRMSE) on calibration data depending on VI and peak/AUC feature for 2016 to 2018. Best models for each year are marked with italics.
ModelR2RMSE, t/haRRMSE, %p-Value
2016
Peak-DVI (data)0.1790.3087.75.61*10−2
Peak-DVI (fitting)0.3320.2787.06.29*10−3
AUC-DVI0.5880.2185.55.02*10−5
Peak-EVI2 (data)0.0560.3308.33.03*10−1
Peak-EVI2 (fitting)0.2820.2887.21.32*10−2
AUC-EVI20.2090.3027.63.71*10−2
Peak-NDVI (data)0.0880.3258.11.92*10−1
Peak-NDVI (fitting)0.4850.2446.14.55*10−4
AUC-NDVI0.0570.3308.32.98*10−1
2017
Peak-DVI (data)0.4220.2477.12.63*10−3
Peak-DVI (fitting)0.4000.2527.23.67*10−3
AUC-DVI0.5890.2086.01.26*10−4
Peak-EVI2 (data)0.4050.2517.23.40*10−3
Peak-EVI2 (fitting)0.3810.2567.34.89*10−3
AUC-EVI20.5700.2136.11.87*10−4
Peak-NDVI (data)0.3880.2547.34.38*10−3
Peak-NDVI (fitting)0.3930.2537.34.05*10−3
AUC-NDVI0.4070.2507.23.28*10−3
2018
Peak-DVI (data)0.5970.1764.74.53*10−4
Peak-DVI (fitting)0.5710.1824.87.08*10−4
AUC-DVI0.6080.1744.63.66*10−4
Peak-EVI2 (data)0.5650.1834.97.81*10−4
Peak-EVI2 (fitting)0.5070.1955.21.97*10−3
AUC-EVI20.5320.1905.11.34*10−3
Peak-NDVI (data)0.4060.2145.77.92*10−3
Peak-NDVI (fitting)0.3490.2246.01.60*10−2
AUC-NDVI0.2020.2486.68.07*10−2
Table
Table 4. Results of the cross-validation of different models. For each group, the best models are marked in italics.
Table 4. Results of the cross-validation of different models. For each group, the best models are marked in italics.
Spatial CVTemporal CV
ModelRMSE, t/haRRMSE, %R2RMSE, t/haRRMSE, %R2
VI-based and AUC
AUC–DVI0.2266.00.650.2576.90.60
AUC–NDVI0.3348.90.240.40810.90.15
AUC–EVI20.2717.20.500.3238.60.45
AUC–NIR0.2266.00.650.2366.30.63
AUC–red0.39610.60.180.47912.80.31
AUC–green0.3689.80.100.40810.90.00
AUC–SWIR10.38810.30.010.45912.30.19
SR-based and AUC
AUC–NIR + red0.2296.10.640.2536.70.60
AUC–NIR + green0.2296.10.640.2446.50.62
AUC–NIR + SWIR10.2286.10.650.2496.60.61
AUC–red + green0.2897.70.430.42711.40.35
AUC–red.+ SWIR10.3379.00.240.3569.50.14
AUC–green + SWIR10.3579.50.150.45712.20.04
AUC–NIR + red + green + SWIR10.2376.30.620.2687.10.53
VI/SR-based, and a0, a1, a2, AUC
AUC, coefficients–DVI0.2175.80.680.2185.80.68
AUC, coefficients–NDVI0.2837.60.450.3288.80.30
AUC, coefficients–EVI20.2727.30.500.33690.38
AUC, coefficients–NIR + red0.2075.50.710.2015.40.73
AUC, coefficients–NIR + green0.2185.80.680.2336.20.63
AUC, coefficients–NIR + SWIR10.2225.90.670.2215.90.67
AUC, coefficients–red + green0.2496.60.580.3669.80.53
AUC, coefficients–red + SWIR10.2837.60.450.2917.80.43
AUC, coefficients–green + SWIR10.3599.60.160.48613.00.03
AUC, coefficients–NIR + red + green + SWIR10.2125.70.700.2185.80.73

Share and Cite

MDPI and ACS Style

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. https://doi.org/10.3390/rs11151768

AMA Style

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 Sensing. 2019; 11(15):1768. https://doi.org/10.3390/rs11151768

Chicago/Turabian Style

Skakun, Sergii, Eric Vermote, Belen Franch, Jean-Claude Roger, Nataliia Kussul, Junchang Ju, and Jeffrey Masek. 2019. "Winter Wheat Yield Assessment from Landsat 8 and Sentinel-2 Data: Incorporating Surface Reflectance, Through Phenological Fitting, into Regression Yield Models" Remote Sensing 11, no. 15: 1768. https://doi.org/10.3390/rs11151768

APA Style

Skakun, S., Vermote, E., Franch, B., Roger, J. -C., Kussul, N., Ju, J., & Masek, J. (2019). Winter Wheat Yield Assessment from Landsat 8 and Sentinel-2 Data: Incorporating Surface Reflectance, Through Phenological Fitting, into Regression Yield Models. Remote Sensing, 11(15), 1768. https://doi.org/10.3390/rs11151768

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop