Mapping Winter Wheat Biomass and Yield Using Time Series Data Blended from PROBA-V 100- and 300-m S1 Products

Abstract

Monitoring crop areas and yields is crucial for food security and agriculture management across the world. In this paper, we mapped the biomass and yield of winter wheat using the new Project for On-Board Autonomy-Vegetation (PROBA-V) products in the North China Plain (NCP). First, the daily 100-m land surface reflectance was generated by fusing the PROBA-V 100-m and 300-m S1 products. Our results show that the blended data exhibited high correlations with the referenced data (0.71 ≤ R² ≤ 0.94 for the red band, 0.50 ≤ R² ≤ 0.95 for the near-infrared band, and 0.88 ≤ R² ≤ 0.97 for the shortwave infrared band). The time-series Normalized Difference Vegetation Index (NDVI) derived from the synthetic reflectance was then clustered for winter wheat identification. The overall classification accuracy was between 78% and 87%, with a kappa coefficient above 0.57, which was 10%–20% higher than the classification accuracy using the 300-m data. Finally, a light use efficiency model was employed to estimate the biomass and yield. The estimation results were closely related to the field-measured biomass and yield, with high R² and low root mean square errors (RMSE) (0.864 ≤ R² ≤ 0.871 and 168 ≤ RMSE ≤ 191 g/m² for biomass; and 0.631 ≤ R² ≤ 0.663 and 41.8 ≤ RMSE ≤ 62.8 g/m² for yield). This paper shows the strong potential of using PROBA-V 100-m data to enhance the spatial resolution of PROBA-V 300-m data and because the proposed framework in this study was based only on the relatively high spatio-temporal resolution PROBA-V data and achieved favorable results, it provides a novel approach for crop areas and yields estimation utilizing the relatively new data set.

1. Introduction

Food is necessary for human survival and development, and access to sufficient food is a basic right every person should enjoy. However, because of constraints imposed by multiple factors, such as population growth and poverty, food security is still one of the most important challenges in many parts of the world [1,2]. The latest FAO report indicates that although the trend in global hunger reduction continues, approximately one in nine people across the world remains suffered from nutrient deficiency in 2012–2014, particularly in developing countries [3]. Wheat is the world’s third largest food crop in terms of production, and it is traded frequently between different countries to satisfy the growing demands of food-deficit countries [4]. Therefore, the accurate and timely estimation of wheat areas and yields is important so that policy-making departments and international communities can make reasonable decisions.

Traditional crop biomass and yield estimation mainly involves time-consuming and laborious field surveys. Recently, with the development of remote sensing techniques, different methods have been developed to estimate biomass and yields based on satellite data. One approach develops statistical relationships between the satellite-derived vegetation indices (VIs) and field measured biomass or yields [5,6,7]. These regression models are widely used because they are simpler and more convenient than multivariate statistical methods; nevertheless, they are limited to local applications and are difficult to extend because they lack a theoretical basis. The other approach applies satellite data to calibrate crop growth models that dynamically simulate the physical processes of crop growth and management, whereby energy and nutrients are converted into biomass [8,9]. However, in consideration of the initial conditions and the numerous required input parameters, a compromise between physical correctness and data unavailability must be accepted for crop growth models. Taking the difficulties of the two above-mentioned approaches into consideration, the theory proposed by Monteith, which is based on absorbed photosynthetically active radiation (APAR) and light use efficiency (LUE), has been widely applied to terrestrial ecosystem net primary productivity and crop production estimations at both regional and global scales [10,11,12].

In the past few decades, moderate-spatial resolution instruments such as the Advanced Very High Resolution Radiometer (AVHRR), SPOT VEGETATION (SPOT-VGT) and Moderate Resolution Imaging Spectroradiometer (MODIS) have provided a large amount of space-sourced observational data that, given their consecutive and near real-time measurements, are both valuable and obligatory for global vegetation and land-cover dynamics monitoring [13,14]. However, after nearly 15 years of successful operations, the European SPOT-VGT stopped working in May 2014, a period in which the subsequent Sentinel-3 platform was still being prepared. To ensure continuity with the SPOT-VGT mission, the Belgian-led Project for On-Board Autonomy-Vegetation (PROBA-V) satellite, which has a finer spatial resolution, was launched in early May 2013 and released data freely from mid-March of 2014 [15]. PROBA-V has been proven to be well consistent with SPOT-VGT and has shown better performance in several areas such as crop identification and cropland mapping, opening up a plethora of new possibilities to scientists and researchers for land surface exploration at more sophisticated levels [16,17,18].

Data fusion is a valuable tool when high spatial and temporal resolution data are insufficient because of their long revisit cycles and the influence of unfavorable weather conditions in regions with fragmented landscapes [19,20,21]. Recently, several spatial and temporal fusion models have been proposed and proven to be practicable. Gao et al. used the Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM) to predict daily surface reflectance data at the Landsat spatial resolution and the MODIS temporal resolution [22]. To eliminate shortcomings of STARFM, several subsequent algorithms based on the algorithm have been presented. For example, taking crop phenology changes into consideration, Meng et al. developed the Spatial and Temporal Adaptive Vegetation index Fusion Model (STAVFM) to blend NDVIs extracted from the China Environment satellite (HJ-1) and MODIS data [23]. To overcome STARFM’s inaccurate predictions of the surface reflectance over heterogeneous landscapes, the Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model (ESTARFM) was proposed and has proven to be extremely effective [24].

The main objective of our study was to map the biomass and yield of winter wheat using PROBA-V 100- and 300-m S1 products. To achieve this, three major procedures were used: (1) generate daily 100-m land surface reflectance data by fusing PROBA-V 100- and 300-m S1 data using the STARFM and ESTARFM algorithms; (2) identify wheat fields based on time-series Normalized Difference Vegetation Index (NDVI) clustering and a priori knowledge; and (3) estimate biomass and yield by inverting time-series remote sensing data into a widely used light use efficiency model.

In this paper, Section 2 contains an introduction to the study areas and data; Section 3 discusses the data pre-processing and fusion, crop classification, biomass and yield estimation, and result evaluation strategies; Section 4 includes the results and accuracy assessments from Section 3; and Section 5 and Section 6 contain the discussion and conclusion of the paper, respectively.

2. Study Area and Data

2.1. Study Area

The study was conducted in areas of two counties in the southern North China Plain (NCP) (Figure 1). One is situated in Yucheng County, Shandong Province (Site 1, centered at 116°39′E, 36°26′N). Yucheng is an alluvial plain in the lower reaches of the Yellow River; it has an average altitude of 28 m above sea level and enjoys irrigated agriculture. It is characterized by a mean annual temperature of 13.2 °C and precipitation of 528 mm, approximately 60% of which occurs in July and August. The other study area is located in Guantao County, Hebei Province (Site 2, centered at 115°23′E, 36°37′N). Here, arable land accounts for approximately 70% of the county’s total area. Guantao is dominated by a typical warm, temperate, semi-humid continental monsoon climate, and the yearly frost-free period and sunshine duration are approximately 200 days and 2557 h, respectively. The Zhang and Yunwei Rivers, which lie across the county, provide irrigation water for agriculture. The study areas were chosen because they are typical planting regions of winter wheat, which ranks third among the most produced food crops in China.

The traditional tillage practice in the study areas is dominated by the two-season cropping of winter wheat and summer maize. In general, wheat is sown in mid or late October and harvested in late May or early June of the following year, and maize is planted in mid or late June and harvested in late September or early October.

2.2. Data Sources

2.2.1. Field Measured Data

Field experiments were conducted three times in 2015 during the crucial growing stages of winter wheat (

Table 1. The observations of winter wheat in the 2014–2015 agricultural year.

Study Sites	Crop Type	Sample Date	Growing Stage	Number (N)
Site 1	Winter WHEAT	16 April 2015	Booting	25
		17 May 2015	Flowering	25
		4 June 2015	Harvest	25
Site 2	Winter Wheat	15 April 2015	Booting	15
		17 May 2015	Flowering	15
		6 June 2015	Harvest	15

). The observations mainly included dry aboveground biomass and yield. Only large and relatively homogeneous fields (>1 km²) were selected, and all measurements were made 200 m or farther from the field boundaries. In addition, a Global Positioning System (GPS, Trimble GeoExplorer 6000 Series GeoXH, Trimble Navigation, Ltd., Sunnyvale, CA, USA) was used to record the coordinates so that they could be accurately located in the satellite images.

The total aboveground biomass (including stems, leaves and ears) was collected from five randomly selected small plots (1 m × 1 m) at each sample site in the center of the field. All the plant samples were heated to 105 °C and oven dried at 80 °C until they reached a constant weight in the laboratory. The final dry weight (DW) of the samples was recorded. The DW was then divided by the sample area and converted to g/m². The yield was measured following similar steps when the wheat was harvested; the difference was that only the grain portion was required for the yield calculation. In addition, we spent approximately two weeks in May conducting field surveys, during which the land use types, crop types and phenology were investigated.

2.2.2. Satellite Data

The PROBA-V satellite has three co-aligned cameras. It can provide data at 300- and 100-m spatial resolution. The 300-m products combine the image stripes of the three cameras into a swath of 2295 km and a near-daily global coverage (approximately 90%), whereas the 100-m products were obtained solely from the central camera, which covers a much smaller swath (approximately 517 km) and has full global coverage every five days. Specific and detailed information of the PROBA-V mission was previously described [25,26,27]. The PROBA-V 100- and 300-m S1 (daily synthesis product) top of canopy (TOC) reflectance products were downloaded from VITO’s product distribution portal [28] in the HDF5 file format and then converted to the ENVI file format using SPIRITS (Software for the Processing and Interpretation of Remotely sensed Image Time Series [29,30,31].

Because of the lower latitude of the study areas, PROBA-V 100-m S1 products are available for approximately every 4–5 days. The quality of the PROBA-V data is largely influenced by severe clouds and haze in the NCP. The corresponding quality status map dataset (a quality state indicator for each pixel) indicates that only approximately 42% of the 300-m images during the whole growing season of winter wheat (October to June in the following year) had more than 90% good pixels, whereas the percentage was much lower for the 100-m data. Therefore, there are a lot of gaps (spatial) over the study areas (same for most of high latitude areas) if only the PROBA-V 100-m imageries are used for the classification and yield estimation. All the 100-m images of the 2014–2015 agricultural season covering the two study sites with good pixels exceeded 96% were selected and processed (

Table 2. The Project for On-Board Autonomy-Vegetation (PROBA-V) 100-m S1 product used in the study.

Crop Types	Site 1		Site 2
	Date	Good Pixels Ratio	Date	Good Pixels Ratio
Winter wheat	2 November 2014	100%	2 November 2014	100%
	13 December 2014	100%	13 December 2014	100%
	18 December 2014	100%	31 December 2014	100%
	22 December 2014	100%	9 January 2015	100%
	31 December 2014	100%	10 February 2015	100%
	19 January 2015	99%	27 March 2015	100%
	10 February 2015	100%	10 April 2015	100%
	28 March 2015	98%	15 April 2015	100%
	10 April 2015	99%	12 May 2015	97%
	15 April 2015	98%	25 May 2015	99%
	7 May 2015	98%	30 May 2015	100%
	16 May 2015	100%	-	-
	30 May 2015	98%	-	-

Note: the images marked in bold were employed as references to evaluate the predicted images.

), which mainly included trimming the study area borders and projecting to the Universal Transverse Mercator (UTM) Zone 50N projection.

2.2.3. Meteorological Data

Both daily mean temperature and radiation are required as input for the light use efficiency model. Daily air temperature and sunshine duration data from the nearby climate station were obtained from the China Meteorological Data Sharing Service System (CMDSSS) [32]. The sunshine duration was used to calculate the photosynthetically active radiation (PAR) based on the FAO equations [33]. The highest daily mean temperature occurred in August and was below 30 °C. In addition, we did not take the influence of altitude changes on temperature into account because the study sites are flat.

3. Methods

Four main steps were performed for the mapping of winter wheat biomass and yield: Generating the time-series 100-m dataset by fusing the PROBA-V 100-m and 300-m data. Identifying the wheat fields based on the time-series NDVI clustering. Estimating biomass and yield combining a light use efficiency model with time-series NDVI and LSWI. Mapping biomass and yield based on the spatial distribution of wheat and the estimation results (Figure 2).

3.1. Daily 100-m Reflectance Dataset Generation

The STARFM and ESTARFM models were employed to blend the PROBA-V 100-m and 300-m S1 data because both of them have been widely and successfully used due to their ease of implementation and reasonable algorithm complexity [34,35,36,37]. We chose ESTARFM as the main fusion model because it was developed based on STARFM but overcomes STARFM’s shortcoming of making inaccurate predictions in heterogeneous landscapes such as that of our study area, which has various field sizes comparable to that of a 300-m pixel. Because the ESTARFM algorithm requires at least two pairs of 100-m/300-m images, STARFM was used when only one pair of images was available, for example, at the start and end of the growing season (i.e., in October 2014 and June 2015; Table 2). To minimize the uncertainty caused by human activity or environmental changes, the temporally closest available data were selected as two pairs in ESTARFM, and the images marked in bold in Table 2 were used as references to assess the predicted images [38]. The implementation of the two data fusion methods required resampling of the 300-m data to the spatial resolution of the 100-m imagery using a nearest neighbor algorithm. No further processing was applied because the provided PROBA-V data were pre-processed for atmospheric and geometric corrections. The main technical basis of the STARFM and the ESTARFM algorithms are briefly described as follows.

The STARFM can generate synthetic images based on a pair image of coarse-spatial-resolution and high-spatial-resolution (base images) and the coarse-spatial-resolution data on the prediction date through establishing a linear relationship between the base images of coarse-spatial-resolution and high-spatial-resolution. Assuming that land cover types keep stable over time, the data at high-spatial-resolution for date ($t_k)$can be predicted through using coarse-spatial-resolution observations for date ($t_k)$and one or several pairs of high-spatial-resolution and coarse-spatial-resolution images acquired at the same date ($t_0$).

$$H\left(x_i,y_j,t_k\right)= C\left(x_i,y_j,t_k\right)+ H\left(x_i,y_j,t_0\right)- C\left(x_i,y_j,t_0\right)$$

(1)

where H and C represent high-spatial-resolution and coarse-spatial-resolution image reflectance, (x, y) is the location of the predicted pixel reflectance, $t_0$ is the date of one input pair data, and $t_k$ is the date of predicted data.

During the implementation of ESTARFM algorithm, four major steps come about. Firstly, two high-spatial-resolution images are used to search for pixels similar to the central pixel in a local window. Secondly, the weight of each similar pixel $w_i$ was calculated by the spectral and spatial distance between similar pixel and its corresponding coarse-spatial-resolution pixel. Thirdly, the conversion coefficient $v_i$ can be computed from the linear regression of the high-spatial-resolution reflectance in the two observed pairs ($t_1$ and $t_2$) against the coarse-spatial-resolution reflectance of the similar pixels. Finally, $w_i$ and $v_i$ are used to calculate the high-spatial-resolution reflectance from the coarse-spatial-resolution image at the desired prediction date $t_p$, predicted as:

$$H\left(x,y,t_p\right)=H\left(x,y,t_0\right)+ {\displaystyle \sum }_{i=1}^N{w}_i\times v_i\times \left(C\left(x_i,y_i,t_p\right)-C\left(x_i,y_i,t_0\right)\right)$$

(2)

where H and C represent high-spatial-resolution image and coarse-spatial-resolution image reflectance, (x, y) is the location of the predicted pixel value while $x_i$ and $y_i$ is the location of ith similar pixel, and $t_0$ is the date of one input pair ($t_1$ or $t_2$). N is the total number of similar pixels of the predicted pixel within a local window. The two algorithms are explained in detail in Gao et al. [23] and Zhu et al. [24].

To derive parameters for the crop identification and yield estimation, the widely used NDVI and land surface water index (LSWI) were employed. NDVI is calculated based on the surface reflectance data of the red (RED) and near infrared (NIR) bands, and LSWI is based on the reflectance of the NIR and short-wave infrared (SWIR) bands [39,40,41]. To avoid spatial coverage gaps resulting from negative atmosphere conditions such as clouds and haze, the values of the pixels that were labeled as “good” during 10 days (10-day syntheses, the starting days were the 1st, 11th, or 21st day of a month) were averaged to composite the NDVI and LSWI. Finally, Savitzky–Golay (S-G) filtering was applied to remove the noise resulting from atmospheric variability and sun zenith angle changes from the two indices that [42,43].

$$\mathrm{NDVI}=\left({\rho }_{\mathrm{NIR}}-{\rho }_{\mathrm{RED}}\right)/\left({\rho }_{\mathrm{NIR}}+{\rho }_{\mathrm{RED}}\right)$$

(3)

$$\mathrm{LSWI}=\left({\rho }_{\mathrm{NIR}}-{\rho }_{\mathrm{SWIR}}\right)/\left({\rho }_{\mathrm{NIR}}+{\rho }_{\mathrm{SWIR}}\right)$$

(4)

where ${\rho }_{\mathrm{RED}}$, ${\rho }_{\mathrm{NIR}}$ and ${\rho }_{\mathrm{SWIR}}$ represent the reflectance of the red, near-infrared and short-wave infrared bands.

3.2. Crop Identification Based on Time-Series NDVI Clustering

According to ground investigations, wheat is the dominant crop in the study areas, accounting for 90% of the summer harvest crops. The human-managed wheat has different seasonality characteristics from single season natural vegetation and man-made features. Therefore, wheat areas can be identified by analyzing various phenological characteristics. Each wheat survey point was plotted as a time-series NDVI profile (Figure 2); the shaded areas indicate the differences between the maximum and minimum NDVIs for wheat. A large amount of research has proven that time-series images can provide improved classification accuracy compared with single images. Herein, a simple unsupervised Iterative Self-organizing Data Analysis Techniques Algorithm (ISODATA) and maximum likelihood (MLC) classification techniques, which have been used in several studies, were adopted for the time-series PROBA-V 300-m and the synthetic 100-m image clustering [16,44]. Firstly, the NDVI time series were created separately for the 300-m and 100-m dataset by stacking all the NDVI images covering the whole growing season of the individual crop. Afterwards, pixels with specific spatial intervals from the images were selected as training samples (a gridded set of samples) and the ISODATA algorithm was applied to those training samples to generate several classes. Each class (a group of many separated pixels) again serves as the training samples for the MLC method which was employed to the whole image to cluster the time-series NDVI. We classified the NDVI time series into 10 classes because only no more than ten land cover types existed and therefore 10 classes are basically enough for the identification of wheat fields. The time-series mean NDVI values of each class were then extracted and displayed as the NDVI temporal profile (10 profiles in total). Finally, the wheat fields were identified based on the classification results and NDVI temporal profiles in combination with the prior maximum and minimum NDVI curves (Figure 3). Specifically, if the profile of a corresponding class falls in the shaded areas of Figure 3, the class is considered to be wheat. Otherwise, it belongs to another land cover type.

3.3. Algorithms for Biomass and Yield Estimation

Based on the Monteith model, the accumulation of aboveground biomass is proportional to accumulated APAR [10,12]:

$$\mathrm{AGB}=\mathrm{R}\times \mathrm{LUE}\times {\displaystyle \sum }_{t=0}^N\left(\mathrm{APAR}\left(\mathrm{\Delta }t\right)\times \mathrm{\Delta }t\right)$$

(5)

$$\mathrm{APAR}={\displaystyle \sum }_{t=0}^N\left(\mathrm{PAR}\times \mathrm{FPAR}\right)\times \mathrm{\Delta }t$$

(6)

$$\mathrm{LUE}\left(x,t\right)= {\epsilon }^{*}\times T_{\epsilon 1}\left(x,t\right)\times T_{\epsilon 2}\left(x,t\right)\times W_{\epsilon }\left(x,t\right)$$

(7)

where AGB is the aboveground biomass during the growing season ($\mathrm{\Delta }t$), R is the proportion of total biomass allocated to aboveground productivity, and N represents the day when the experiment was carried out. Here, ε* is the maximum LUE for green plants under optimal conditions. $T_{\epsilon 1}$, $T_{\epsilon 2}$ and $W_{\epsilon }$ are environmental stressors that limit ε* under unfavorable environmental conditions, $T_{\epsilon 1}$ and $T_{\epsilon 2}$represent temperature stress and have been described in detail previously [11], and $W_{\epsilon }$ reflects the water stress on the photosynthesis of crops. In general, $W_{\epsilon }$is calculated based on potential evapotranspiration (PET) and estimated evapotranspiration (EET). Our study used the approach presented by Xiao et al. for the vegetation photosynthesis model (VPM) to estimate $W_{\epsilon }$ [39,40]:

$$W_{\epsilon }\left(x,t\right) = \left(1+\mathrm{LSWI}\right)/\left(1+{\mathrm{LSWI}}_{\max }\right)$$

(8)

where ${\mathrm{LSWI}}_{\max }$ is the maximum LSWI within the crop growing season for individual pixels.

The fraction of PAR intercepted by the canopy (FPAR) is a key parameter that is usually estimated based on different VIs such as NDVI and simple ratio (SR). In this research, FPAR was calculated as a linear function of SR, following Sellers et al. [45]:

$$\mathrm{FPAR}= \frac{\left(\mathrm{SR}-{\mathrm{SR}}_{\min }\right)\times \left({\mathrm{FPAR}}_{\max }-{\mathrm{FPAR}}_{\min }\right)}{{\mathrm{SR}}_{\max }-{\mathrm{SR}}_{\min }}+{\mathrm{FPAR}}_{\min }$$

(9)

$$\mathrm{SR}=\mathrm{NIR}/\mathrm{RED}=\left(1+\mathrm{NDVI}\right)/\left(1-\mathrm{NDVI}\right)$$

(10)

where ${\mathrm{SR}}_{\max }$and ${\mathrm{SR}}_{\min }$ correspond to the 95th and 5th percentile of SR for all cropland, respectively, and ${\mathrm{FPAR}}_{\max }$ = 0.95, and ${\mathrm{FPAR}}_{\min }$ = 0.001.

Previous studies have found that biomass can be converted to crop yields based on a common method (Equation (11)) [46]. Despite the variations in the harvest index (HI) during different crop growing stages, we used the average value of HI during the entire season as the final HI. The averaged HI was defined as 0.45 for wheat based on former experiments in the study areas [23,47]. The specific crop-type parameters in the model are defined in

Table 3. The parameters used to estimate wheat biomass and yield.

Parameter	Description	Value	Unit
${\epsilon }^{*}$	the maximum light use efficiency	2.54	g·MJ−1 PAR
R	proportion of aboveground productivity	0.90	dimensionless
HI	the harvest index	0.45	dimensionless

.

$$\mathrm{Yield} ={\displaystyle \sum }_{n=N_0}^{n=N_1}{\mathrm{AGB}}_t\times \mathrm{HI}$$

(11)

where $N_0$ and $N_1$ are the arbitrary growing season start and end dates, respectively.

3.4. Results Evaluation Strategy

3.4.1. Evaluation of the Data Fusion Result

The accuracy assessment of the synthetic images was undertaken using a set of images independent from those used as ESTARFM inputs, and the five statistical criteria shown in

Table 4. Statistical indices for the assessment of the quality of the synthetic images.

Indicators	Formula
Determination coefficient (R2)	$1-{\displaystyle \sum }{\left(Y-Y_f\right)}^2/{\displaystyle \sum }{\left(Y-Y_m\right)}^2$
Root mean square errors (RMSE)	$\sqrt{{\displaystyle \sum }{\left(Y-Y_p\right)}^2/n}$
Relative RMSE (RRMSE)	$\frac{1}{Y_m}\times \sqrt{{\displaystyle \sum }{\left(Y-Y_p\right)}^2/n}$
Average absolute deviation (AAD)	${\displaystyle \sum }\left|Y-Y_p\right|/n$
Average deviation (AD)	${\displaystyle \sum }\left(Y-Y_p\right)/n$

Note: $Y$ represents the pixel values of the referenced images, $Y_m$ represents the mean values of the referenced images, $Y_p$ represents the values of the predicted images, and $Y_f$ represents the values of the linear fitting.

were chosen as the evaluation indicators. Since there are many negative values in the images of NDVI and LSWI, only the R², root mean square errors (RMSE), average absolute deviation (AAD) and average deviation (AD) were selected to assess the prediction accuracy of the two indices.

3.4.2. Assessment of Biomass and Yield Estimation

Due to the small amount of observational data (N < 100), it was quite inefficient to withhold part of the data to validate the crop biomass and yield estimation results; therefore, the leave-one-out cross-validation (LOOCV) approach was chosen to examine and validate the accuracy of the estimation models [48]. The LOOCV method involves using one observation as the validation sample and the remaining observations as the training samples. This procedure was repeated N times. The R², RMSE and relative RMSE (RRMSE) of the algorithm were then estimated by averaging the values obtained from the N iterations.

$$\mathrm{RMSE}=\sqrt{\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\left(O_i-E_i\right)}^2}$$

(12)

$$\mathrm{RRMSE}=\frac{\mathrm{RMSE}}{\overline{O}}\times 100\%$$

(13)

where n is the number of observations, $O_i$ is the observed value, $E_i$ is the estimated value, and $\overline{O}$ is the mean value of all the observed values.

4. Results

4.1. The ESTARFM Prediction Results

Figure 4 shows the relationships between the observed PROBA-V 100-m and ESTARFM predicted reflectance values. All the slopes of the scatter plots were close to 1, and the intercepts were small, proving that ESTARFM accurately predicted the reflectance for each image. The R² values ranged from 0.712 to 0.931 for the RED band, from 0.503 to 0.941 for the NIR band and from 0.880 to 0.964 for the SWIR band, and the average deviations were almost negligible, as shown in

Table 5. Results for the pixel-based, band-by-band accuracy assessment of the referenced PROBA-V 100-m reflectance compared to the ESTARFM predicted reflectance.

Sites	Date (100-m)	Date (300-m)	Band	R2	RMSE	RRMSE	AAD	AD
Site 1	18 December 2014	18 December 2014	RED	0.804	0.0083	0.0775	0.0650	9.5 × 10−8
			NIR	0.899	0.0132	0.0536	0.0099	4.8 × 10−7
			SWIR	0.895	0.0062	0.0278	0.0047	3.6 × 10−7
Site 1	22 December 2014	22 December 2014	RED	0.760	0.0089	0.0838	0.0068	1.6 × 10−7
			NIR	0.860	0.0145	0.0624	0.0108	3.9 × 10−7
			SWIR	0.880	0.0061	0.0281	0.0046	9.6 × 10−8
Site 1	15 April 2015	15 April 2015	RED	0.893	0.0109	0.1238	0.0082	1.8 × 10−7
			NIR	0.917	0.0172	0.0482	0.0130	4.7 × 10−7
			SWIR	0.939	0.0082	0.0412	0.0056	1.2 × 10−6
Site 1	7 May 2015	7 May 2015	RED	0.931	0.0111	0.1432	0.0082	2.0 × 10−7
			NIR	0.915	0.0175	0.0464	0.0130	9.4 × 10−7
			SWIR	0.923	0.0086	0.0595	0.0061	3.0 × 10−7
Site 2	31 December 2014	31 December 2014	RED	0.900	0.0055	0.0522	0.0042	8.8 × 10−8
			NIR	0.941	0.0094	0.0308	0.0070	1.4 × 10−7
			SWIR	0.948	0.0097	0.0224	0.0067	3.5 × 10−7
Site 2	15 April 2015	15 April 2015	RED	0.930	0.0112	0.1365	0.0081	5.3 × 10−8
			NIR	0.932	0.0181	0.0443	0.0134	3.9 × 10−8
			SWIR	0.964	0.0074	0.0356	0.0051	9.4 × 10−8
Site 2	25 May 2015	25 May 2015	RED	0.712	0.0175	0.1813	0.0122	3.1 × 10−7
			NIR	0.503	0.0190	0.0592	0.0137	4.3 × 10−6
			SWIR	0.910	0.0115	0.0629	0.0082	3.8 × 10−7

. The higher R² values were achieved during the middle of the growing season (April), whereas the smaller R² values were observed at the early part of the growing season and at its end (in December and May, respectively). Notably, the predictions were generally best for the SWIR band and worst for the RED band, except for the predictions for 7 May at Site 1 and 25 May at Site 2. This result is consistent with numerous studies that the ESTARFM have shown better performance at the longer wavelength [49]. For the predicted images on 7 and 25 May, the accuracy was lower for the NIR than for the RED band, which was probably a result of the changes in land-cover from seasonal farming activities in the cropped areas at the end of the growing season, and the NIR band was more sensitive to these changes than the RED and SWIR bands.

Comparisons between the blended NDVI/LSWI and the observed 100-m NDVI/LSWI are provided in Figure 5. The accuracy assessment shows that most of the blended data were also closer to the 1:1 line. The R² between the predicted NDVI and the referenced NDVI was higher than 0.84 for Site 1 and 0.69 for Site 2, and the R² between the predicted LSWI and the referenced LSWI was larger than 0.84 for Site 1 and 0.86 for Site 2. The average absolute deviation was less than 0.06 for NDVI and 0.03 for LSWI (

Table 6. Results for the NDVI and LSWI accuracy assessment of the referenced PROBA-V 100-m images compared to the ESTARFM predicted images.

Sites	Date (100-m)	Date (300-m)	Indices	R2	RMSE	AAD	AD
Site 1	18 December 2014	18 December 2014	NDVI	0.878	0.0447	0.0344	5.9 × 10−7
			LSWI	0.892	0.0285	0.0203	4.3 × 10−8
Site 1	22 December 2014	22 December 2014	NDVI	0.846	0.0379	0.0496	2.4 × 10−7
			LSWI	0.840	0.0329	0.0238	5.1 × 10−8
Site 1	15 April 2015	15 April 2015	NDVI	0.924	0.0473	0.0361	6.1 × 10−7
			LSWI	0.957	0.0293	0.0215	1.7 × 10−7
Site 1	7 May 2015	7 May 2015	NDVI	0.937	0.0423	0.0316	1.6 × 10−6
			LSWI	0.955	0.0281	0.0204	5.7 × 10−7
Site 2	31 December 2014	31 December 2014	NDVI	0.950	0.0274	0.0270	3.5 × 10−7
			LSWI	0.957	0.0213	0.0156	2.9 × 10−8
Site 2	15 April 2015	15 April 2015	NDVI	0.944	0.0453	0.0335	1.1 × 10−6
			LSWI	0.973	0.0256	0.0181	2.5 × 10−7
Site 2	25 May 2015	25 May 2015	NDVI	0.699	0.0703	0.0509	3.7 × 10−7
			LSWI	0.860	0.0407	0.0298	2.0 × 10−8

). The lowest R² value for the predicted and referenced NDVI was observed on 25 May at Site 2 because of the inaccurate prediction of the NIR band. Generally speaking, the predictions were better for the LSWI than for the NDVI due to the more accurate predictions for the SWIR band than for the RED band. Furthermore, although data from similar sensors may have some advantages in terms of fusion, the blended images were not completely consistent with the actual images (as shown in Figure 4 and Figure 5), which possibly resulted from different viewing angles and land-surface changes. Overall, the predicted NDVI and LSWI were fairly reliable and could be used for the research discussed below.

4.2. Generation of Winter Wheat Maps

To assess the classification accuracy at different spatial resolutions, we classified the time-series PROBA-V 300-m and the predicted 100-m images. The NDVI profiles of the 10 classes are shown in Figure 6. Combined with Figure 3 and Figure 6, Classes 1–5 were considered to be wheat, whereas the remaining classes were considered to be other types. Classes 1–5 and Classes 6–10 were then merged into two larger classes to evaluate the accuracy of the classification. An interesting finding that is inconsistent with the previous assumption is that almost all 10 profiles followed a similar trend though the shapes were different, mixed pixels of the images and misclassifications of the clustering methods can well explain this phenomenon. The classification results are provided in Figure 7 and Figure 8. Some relatively small objects such as roads and rivers can be clearly identified in Figure 7 when compared with the classification shown in Figure 8 due to the higher spatial resolution of the 100-m data. The classification accuracy assessments were conducted using the ground investigation data (92 points for wheat and 104 points for the other types at Site 1; and 126 points for wheat and 154 points for the other types at Site 2).

The accuracy assessments, including the overall classification accuracies, producer’s accuracies, user’s accuracies and kappa coefficients, are presented in

Table 7. Classification accuracies of the 300-m and 100-m time-series data at Site 1.

Class	100-m		300-m
	Producer’s Accuracy	User’s Accuracy	Producer’s Accuracy	User’s Accuracy
Wheat	86.96%	81.63%	72.83%	65.05%
Others	82.69%	87.76%	65.38%	73.91%
Overall Accuracy: 84.69%; Kappa: 0.7198			Overall Accuracy: 68.88%; Kappa: 0.3795

and

Table 8. Classification accuracies of the 300-m and 100-m time-series data at Site 2.

Class	100-m		300-m
	Producer’s Accuracy	User’s Accuracy	Producer’s Accuracy	User’s Accuracy
Wheat	80.95%	73.91%	69.05%	63.04%
Others	76.62%	83.1%	66.88%	72.54%
Overall Accuracy: 78.57%; Kappa: 0.5708			Overall Accuracy: 67.86%; Kappa: 0.3562

. The overall accuracies and kappa coefficients of classification for the time-series blended data were 84.69% and 0.7198 at Site 1, respectively, and 78.57% and 0.5708 at Site 2, respectively. These values were noticeably higher than those of the 300-m data (for which the overall accuracies and kappa coefficients were 68.88% and 0.3795 at Site 1, respectively, and 67.86% and 0.3562 at Site 2, respectively), the existence of more mixed pixels in the 300 m data might explain the lower accuracies. Nevertheless, the classification accuracies were acceptable no matter for the 100-m or for the 300-m data. The producer’s accuracy of wheat at Site 1 was higher than that at Site 2 because a mountainous region at Site 2 led to many more fragmentized fields that were prone to misclassification. Moreover, the results also showed that the misclassification of wheat was lower than that of the other land-cover features, which probably occurred because the farmlands in the study areas were typically much larger and more homogeneous. This probably because the BRDF effects in the mountainous regions cause additional variations on the slopes. However, wheat fields are not inclined, so the BRDF effects are smaller there.

4.3. Mapping the Biomass and Yield

The winter wheat biomass and yield in Yucheng and Guantao were estimated using the generated time-series NDVI and LSWI datasets. Comparisons between the simulated and measured biomass are presented in Figure 9. The estimated biomass from the blended 100-m data shown in Figure 9a,b were generally in good agreement with the observed biomass (R² = 0.864, RMSE = 191 g/m² and RRMSE = 16.7% at Site 1; R² = 0.871, RMSE = 168 g/m² and RRMSE = 12.7% at Site 2); most of the scatter points are distributed along the fitting line, and the slopes are between 0.888 and 0.916. Figure 10 shows that the simulated yields were also correlated with the measured yields (R² = 0.663, RMSE = 62.8 g/m² and RRMSE = 7.27% for Site 1; R² = 0.631, RMSE = 41.8 g/m² and RRMSE= 4.88% for Site 2), which indicates that the fusion dataset can be used for crop biomass and yield estimations when the original PROBA-V 100-m data are limited. It also can be observed from Figure 9 that the measured biomass was more correlated with the estimations that involved the blended data (they had higher R², and lower RMSE and RRMSE) than were the 300-m products. This is because the 100-m data, with its higher spatial resolution, can detect spatial details that ignored by the 300-m data. Because the agricultural landscape in the NCP is extremely complex, monitoring crop yields using higher resolution remote sensing data is appropriate on a local scale.

The spatial distributions of the wheat biomass and yield (Figure 11 and Figure 12, respectively) can be obtained by combining the field maps derived in Section 4.2 and the estimation results. The biomass was greater than 1000 g/m² in mid-May, whereas the yields ranged from 750 to 1000 g/m² in most parts of the study areas. Figure 11a and Figure 12a indicate that the biomass and yield were much lower in the northern part of Yucheng, which has infertile farmland, whereas the southern region is a major grain-producing area in the small county due to its fertile soil and convenient agricultural facilities. In contrast, Figure 11b and Figure 12b demonstrate that the biomass and yield were relatively uniform in Guantao because of its excellent irrigation facilities. In addition, wheat production can be easily acquired from the yield distribution maps.

5. Discussion

In this paper, based on the existing PROBA-V 100- and 300-m S1 data, we were able to adequately map the biomass and yield of winter wheat. However, there are still several limitations and suggestions for this study in the near future, which will be discussed below.

5.1. Data Fusion Methods

Although the data fusion methods were appropriate to generate time-series 100-m data in this study, several limitations are still unavoidable, such as the STARFM algorithm may not perform so excellent in the complex heterogeneous areas and predictions from the ESTARFM may be inaccurate during a long term because the ratio of reflectance linear change is constant might be not suitable under this situation [22,24]. In addition, the proposed framework would be less effective in the regions with successive clouds because if the base high-spatial-resolution or coarse-spatial-resolution images are heavily influenced by unfavorable weather conditions, the predicted high-spatial-resolution data might not be so accurate and ultimately the crop classification and yield estimation results will be affected. Recently, remote sensing images from multiple satellite sensors is widely used to construct data time-series for agriculture monitoring. For example, Battude et al. used the Formosat-2, Landsat-8, Deimos-1 and SPOT-4 images to construct high temporal resolution data for the maize biomass and yield estimation in the southwest of France [50]. Hao et al. demonstrated the strong potential of time series data merged from Landsat-5 TM and HJ-1 CCD for crop classification in northern Xinjiang, China [51]. Siachalou et al. combined images from RapidEye and Landsat satellites to identify crop types in an irrigated agricultural region of the northern Greece [52]. Since this state-of-the-art method also performs very well in numerous applications and studies, it is necessary to compare performance of the two methods in the future studies to find which method is more suitable for the time series data construction.

5.2. Mixed Pixels

The potential and advantages of using high spatial and temporal remote sensing data to describe the spatiotemporal variability in crop biophysical variables have already been demonstrated [53,54]. Although the spatial resolution of the PROBA-V 100-m data improved when compared with the MODIS and SPOT-VGT data, mixed pixels still seem inevitable. For example, the effect of mixed pixels is a main reason that most of the NDVI profiles in Figure 5 follow a similar trend. In our study, this effect was not as serious because only larger and relatively homogeneous fields were chosen to evaluate the classification and estimation results. However, in regions with high farmland fragmentation, less accurate estimation results may be obtained because each pixel of PROBA-V reflectance data represents the averaged values of different crops. Fortunately, new satellites that combine higher spatial resolution with much shorter revisit cycles (such as the recently launched Sentinel-2A satellite and subsequent Sentinels), will weaken the influence of mixed pixels and, thus provide improved opportunities for crop classifications and yield estimations at field to sub-field scales [55,56]. Currently, the unmanned aerial vehicle (UAV) data is also an important source of high spatial resolution data and it will be tested to fuse with satellite data to gain a better monitoring of agriculture in the further studies [57].

5.3. LUE

Many researchers have demonstrated that LUE is positively related to crop growth and, ultimately, to biomass production [54,58]. LUE is approximately constant for forests and natural ecosystems, particularly for crops, when their growth is not limited by water or nutrient shortages or adverse climatic conditions [10,59]. However, in practical cases, water, nitrogen, and climate conditions are the main factors that impact wheat growth, all of which finally influence LUE. Sinclair and Muchow reviewed several studies and reported that the seasonal wheat LUE ranges from 0.73 to 1.62 g·MJ⁻¹ [60]. The LUE for wheat has been found to have tremendous variations. O’Connell et al. and Kiniry et al. revealed LUE values of 2.38 and 1.81 ± 0.05 g·MJ⁻¹ in Western Australia and New Zealand, respectively [61,62]. Duchemin et al. and Dong et al. adopted the SCEM-UA algorithm to calibrate LUE, which ranged between 1.5 and 3.5 g·MJ⁻¹ for winter wheat in central Morocco and Southern Canada, respectively [54,63]. Lobell et al. found that the LUE values in northwestern Mexico were 2.2–2.4 g·MJ⁻¹ and were similar to those determined in previous studies of wheat, which generally discovered LUEs in the range of 2–3 g·MJ⁻¹ [46,62,64]. Tao et al. calibrated LUE as 4.23 g·MJ⁻¹ with a temperature and moisture down-regulator and 3.92 g·MJ⁻¹ without a moisture down-regulator to estimate crop production across China [65]. Some studies have proven the significant dependence of LUE on photosynthesis pathways, such as those of C3 and C4 crops [46,66].

In this study, the adoption of a constant LUE value might be a potential source of error in the biomass and yield estimations because LUE has been proven to vary with different nutritional conditions and growth stages. Therefore, more field experiments over all the growth stages of wheat should be carried out to identify appropriate LUE values. Generally, wheat LUE can either be calculated as the ratio of the total aboveground dry weight to the intercepted radiation or can be measured using a portable gas exchange system, such as an LI-6200 (LI-COR, Inc., Lincoln, NE, USA) [59,67,68].

5.4. FPAR

FPAR is a crucial input of the light use efficiency model that can usually be obtained with the help of remote sensing data [69,70]. Because many studies have found that several VIs and FPAR dominated linear or nearly linear relationships, a simpler and common method was adopted to estimate the parameter in this study. To evaluate our FPAR estimation method, FPAR products from Copernicus (a European system for monitoring the earth [71] were chosen as reference data because ground-measured data were not available. We selected three key periods (mid-March, mid-April and mid-May), corresponding to the early, middle and late wheat growing stages, respectively, to validate the relationship between the Copernicus FPAR and our estimated FPAR. The results (Figure 13) indicated that there was a good agreement between the two products (0.51 ≤ R² ≤ 0.75, 0.05 ≤ RMSE ≤ 0.06 and 8.84% ≤ RRMSE ≤ 16.64%), which demonstrates the feasibility of this approach for FPAR estimation.

Despite acquiring satisfactory results, several limitations also exist, including the varying correlations between the VIs and FPAR for different crop types and growing stages. Meanwhile, VIs based on red and near infrared wavelengths such as NDVI and SR tend to gradually saturate at moderate-to-high vegetation coverages because light absorption and scattering are very strong in the red and near infrared regions [42,72]. Recently, many researchers have attempted to modify existing VIs or develop new VIs to reduce saturation and improve accuracy in FPAR estimations. The wide dynamic range vegetation index (WDRVI) and green chlorophyll index (CIgreen) have shown higher sensitivity to FPAR than NDVI over the entire growing season [73,74]. VIs that incorporate red-edge bands such as red-edge NDVI and MERIS terrestrial chlorophyll index (MTCI) also have been proven to be more sensitive to FPAR at high LAI stages [54]. In future research, FPAR can be calculated based on VIs that are more sensitive, and additional field experiments should be carried out to validate the results. Moreover, estimations based on radiative transfer mechanism models are another potential solution [75].

5.5. Meteorological Data

The final limitation of this research arises from using the low-spatial-resolution, site-observed meteorological data to represent the entire study area, which is a common limitation in almost all modeling approaches [8,54]. Generally, the representativeness of weather data is sufficient in flat and homogeneous regions. However, it must be carefully considered for the representativeness of the observation data in topographically irregular and complex areas. Nonetheless, this problem could be largely attenuated by the development of extrapolation or interpolation techniques for meteorological data [76]. Remote sensing techniques might also potentially be used to acquire meteorological variables at higher spatial resolution [77].

6. Conclusions

This study proposes a framework for biomass and yield mapping based on the new PROBA-V S1 products at different spatial resolutions. First, the STARFM and ESTARFM algorithms were employed to fuse PROBA-V 100-m and 300-m data. The results indicated that this fusion had excellent abilities because the blended reflectance was closely correlated with the observed reflectance (R² was higher than 0.71, 0.50 and 0.88 for the RED, NIR and SWIR bands, respectively). The NDVI and LSWI based on the predicted surface reflectance also had strong correlations with the referenced images (R² between 0.69 and 0.95 for NDVI and R² between 0.84 and 0.98 for LSWI). Then, winter wheat fields were identified based on time-series NDVI clustering, and favorable results were acquired. In comparison with the PROBA-V 300-m images, the overall accuracy of the classification improved by 10%–20% for the blended 100-m data. Finally, the biomass and yield were simulated using a simple light use efficiency model. The estimation results agreed well with the field measured values (0.864 ≤ R² ≤ 0.871, 168 ≤ RMSE ≤ 191 g/m² for biomass; 0.631 ≤ R² ≤ 0663, 41.8 ≤ RMSE ≤ 62.8 g/m² for yield). Considering the relatively high spatial and temporal resolutions of the PROBA-V images, our approach will be useful for the timely and accurate monitoring of crop planting areas and yields based on the new data.