Sentinel-1 Data for Winter Wheat Phenology Monitoring and Mapping

: The ability of Synthetic Aperture Radar (SAR) Sentinel-1 data to detect the main wheat phenological phases was investigated in the Bekaa plain of Lebanon. Accordingly, the temporal variation of Sentinel-1 (S1) signal was analyzed as a function of the phenological phases’ dates observed in situ (germination; heading and soft dough), and harvesting. Results showed that S1 data, unlike the Normalized Di ﬀ erence Vegetation Index (NDVI) data, were able to estimate the dates of theses phenological phases due to signiﬁcant variations in S1 temporal series at the dates of germination, heading, soft dough, and harvesting. Particularly, the ratio VV / VH at low incidence angle (32–34 ◦ ) was able to detect the germination and harvesting dates. VV polarization at low incidence angle (32–34 ◦ ) was able to detect the heading phase, while VH polarization at high incidence angle (43–45 ◦ ) was better than that at low incidence angle (32–34 ◦ ), in detecting the soft dough phase. An automated approach for main wheat phenological phases’ determination was then developed on the western part of the Bekaa plain. This approach modelled the S1 SAR temporal series by smoothing and ﬁtting the temporal series with Gaussian functions (up to three Gaussians) allowing thus to automatically detect the main wheat phenological phases from the sum of these Gaussians. To test its robustness, the automated method was applied on the northern part of the Bekaa plain, in which winter wheat is harvested usually earlier because of the di ﬀ erent weather conditions. The Root Mean Square Error (RMSE) of the estimation of the phenological phases’ dates was 2.9 days for germination, 5.5 days for heading, 5.1 days soft dough, 3.0 days for West Bekaa’s harvesting, and 4.5 days for North Bekaa’s harvesting. In addition, a slight underestimation was observed for germination and heading of West Bekaa ( − 0.2 and − 1.1 days, respectively) while an overestimation was observed for soft dough of West Bekaa and harvesting for both West and North Bekaa (3.1, 0.6, and 3.6 days, respectively). These results are encouraging, and thus prove that S1 data are powerful as a tool for crop monitoring, to serve enhanced crop management and production handling.


Introduction
Demanding and achieving sufficient and stable crop productivity is a global desire. Crop management must be conscientious to climatic variability and adapt its practices to the given conditions, whether socio-economic oriented, or environmentally related. From this perspective, winter wheat and potato occupy the largest areas of cultivated arable lands [57]. In November, after the first effective rain, wheat farmers launch their sowing procedure to eventually harvest in June-July the year after, occupying up to 12,000 ha annually of the plain. Over 80% of the wheat plots are supplementary irrigated. The frequency of the irrigation occurrence is linked directly with the rainfall received during the season. However, as normally the rain stops in February-March, farmers tend to irrigate for the sake of better grain yield at the end of the season. Nitrogen is the other main input to wheat. The most growth driving nutrient (Nitrogen) is supplied by wheat farmers as ammonium sulfate in amounts of up to 230 Kg/ha. However, the amounts and timing vary among farmers because of different fertilization management practices, as well as the crop rotation types.
Remote Sens. 2019, 11, x FOR PEER REVIEW 4 of 30 eventually harvest in June-July the year after, occupying up to 12,000 ha annually of the plain. Over 80% of the wheat plots are supplementary irrigated. The frequency of the irrigation occurrence is linked directly with the rainfall received during the season. However, as normally the rain stops in February-March, farmers tend to irrigate for the sake of better grain yield at the end of the season. Nitrogen is the other main input to wheat. The most growth driving nutrient (Nitrogen) is supplied by wheat farmers as ammonium sulfate in amounts of up to 230 Kg/ha. However, the amounts and timing vary among farmers because of different fertilization management practices, as well as the crop rotation types.

Optical Data
Fifty-eight optical images were acquired by the two twin satellites of Sentinel-2 (Sentinel-2A and Sentinel-2B). The acquired images cover a period from November 2017 to August 2018. The period corresponds to some period before sowing, the whole crop (i.e., winter wheat) cycle, and some period after harvesting the crop. The reason behind this is to make sure that the whole crop cycle is tracked on one hand, and to see whether there is any contribution of a previous and/or following crop, on the other. The pre-processing of L1C (Top of Atmosphere or TOA reflectance) Sentinel-2 images, which includes ortho-rectification, cloud removal (using cloud mask produced by Sen2Cor/SNAP), radiometric calibration, and atmospheric correction, was produced using SNAP/Sentinel-2 toolbox. The output of the pre-processing corresponded to L2A (bottom of atmosphere or BOA reflectance). The reason behind using the optical images (Sentinel-2) is to calculate the Normalized Difference Vegetation Index (NDVI). It has been shown that the NDVI is positively correlated with the above ground green biomass, meaning that as leafy green biomass increases, the NDVI gets closer to 1 [58].

Optical Data
Fifty-eight optical images were acquired by the two twin satellites of Sentinel-2 (Sentinel-2A and Sentinel-2B). The acquired images cover a period from November 2017 to August 2018. The period corresponds to some period before sowing, the whole crop (i.e., winter wheat) cycle, and some period after harvesting the crop. The reason behind this is to make sure that the whole crop cycle is tracked on one hand, and to see whether there is any contribution of a previous and/or following crop, on the other. The pre-processing of L1C (Top of Atmosphere or TOA reflectance) Sentinel-2 images, which includes ortho-rectification, cloud removal (using cloud mask produced by Sen2Cor/SNAP), radiometric calibration, and atmospheric correction, was produced using SNAP/Sentinel-2 toolbox. The output of the pre-processing corresponded to L2A (bottom of atmosphere or BOA reflectance). The reason behind using the optical images (Sentinel-2) is to calculate the Normalized Difference Vegetation Index (NDVI). It has been shown that the NDVI is positively correlated with the above Remote Sens. 2019, 11, 2228 5 of 29 ground green biomass, meaning that as leafy green biomass increases, the NDVI gets closer to 1 [58]. Hence, it is very useful to compute the NDVI for the crop to understand the temporal behavior with respect to the crop growth.
Consequently, the NDVI was computed. The NDVI is capable in predicting plants' photosynthetic activity by being determined from the red (ρ RED ) and near infrared (ρ NIR ) reflectance values, producing an index ranging from −1 to +1 [59], as follows:

Sentinel-1 Data
One hundred Sentinel-1A and Sentinel-1B images acquired at the ascending overpass, were downloaded between November, 2017 and August, 2018. These images were acquired at two incidence angles ranges, over the Bekaa plain of Lebanon: 32-34 • (50 images at ascending overpass) and 43-45 • (50 images at ascending overpass). The interferometric wideswath (IW) is a predefined mode by the ESA Sentinel-1 observation strategy, providing dual polarization (VV and VH) imageries at a 10-m spatial resolution with a revisit time of 6 days [60][61][62]. In this study, all the images were generated from the high-resolution, Level-1 Ground Range Detected (GRD) product. The data were pre-processed using the Sentinel-1 Toolbox in SNAP (developed by ESA). The preprocessing procedure included thermal noise removal, radiometric and geometric calibration, and speckle filtering [63,64]. Through the radiometric calibration, the raw image digital values were converted into backscattering coefficients in linear unit, whereas through the geometric calibration, the images were ortho-rectified using a Digital Elevation Model (DEM) of the Shuttle Radar Topography Mission (SRTM) at 30 m spatial resolution. For each reference plot, the mean backscattering coefficient (σ 0 ) was calculated from each calibrated Sentinel-1 image by averaging all pixel values within that plot.
The expression of the Sentinel-1 (S1) temporal behavior for the study site was through VV (dB), VH (dB), and the ratio VV/VH (dB) at the two ranges of incidence angles.
As synthesized by several studies, there are different factors affecting the radar backscatter, as the C-band represents an integration of the backscatter from the ground attenuated by the canopy layer, volume scattering, and stem-ground interactions (usually negligible for wheat) [54,65,66].
In addition to SAR configurations (incidence angle and polarization), the signal is affected by soil parameters (soil moisture and surface roughness) [67,68] and the vegetation, which affect the signal by two pathways, the attenuation of ground backscatter and direct volume scattering. This vegetation contribution is basically driven by the structure, dry matter, and vegetation water content [69].

Gaussian Decomposition of SAR Temporal Profiles
The time-series temporal profiles of the radar backscattering coefficient (σ 0 ) were fitted with Gaussian functions. Initially, VV, VH, and VV/VH (dB) were normalized using the unity-based normalization in Equation (2), to adjust the values at a common scale (between 0 and 1): where Y is the normalized σ 0 -value of the radar backscattering coefficient (in VV, VH, and VV/VH), Y is the initial value of the time series, Y min is the minimum value of the time series, and Y max is the maximum value of the time series.
The normalized profiles (in VV, VH, and VV/VH) were next smoothed using a Gaussian filter and then modeled with the sum of Gaussian functions (up to three Gaussians): where, i is the order of Gaussian, a i is the amplitude of the Gaussian, b i is the Gaussian peak position in Julian days, and c i is the standard deviation of the Gaussian in Julian days. The least square method was applied to determine the values of a i , b i , c i separately for VV, VH, and VV/VH (corresponding to the lowest difference between the smoothed profile and the function sum of Gaussians). To operate the least square method, initial values of (a i ,b i ,c i ) are chosen according to a preliminary analysis of the σ 0 -profiles.
The number of Gaussian functions necessary to model the smoothed profiles is equal to the number of times that the derivative of the smoothed profile changes from positive to negative (known local maximums).

In Situ Observations (Reference Plots)
In the West Bekaa plain of Lebanon (Figure 1), field campaigns were conducted on 21 reference winter wheat plots. For each winter wheat plot, the dates of the three phenological phases and the harvesting were noted (Table 1). In addition, as harvesting was expected to be sooner in the North Bekaa, eight reference winter wheat plots were selected and their harvesting dates were also observed ( Table 1).

Meteorological Data
Rainfall products with 1-day revisit time are obtained from Global Precipitation Measurement (GPM) sensor, which is the most recent sensor providing precipitation readings of up to 30-min revisit Remote Sens. 2019, 11, 2228 7 of 29 time [70]. The GPM mission is used in this study to retrieve rainfall data at 0.1 • latitude/longitude (~10 km × 10 km) over each part (West and North Bekaa) of the study site ( Figure 1). As for the relative humidity and air temperature, daily data were obtained from two weather stations (Ammiq station for West Bekaa and Doures station for North Bekaa), each is located in one of the two parts studied. The air relative humidity and temperature data were obtained for the three months April, May, and June of the corresponding season (2017-2018) for the sake of explaining the climatic difference among both sites affecting the date of the last phenological phase of wheat (i.e., physiological maturity).

Software Employed and Statistical Analysis
Since the proposed approach is operational, the method was performed through free open access software. For instance, the radiometric and geometric calibration of the Sentinel-1 images were performed using the Sentinel toolbox in SNAP. In addition, the Gaussian modelling of the SAR signal was performed using Python scripting 3.6.
Regarding the statistical analysis, the Root Mean Square Error (RMSE) and the difference between the estimated and observed dates (bias) were calculated for each phase as follows: where P i is the predicted date and O i is the observed date, for each reference plot (i) at each phenological phase, and the harvesting n is the total number of reference plots.

Methodological Approach
The methodological approach toward mapping the highly important three phenological phases in addition to the harvesting, throughout the cropping season 2017-2018, consists of four main steps ( Figure 2).
First, (1) calibration and terrain correction were performed on the downloaded S1 images to obtain the backscatter in both polarizations (VV and VH) in decibels (dB) in addition to the calculation of the ratio VV/VH (dB). Consequently, reference plots were selected and the key phenological dates were observed in situ (Section 2.3). Second, (2) together with the backscattering data, temporal profiles corresponding to the wheat reference plots were constructed to determine the best S1 configuration (optimum polarization and incidence angle), allowing the detection of each of the phases. After that, using Python scripting, smoothing of the radar backscatter and Gaussian fitting of the temporal series were applied on the reference plots in each polarization (VV and VH) and for the ratio VV/VH at both, low and high incidence angles (32-34 • and 43-45 • ), to come up with the most appropriate configuration for mapping each of the phases, and the harvesting.
Third, (3) on the already classified winter wheat plots in the West Bekaa plain of Lebanon, following a proposed classification approach by Nasrallah et al. [24], smoothing and Gaussian fitting were applied on all the winter wheat classified plots (2017-2018 season). Eventually, in step four (4), after the analysis of the temporal profiles of the reference plots, each phenological phase in addition to harvesting (in West Bekaa) was mapped using the most appropriate configuration (polarization and incidence angle). In Section 3.3, the results of the Gaussian fitting are shown, and each of the phenological phases in addition to harvesting are displayed with their corresponding utilized S1 configurations.
Eventually, to prove the robustness of the approach, since wheat is harvested earlier in the northern part of the plain than its western part, because of weather conditions and the wheat varieties, the same approach was performed on the northern part to estimate and map the harvesting time, as in West Bekaa. The behavior of S1 temporal profiles in VV, VH, and VV/VH at both incidence angles (32°-34° and 43°-45°) were analyzed using the reference plots located in the West Bekaa plain of Lebanon for which each winter wheat phenological phase's date was observed, in addition to the harvesting date. This was done to know the most appropriate S1 configuration in terms of polarization and incidence angle, in order to map each of the phenological phases, in addition to harvesting.
The NDVI temporal behavior was first interpreted, then compared with the analysis of the S1 temporal profiles to show the need of using the S1 data to map the important phenological phases and the harvesting dates of wheat. Then, the fitting results were illustrated for the S1 data, with the identification of the Gaussians' positions that were needed to perform the mapping. Afterward, since the West and North Bekaa plain have different harvesting time, weather data (air temperature and relative humidity) corresponding to April, May, and June of 2018 were used in order to analyze this different timing in achieving harvesting, between the two parts of the plain. Eventually, the maps of The behavior of S1 temporal profiles in VV, VH, and VV/VH at both incidence angles (32-34 • and 43-45 • ) were analyzed using the reference plots located in the West Bekaa plain of Lebanon for which each winter wheat phenological phase's date was observed, in addition to the harvesting date. This was done to know the most appropriate S1 configuration in terms of polarization and incidence angle, in order to map each of the phenological phases, in addition to harvesting.
The NDVI temporal behavior was first interpreted, then compared with the analysis of the S1 temporal profiles to show the need of using the S1 data to map the important phenological phases and the harvesting dates of wheat. Then, the fitting results were illustrated for the S1 data, with the identification of the Gaussians' positions that were needed to perform the mapping. Afterward, since the West and North Bekaa plain have different harvesting time, weather data (air temperature and relative humidity) corresponding to April, May, and June of 2018 were used in order to analyze this different timing in achieving harvesting, between the two parts of the plain. Eventually, the maps of the three phenological phase dates (i.e., germination, heading, and soft dough) plus the harvesting date were generated for West Bekaa in addition to the harvesting date map of North Bekaa.

NDVI Temporal Profiles
In this section, the NDVI temporal behavior of winter wheat within each part of the Bekaa plain (north and west) are presented ( Figure 3). The profiles include on each date, the average NDVI ± the standard deviation of the reference plots. Through the field work, the dates of sowing and harvesting, in addition to the three phenological phases (i.e., germination, heading, and soft dough) that are the focus of this study, are observed and represented on the graphs with vertical lines.
Remote Sens. 2019, 11, x FOR PEER REVIEW 9 of 30 the three phenological phase dates (i.e., germination, heading, and soft dough) plus the harvesting date were generated for West Bekaa in addition to the harvesting date map of North Bekaa.

NDVI Temporal Profiles
In this section, the NDVI temporal behavior of winter wheat within each part of the Bekaa plain (north and west) are presented ( Figure 3). The profiles include on each date, the average NDVI ± the standard deviation of the reference plots. Through the field work, the dates of sowing and harvesting, in addition to the three phenological phases (i.e., germination, heading, and soft dough) that are the focus of this study, are observed and represented on the graphs with vertical lines. On the temporal profiles, the daily precipitation is shown (right y-axis). On each profile, the three main phenological phases are denoted each by a vertical line, in addition to sowing and harvesting dates.
In the area of interest, sowing (red vertical bar) had taken place during the third week of November, after some consecutive rain events, ensuring a sufficient soil moisture for a good germination. After germination (pink vertical line), winter wheat goes through the vegetative cycle. The vegetative cycle includes seedling, tillering, stem elongation (jointing), booting, and heading (green vertical bar). After heading, flowering immediately occurs to kick off the reproductive cycle of winter wheat. After flowering, milky stage follows, followed by soft dough (blue vertical bar), hard In the area of interest, sowing (red vertical bar) had taken place during the third week of November, after some consecutive rain events, ensuring a sufficient soil moisture for a good germination. After germination (pink vertical line), winter wheat goes through the vegetative cycle. The vegetative cycle includes seedling, tillering, stem elongation (jointing), booting, and heading (green vertical bar). After heading, flowering immediately occurs to kick off the reproductive cycle of winter wheat.
After flowering, milky stage follows, followed by soft dough (blue vertical bar), hard dough and then ripening. When the crop reaches physiological maturity with no more than 15% of grain moisture, harvesting can take place (yellow vertical bar).
The NDVI has shown maximum values during April, 105 to 145 days after sowing (DAS). Except for heading, at which the NDVI values were at saturation level for a period of one month (Figure 3a), there was no clear significant indicators for the rest of the phases, and the harvesting, appearing in the temporal profile. At germination date, the NDVI had already started increasing, the same was observed for the soft dough phase. At harvesting, the NDVI values were already at minimum. Similar results were pinpointed in a recent study [55].

Sentinel-1 Temporal Profiles
In this section, the S1 polarizations (VV, VH, VV/VH) at a given incidence angle are illustrated to eventually allow the detection of each phenological phase and harvesting. After the analysis of the available configurations, the optimal configurations to map each phenological phase and harvesting are illustrated afterward.
The temporal profiles of VV, VH, and VV/VH (dB) at the two incidence angles (32-34 • and 43-45 • ) are shown in Figure 4. dough and then ripening. When the crop reaches physiological maturity with no more than 15% of grain moisture, harvesting can take place (yellow vertical bar). The NDVI has shown maximum values during April, 105 to 145 days after sowing (DAS). Except for heading, at which the NDVI values were at saturation level for a period of one month (Figure 3a), there was no clear significant indicators for the rest of the phases, and the harvesting, appearing in the temporal profile. At germination date, the NDVI had already started increasing, the same was observed for the soft dough phase. At harvesting, the NDVI values were already at minimum. Similar results were pinpointed in a recent study [55].

Sentinel-1 Temporal Profiles
In this section, the S1 polarizations (VV, VH, VV/VH) at a given incidence angle are illustrated to eventually allow the detection of each phenological phase and harvesting. After the analysis of the available configurations, the optimal configurations to map each phenological phase and harvesting are illustrated afterward.
The temporal profiles of VV, VH, and VV/VH (dB) at the two incidence angles (32°-34° and 43°-45°) are shown in Figure 4.  In Table 2, descriptive statistics in terms of NDVI and SAR backscatter are provided for each of the phenological phases. The statistics are thus provided for West Bekaa including the mean and the standard deviation of wheat for NDVI and SAR configurations (polarizations and incidence angles).  In West Bekaa, winter wheat was sown during the third week of November, 2017. By looking at the VV/VH (dB) (Figure 4c), germination (pink vertical line) was observed three weeks after sowing. This phase was well observed as the ratio VV/VH (dB) was maximum. This means that at germination, VV/VH (dB) showed a maximum value before starting decreasing. This decrease was sharper at a low incidence angle (32-34 • ) than at the higher one. For this reason, VV/VH at 32-34 • was selected for mapping the germination.
By looking at the VV polarization curve (Figure 4a), the signal backscatter fluctuated between −12.8 dB and −6.2 dB until the start of March, 2018 (DAS 90) because of continued events of rain [1]. The signal was basically affected by the soil contribution [71] as until then, the canopy height reached a maximum of 35 cm at DAS 90. After the jointing phase started, the canopies witnessed a growing boot at which they began elongating when the nodes started developing above the crown node. During stem elongation, the S1 signal in VV polarization decreased as a result of soil scattering attenuation caused by vertical stems and leaves [72]. Heading (green vertical line) started in West Bekaa at 144 DAS. At heading, the spike (ear) started emerging out from the leaf sheath after the upper-most leaf swelled out into flag holding the spike into it. At the heading phase (DAS between 144 and 150), where the vegetative attenuation was at its extreme [73,74], the VV backscatter was between −15.1 dB and −13.1 dB (lowest level), leading to the best observation of this phase (i.e., heading) in VV polarization at low incidence angle (32-34 • ) (Figure 4a). Just after heading, anthesis of florets (flowering) and fertilization of ovaries occurred.
After flowering (post-anthesis phase) and when pollination was complete (DAS between 155 and 180), the ovaries started their transformation into seeds (grain filling). The process started by increasing moisture levels into the spikelet, passing through milky phase and then soft dough phase (blue vertical bar). At soft dough phase, the average VH backscatter at high incidence angle (43-45 • ) showed a maximum backscatter of −15.9 dB at 186 DAS (Figure 4e). Similar to high incidence angle (43-45 • ), at low incidence angle (32-34 • ) (Figure 4b), the signal also showed an increase from heading to soft dough phase with a similar maximum signal of −16.0 dB. However, the increase from heading to soft dough was sharper at high incidence angle (43-45 • ) (ranging from −20 dB to −15.9 dB) than at low incidence angle (ranging from−18.4 dB to −16.0 dB). This shaper increase at 43-45 • was seen because at heading, the VH reached a lower value at the high incidence angle, compared to that at low incidence angle. Thus, since the increase from heading to soft dough phase at high incidence angle was more significant (sharper) than that at low incidence angle, high incidence angle images, when available, are more recommended to be used for estimating the date of this phase (i.e., soft dough).

Optimal S1 Configuration for Mapping Harvesting (West and North Bekaa)
After the soft dough phase was reached, the grains moisture level started decreasing reaching the hard dough phase and then maturation. At maturation phase, the kernels became hard and moisture percentage was gradually reduced and the plant could be harvested when physiological maturity was reached (yellow vertical bar). This event was clearly shown by referring to the VV/VH ratio, as similar to germination phase, the ratio VV/VH reached maximum again (Figure 4c). The VV/VH ratio reached a maximum value (DAS 228) because of the sharp decrease in the VH polarization because of the absence of volume and multiple scattering after harvesting event.
The harvesting of winter wheat was during the third week of July, 2018 for the plots in West Bekaa, while in the first 10 days of July for those in the North Bekaa, this was mainly due to weather-related conditions ( Figure 5). As the reproductive phase started (after DAS 150), the canopies dried down and canopy water content was reduced. As the temperature was higher and the air relative humidity was lower, the rate of canopy water content reduction was higher, leading to an early harvest [75]. Regarding the ratio VV/VH, throughout the whole season, the effect of rain was reduced for the two regions and therefore the signal fluctuations because of rain were not observed as for VV and VH alone. The VV/VH (dB) profile decreased (Figure 4c) until heading. After heading, the VV/VH increased until harvesting took place. As the sharpest increase from soft dough to harvesting was seen in VV/VH (dB) at a low incidence angle (32-34 • ) (from around 5.8 until 9 dB, Figure 4c), while at high incidence angle (43-45 • ), the increase seen from soft dough to harvesting was from around 6 to 7 dB. Thus, VV/VH (in dB) at low incidence angle was chosen to map harvesting in both parts (West and North Bekaa plain).
By referring to Figure 5, the difference in air temperature (Ta) and air relative humidity (RH) among both parts was seen over April, May, and June of the corresponding year (2018). In North Bekaa, the monthly averages RH for the three months were 48.7, 50.6, and 45. By referring to Figure 5, the difference in air temperature (Ta) and air relative humidity (RH) among both parts was seen over April, May, and June of the corresponding year (2018). In North Bekaa, the monthly averages RH for the three months were 48.7, 50.6, and 45.1, respectively. As for the west part, the RH was higher, showing values of 58.2, 58.5, and 54.1 for the same three months, respectively. Regarding the Ta, both parts witnessed similar average temperature in April of 15 °C, but slightly lower Ta in the west part for May and June (19°C and 21 °C) than the Northern part (19.3 °C and 22 °C), for the same months.

Smoothing and Gaussian Fitting
The automated method for phenology phase dates estimation was applied on S1 data. After performing the smoothing and Gaussian fitting on S1 data, the fitting output parameters were obtained automatically and used to map the three phenological phase dates as well as the harvesting date. Figures 6, 7, and 8 illustrate an example of the smoothed profile and the Gaussian fitting used

Smoothing and Gaussian Fitting
The automated method for phenology phase dates estimation was applied on S1 data. After performing the smoothing and Gaussian fitting on S1 data, the fitting output parameters were obtained automatically and used to map the three phenological phase dates as well as the harvesting date. Figures 6-8 illustrate an example of the smoothed profile and the Gaussian fitting used to model the S1 data to estimate the dates of each of the phenological phases. Figure 6 corresponds to the output of the smoothing and Gaussian fitting of the VV/VH at low incidence angle (32-34 • ) for estimating germination and harvesting dates, Figure 7 corresponds to the output of the smoothing and Gaussian fitting of the VV polarization at low incidence angle for estimating the heading date, and Figure 8 corresponds to the VH polarization smoothing and Gaussian fitting to estimate the soft dough date, at a high incidence angle (43-45 • ). Table 3 explains the way in which the phenological phases (germination, heading, soft dough, and harvesting) were automatically detected from the Gaussian fitting. The detection of each phase is accompanied with the best S1 configuration (i.e., incidence angle and polarization). For each phase, the best SAR configurations already determined through the analysis of S1 temporal variations (Section 3.2), are considered.
In practice, after the first peak in the sum of Gaussians fitting (positive derivative) was found (using VV/VH) and the germination dates for each of the plots were assigned, the next step was to find the heading dates, which is the phase that follows. On the VV Gaussian fitting of all the wheat plots, the first negative derivative peak in the sum of the Gaussians fitting, which comes after the last recorded date of germination, was assigned as the heading date, for each plot. On the VH Gaussian fitting of all the wheat plots, the first positive derivative peak in the sum of Gaussians fitting, which is located after the last recorded heading date, was assigned as a soft dough date, for each of the plots. Eventually, relying on the VV/VH, the harvesting date was found as the last positive derivative peak in the sum of Gaussians fitting. The harvesting date was determined after the last recorded date for soft dough.
to model the S1 data to estimate the dates of each of the phenological phases. Figure 6 corresponds to the output of the smoothing and Gaussian fitting of the VV/VH at low incidence angle (32°-34°) for estimating germination and harvesting dates, Figure 7 corresponds to the output of the smoothing and Gaussian fitting of the VV polarization at low incidence angle for estimating the heading date, and Figure 8 corresponds to the VH polarization smoothing and Gaussian fitting to estimate the soft dough date, at a high incidence angle (43°-45°).

Figure 6.
Gaussian fitting of the S1-VV/VH temporal series of a winter wheat plot at the West Bekaa plain of Lebanon. The incidence angle is 32°-34°. The x-axis represents the date and the y-axis refers to the normalized signal. The black arrows indicate the inflection points where germination and harvesting had taken place.   Gaussian fitting of the S1-VV temporal series of a winter wheat plot at the West Bekaa plain of Lebanon. The incidence angle is 32°-34°. The x-axis represents the date and the y-axis refers to the normalized signal. The black arrow indicates the inflection point where heading had taken place.   Table 3. The best S1 configuration (polarization and incidence angle) for phenological phases (i.e., germination, heading and soft dough) and harvesting mapping in addition to the way of determination for each.  Figure 9 shows for each wheat plot the date of each phenological phase for West Bekaa (WB) (maps corresponding to germination, heading, soft dough, and harvesting). For North Bekaa (NB), the produced map corresponds to harvesting only as for the other phases the timing is similar as the West Bekaa. Each phenological phase stretched up to up to few weeks. For this reason, in each generated map (Figure 9), a color bar was added in which a unique color code is assigned for each 6 days. It is important to note that the whole time period for each phenological phase, in addition to harvesting, does not mean that the winter wheat in the Bekaa plain needs this relatively long time (e.g., 2 months for harvesting in North Bekaa) to achieve a particular phase. The color code bars on the maps reflect the first estimated date until the last estimated date. In the following section, the density of each 6-day period (percentage of wheat plots per each 6 days) is demonstrated.

Accuracy Assessment and Quantitative Analysis
After detection each of the three phenological phases for West Bekaa and harvesting event for both West and North Bekaa from S1 Gaussians fitting, the obtained (estimated) dates of each phenological phase and the in situ observations (observed dates) were compared in order to evaluate the precision of the phenological phases detecting approach, using S1 data. The RMSE and the difference between the estimated and observed dates (bias) were calculated for each phase ( Figure  10). It is important to note that the whole time period for each phenological phase, in addition to harvesting, does not mean that the winter wheat in the Bekaa plain needs this relatively long time (e.g., 2 months for harvesting in North Bekaa) to achieve a particular phase. The color code bars on the maps reflect the first estimated date until the last estimated date. In the following section, the density of each 6-day period (percentage of wheat plots per each 6 days) is demonstrated.

Accuracy Assessment and Quantitative Analysis
After detection each of the three phenological phases for West Bekaa and harvesting event for both West and North Bekaa from S1 Gaussians fitting, the obtained (estimated) dates of each phenological phase and the in situ observations (observed dates) were compared in order to evaluate the precision of the phenological phases detecting approach, using S1 data. The RMSE and the difference between the estimated and observed dates (bias) were calculated for each phase (Figure 10).
By analyzing the RMSE values, the germination and harvesting for West Bekaa have shown the least RMSE with 2.9 and 3.0 days, respectively. Heading and soft dough phases, for the same part of the plain, have shown higher RMSE of 5.5 and 5.1 days, respectively. As for the North Bekaa, 4.5 days of RMSE was observed for harvesting event. In addition, a slight negative bias corresponding to an underestimation was observed for germination and heading with −0.2 and −1.1 days, respectively. As for the soft dough and harvesting for West and North Bekaa, an overestimation was observed with 3.1, 0.6, and 3.6 days, respectively. Table 4 shows the percentage of the winter wheat plots within each 6-day period. For better comparison regarding the harvesting of the West and the North Bekaa, same periods (6-day period) are assigned for both.
According to Table 4, during germination, around 79% of the plots had completed germination between 3 December and 2 January 2018, which is very logical due to the already known period of sowing, which was during the third week of November, of the same year. As for heading phase, around 73% of the plots had completed heading between 2 and 26 April 2018. For soft dough phase, around 88% of the plots had achieved this phase between 16 May and 9 June 2018. As for harvesting, the difference in timing between the two parts of the Bekaa plain was observed. The average date of harvesting in West Bekaa is 20 July 2018, while for the North Bekaa, 13 July 2018. In addition, 85% of the wheat in the North Bekaa was harvested between 2 and 26 July 2018, while for the West Bekaa, 92% of wheat was harvested between 8 July and 1 August 2018. These results are very typical for the Bekaa plain and explain the difference in harvesting time among the West and North Bekaa, because of the weather conditions as explained in Section 3.2.2.

Toward Near-Real Time Phenology Monitoring
The proposed method is built on the basis of understanding the crop's phenology by modelling the temporal series of the crop (winter wheat for our case) with Gaussian functions to estimate the phenological phases' dates. However, for real time and near-real time operational phenology monitoring, we will assume in this section, for each phenological phase (in West Bekaa), that the whole temporal series was not available. This is to test with the available Sentinel-1 images and reference plots the error of estimation (RMSE) when the temporal series is considered until the date of the phenological phase, 6, 12, and 18 days after the phenological phase is achieved. In Table 5, we demonstrate for each phenological phase the new difference between estimated and observed dates (bias) and RMSE together with the percentage of plots for which the dates of each phenological phase could be estimated. Table 5. The difference between the estimated and observed dates (bias), the RMSE of the estimated phenological date, and the percent (%) of reference plots for which the date of the different phenological phases could be estimated when the time series ended on the date of the phenological phase, 6, 12, and 18 days later.

S1 Versus NDVI Temporal Behavior
The winter wheat Sentinel-2 (optical) and Sentinel-1 (SAR C-band) remote sensing temporal profiles are shown in Figures 3 and 4. One of the most noticeable aspect is the high sensitivity of SAR temporal behavior to winter wheat growth cycle, in comparison to the NDVI temporal series [55]. The crop seasonal cycle is broken down into five periods and each is investigated and discussed.
From day 0 to day 90 after sowing, the NDVI ( Figure 3) started increasing sharply just after the crop emergence, as it is therefore correlated to the greenness of the canopies [76][77][78]. Consequently, the VV and VH polarizations temporal profiles witnessed a fluctuation over this period. This is because the S1 backscatter is highly sensitive to soil moisture, as in this period continued rainfall events were recorded and thus explaining that both VV and VH polarizations are affected mostly by variations in the soil backscatter driven by soil water content (SWC) only, as surface roughness does not change significantly from sowing to harvesting date [65,68,69]. The VV/VH temporal profiles (Figure 4) reduce the dependence of the radar signal on soil moisture especially after a rain event. Indeed, the effect of soil moisture is not completely eliminated with VV/VH since the sensitivity of the radar signal in the C-band to soil water content is slightly different in VV and in VH polarizations [79]. In addition, VV/VH showed a mild decreasing trend, from after germination until heading. This decrease of VV/VH with time, was mainly due to the reduction of the difference between the signal radar at VV and that at VH, caused by increased vegetation volume [80,81].
From DAS 90 to DAS 105, the NDVI (Figure 3) profile kept on increasing steadily as the photosynthetic activity, LAI and canopies stems number and their height continuously increased in this period. In regard to the radar backscatter, fluctuations in both VV and VH polarizations (Figure 4) faded away and a decrease in the backscatter was observed for two main reasons. First, soil moisture decreased because of no rain events that led to weaker soil contribution, and second, in this period, wheat canopies were tillering and jointing, leading to increased LAI and thus causing more soil backscattering attenuation [82].
From DAS 105 until heading was achieved (DAS 138-144), a maximum saturation in the NDVI (Figure 3) was observed, which is an already known limitation that has been already described previously [42,69]. This is attributed to the highly dense vegetation (maximum LAI). Per contra, both VV and VH polarizations (Figure 4) kept on their continuous sharp decrease to reach the minimum observed backscatter. In this period (DAS 105 until heading), VV decreased faster and sharper than the VH. At the heading period (appearance of ears), the direct vegetation scattering was low and the attenuation of the soil contribution was high. This is explained by the ratio (VV/VH) as the minimum value was reached at heading time ( Figure 4).
From 144 DAS (heading) till 186 DAS (soft dough), the NDVI sharply decreased as the canopies dried up and the photosynthetic activity was over. The green color started turning into yellowish and retrenchment of leaves and stems was observed. While on the other hand, in the VV and VH polarizations (Figure 4), an increase in the backscatter was seen to reach a maximum point. Regardless the incidence angle and the polarization, the main mechanisms taking place in the backscatter are attributed to leaves and ears backscatter (volume scattering) which is in agreement with previous studies [82][83][84]. In addition, as demonstrated in the previous section, in this period, grain filling process was taking place and moisture level was increasing in the wheat ears during the post-pollination stage. The canopy was well developed and the soil contribution was very low (also the stem-ground interaction), meaning that the VV and VH backscatter thus could be also increased by the rising ears moisture level and not by soil contribution, contradicting earlier study [84]. Consequently, the ratio (VV/VH) was constant with a slight increase as VH backscatter was more attenuated by the vegetation [55,84].
From soft dough to harvesting, the sharp decrease in the NDVI turned into a softer one as the rate of "greenness" fading away had sharply decreased and physiological maturity was being achieved. As for the SAR behavior, a soft decrease was seen in the VV polarization, because of decreasing moisture levels in the canopies heads and the moisture within the grains had dried up. As for the VH, which is more affected by vegetation contribution, the decrease was sharper. Consequently, the ratio (VV/VH) continued on increasing as fresh biomass was decreasing until physiological maturity (harvesting time).

Influence of S1 Incidence Angle
Through the examination of the different incidence angles, three phases are discussed, namely germination, heading, and soft dough in addition to physiological maturity (harvesting), in both VV and VH polarizations and the ratio VV/VH: (1) In VV polarization: From the start of jointing till heading (84 DAS until 144 DAS), the decrease in the signal at 32-34 • incidence angle was steeper and sharper than at the higher incidence angle (43-45 • ) because at high incidence angle, in addition to the attenuation, there is the smaller direct vegetation contribution [85]. In addition, at 40 • of incidence angle and beyond, the direct vegetation volume scattering appears [84]. From 144 DAS until 186 DAS (soft dough), the signal increased at the two incidence angles. However, at soft dough phase (186 DAS), different behaviors were significantly recorded among the two ranges of incidence angles. On this date (186 DAS), the backscatter at 43-45 • was slightly higher than that at 32-34 • (Figure 4a,d) by around 1.5 dB, meaning that at the soft dough phase, at the 43-45 • incidence angle, the signal held more canopy contribution than the radar signal did at lower incidence angle (32-34 • ). Such a behavior can be explained by the fact that after heading had occurred, the soil contribution, which was a dominant backscattering mechanism was considerably reduced and the canopy backscatter became more significant (at 43-45 • ). This finding is noted by different previous studies [82][83][84][85]. At harvest, the two incidence angles showed similar σ • . (2) In VH polarization: As wheat canopies reached heading phase, the S1 backscatter at high incidence angle reached lower levels than at low incidence angle (Figure 4b,e). This showed a sharper increase of the S1 signal from heading to soft dough at high incidence angle than at low incidence angle. This is the reason why mapping this phase (soft dough) using high incidence angle (43-45 • ) was more appropriate. Hence, as stated before (Section 3.3), the VH polarization has been seen as a better configuration through the analysis of the S1 temporal profiles for estimating the soft dough date. However, for easier operational application, VH polarization at low incidence angle (32-34 • ) could still be used. Nevertheless, when VH at 32-34 • was used to map the soft dough phase, the soft dough could not be detected for around 10% of the wheat plots. In addition, the detection of the soft dough phase using 32-34 • showed that for about 18% of the wheat plots, a different soft dough date estimation was observed of at least 6 days, in comparison to the one estimated at high incidence angle. Figure 11 shows an example to detect the soft dough phase for the same wheat plot, one at 32-34 • and one at 43-45 • incidence angle. It highlights the difference in the fitting (Gaussian) of a VH temporal profile to detect the soft dough phase.
at low incidence angle. This is the reason why mapping this phase (soft dough) using high incidence angle (43°-45°) was more appropriate. Hence, as stated before (Section 3.3), the VH polarization has been seen as a better configuration through the analysis of the S1 temporal profiles for estimating the soft dough date. However, for easier operational application, VH polarization at low incidence angle (32°-34°) could still be used. Nevertheless, when VH at 32°-34° was used to map the soft dough phase, the soft dough could not be detected for around 10% of the wheat plots. In addition, the detection of the soft dough phase using 32°-34° showed that for about 18% of the wheat plots, a different soft dough date estimation was observed of at least 6 days, in comparison to the one estimated at high incidence angle. Figure 11 shows an example to detect the soft dough phase for the same wheat plot, one at 32°-34° and one at 43°-45° incidence angle. It highlights the difference in the fitting (Gaussian) of a VH temporal profile to detect the soft dough phase. In the VH polarization at low incidence angle (Figure 11a), the soft dough phase could not be detected through the Gaussian fitting, unlike Figure 11b. As discussed before, through the analysis of the temporal profiles (Figure 4), the increase from heading to soft dough was sharper at higher incidence angle (43°-45°). Indeed, for some wheat plots, the fitting process did not fit a Gaussian at the level of soft dough because the increase from heading to soft dough was not significant enough.
3) In VV/VH ratio: The steady mild decrease from sowing to heading in the ratio VV/VH (dB), which was seen at the two incidence angles (32°-34° and 43°-45), is mainly related to the slight increase in the VH (mainly from sowing till the beginning of March). The VH backscatter is dominated by the volume scattering mechanisms, increased as reported previously [55,80,81] while however, the VV backscatter, which is dominated by the direct contribution from the ground and the canopy decreased because of the rising attenuation In the VH polarization at low incidence angle (Figure 11a), the soft dough phase could not be detected through the Gaussian fitting, unlike Figure 11b. As discussed before, through the analysis of the temporal profiles (Figure 4), the increase from heading to soft dough was sharper at higher incidence angle (43-45 • ). Indeed, for some wheat plots, the fitting process did not fit a Gaussian at the level of soft dough because the increase from heading to soft dough was not significant enough.
(3) In VV/VH ratio: The steady mild decrease from sowing to heading in the ratio VV/VH (dB), which was seen at the two incidence angles (32-34 • and 43-45), is mainly related to the slight increase in the VH (mainly from sowing till the beginning of March). The VH backscatter is dominated by the volume scattering mechanisms, increased as reported previously [55,80,81] while however, the VV backscatter, which is dominated by the direct contribution from the ground and the canopy decreased because of the rising attenuation from the predominantly vertical structure of the wheat stems [82], especially from March (96 DAS) through heading, during the stem elongation. From heading to soft dough phase, the ratio (VV/VH) was more or less constant at the two incidence angles (32-34 • and 43-45 • ) as the signal in both VV and VH equally increased throughout this period as explained before (Section 3.3). From soft dough to harvesting, VV/VH (dB) at low incidence angle had shown a more significant dynamics than at high incidence angle. Thus, harvesting was observed when VV/VH (dB) at low incidence angle increased to reach the maximum.

Wheat Phenology Mapping
The germination map contains ten 6-day periods of wheat plots' germination dates. The different timing of germination was strongly dependent on the time of sowing. Almost 90% of the plots in both West and North Bekaa achieved germination throughout December, 2017 (Table 4). This result is reasonable and logical as wheat plots were mainly sown between mid-November and the first week of December of 2017, knowing that wheat seeds start their germination 10 days after sowing [86], up to one month depending on the weather, soil and variety. As for heading time (Figure 9), the estimated dates showed that around 73% of the wheat plots (Table 4)  2018. In the area of interest (the Bekaa plain of Lebanon), this period is a typical heading period [24]. The third mapped phenological phase after heading was the soft dough phase. By looking at Figure 9 and Table 4, around 88% of wheat plots achieved soft dough between 16 May and 9 June 2018. As for harvesting event, in West Bekaa (WB), harvesting took place between 2nd and the 3rd weeks of July, 2018. As for the North Bekaa, the same event took place around one week earlier (Table 4). This difference is explained by weather and wheat varieties. According to Figure 5, the higher air relative humidity in the West Bekaa compared to North Bekaa could be responsible for the late harvesting of wheat. In fact, at harvesting, it is well recommended to reach a grain moisture content of 11-14% [17], a level that is achieved faster in areas with lower air relative humidity as North Bekaa. In addition to that, it was also seen from Figure 5 that the temperature over the last few months of the season was slightly higher for North Bekaa, which could also help in faster grain moisture level reduction. Apart from weather conditions, wheat varieties do affect the length of the season. According to the survey conducted at the plain, farmers in the north part relied on the landraces distributed by the Lebanese Agricultural Research Institute (LARI), which have a relatively shorter season than those varieties cultivated in the west part (e.g., Saragolla and Normano).
Our results have shown that indeed, winter wheat phenological phases are directly interrelated. Those plots which have witnessed early germination, achieved heading early (10 April 2018), while those which have witnessed a late germination, witnessed a late heading (20 April 2018). Plots that achieved early heading reached soft dough phases earlier than those which achieved late heading (25 May and 1 June, 2018, respectively). Eventually, winter wheat plots that finished soft dough early were harvested earlier than those finished soft dough late (14 July and 25 July, 2018, respectively).

Quality Indicator, Strengths, and Perspectives
The quality indicators used in this study are the Root Mean Square Error (RMSE) and the difference between the estimated and the observed dates (bias) for each phenological phase. The RMSE and the bias results were satisfactory (less than 6 days for all the phenological phases). Each of germination, heading and soft dough lasted for up to few days, thus achieving an estimated date with an RMSE of maximum 5 days was satisfactory for the needed interventions. As for harvesting, farmers take few days for storing their grain yield and moving them to warehouses. In addition, the maximum RMSE was less than the Sentinel-1 revisit time (6 days). As for the bias, a slight underestimation (negative bias) was observed for germination and heading (−0.2 and −1.1 days, respectively). This is basically because the field campaigns could not be scheduled every day, leading to a relatively late observation for some plots of an already achieved phase. As for the soft dough and harvesting for West and North Bekaa, the bias was positive (3.1, 0.6, and 3.6, respectively), thus an overestimation was noted. For soft dough, the reason could be that because soft dough phase normally lasts for some days before hard dough starts, the date was estimated by our approach when the dough-like grain was more dry than when observed in situ. As for harvesting, the dates could be overestimated because of the date reporting-error by farmers, as some of them were contacted to ask for harvesting date. However, recording the real dates on the field is challenging and leads to bias when comparing it with the estimated dates. Wheat plots, especially big ones, do not have a unified time in which canopies within the plot reach the same phenological phase at one time, which highly depends on the plot slope, sunlight, water availability, and agricultural practices. Still, it is important to note that organizing campaigns to the study area at the perfect timing could be complicated, thus, the recording of a particular phase could be late by few days.
The proposed automated phenology mapping approach has been proven to be robust. Since in the North Bekaa the wheat was harvested earlier, the approach was automatically applied on the North Bekaa and the results were satisfactory by looking at the low RMSE. In addition, the intra-relationships between the phases along the phenology as a whole were well seen and strongly related, since the timing of each phenological phase affected the timing of the one that follows (Section 3.4.1).
It is clear that when fitting the whole (full) temporal series with Gaussian functions, the RMSE is significantly lower than when fitting until the date (or a close date) of the phenological phase. However the RMSE of the phenological phases and harvesting when fitted until a date close to the phenological phase is promising. In addition, the bias (difference between estimated and observed dates) is significantly higher when the temporal series is fitted with Gaussian functions to estimate the dates of the phenological phases. Both heading and soft dough phases' dates were underestimated unlike germination and harvesting (overestimated).
Actually, it is a matter of the type of operational use and how much one is willing to compensate for the error for having a near-real time monitoring. In addition, for future work, putting extra efforts on this matter is recommended, by having a higher number of reference plots for a better insight on the achievable RMSE for the sake of real time and near-real time phenology monitoring.
Wheat is a very important crop and exists in different species (e.g., durum, spelt, emmer, einkorn . . . etc.). Since the radar signal depends mainly on the crop morphology, different species sharing the same season (having a similar morphology as the studied case) have a strong potential to be monitored phenologically, following the proposed approach. However, for the wheat species that are cultivated over a different season (e.g., spring wheat), further investigation will be carried out. Since the radar signal is basically sensitive to soil moisture (when the NDVI is less than 0.7), the rainfall should be closely observed in parallel with observing the dates of each of the phenological phases. In these cases, the dependency on the NDVI might increase when precipitation takes place

Conclusions
Mapping winter wheat phenology using S1 data was proposed in this study. Using SAR S1 (C-band data) with field observations allowed us to estimate and map the dates of three important phenological phases of wheat plus the harvesting event, because of their importance and relevance to decision making and farmers' interventions. The Root Mean Square Error (RMSE) calculated for each of the phases in West Bekaa revealed the values of 2.9, 5.5, 5.1, 3.0, and 4.5 days for germination, heading, soft dough, harvesting for West Bekaa, and harvesting for North Bekaa, respectively. In addition, a slight underestimation was observed for germination and heading dates of West Bekaa (−0.2 and −1.1 days, respectively) while an overestimation was observed for soft dough of West Bekaa and harvesting for both West and North Bekaa (3.1, 0.6 and 3.6 days, respectively). Moreover, we have also observed that early germination had led to early heading, early heading had led to early soft dough, early soft dough had led to early harvesting, and vice versa.
The three main strengths of the mapping approach proposed in this study are: (1) Automation through Python scripting that has led to fast results generation, (2) robustness when mapping harvesting on a different site with different harvesting dates, and (3) adaptation to the strong intra-relationships between winter wheat phenological phases (germination, heading, soft dough, and harvesting) as the date of achieving each phase has affected the following phases.
To test to what extent the approach can be reliable when a non-full temporal series was fitted, a sensitivity analysis was conducted. It was found that the RMSE and bias significantly increased when the temporal series were fitted until the phenological phase's date, 6, 12, and 18 days after. However, the values are still promising, yet more efforts are required on this matter.
In future work, extending the methodology to other important crops (e.g., potato, maize, other cereals and vineyards) on different sites in the Mediterranean area is the main target to be achieved, to be eventually merged with crop modelling systems, especially when on-field verifications are not easily implemented and/or when they are very costly.
Moreover, following our approach, if there will be a shift in the standard dates of the phenological phases (because of weather conditions), regional mapping of this shifting will be possible with reduced budgets for physical check. Nevertheless, coupling our automated approach with simulation crop models, especially those which simulate crops' growth depending on the phenology (e.g., CropSyst) could be very beneficial, and would nonetheless lead to decreased simulations' uncertainties. Funding: This work was financially supported by the grant research programme project provided by the Conseil National de la Recherche Scientifique (CNRS-Liban) and CIHEAM-IAMM and implemented in collaboration with the National Center for Remote Sensing (NCRS) and Irstea (Montpellier, France).