Using UAV-Based SOPC Derived LAI and SAFY Model for Biomass and Yield Estimation of Winter Wheat

: Knowledge of sub-ﬁeld yield potential is critical for guiding precision farming. The recently developed simulated observation of point cloud (SOPC) method can generate high spatial resolution winter wheat e ﬀ ective leaf area index (SOPC-LAIe) maps from the unmanned aerial vehicle (UAV)-based point cloud data without ground-based measurements. In this study, the SOPC-LAIe maps, for the ﬁrst time, were applied to the simple algorithm for yield estimation (SAFY) to generate the sub-ﬁeld biomass and yield maps. First, the dry aboveground biomass (DAM) measurements were used to determine the crop cultivar-speciﬁc parameters and simulated green leaf area index (LAI) in the SAFY model. Then, the SOPC-LAIe maps were converted to green LAI using a normalization approach. Finally, the multiple SOPC-LAIe maps were applied to the SAFY model to generate the ﬁnal DAM and yield maps. The root mean square error (RMSE) between the estimated and measured yield is 88 g / m2, and the relative root mean squire error (RRMSE) is 15.2%. The pixel-based DAM and yield map generated in this study revealed clearly the within-ﬁeld yield variation. This framework using the UAV-based SOPC-LAIe maps and SAFY model could be a simple and low-cost alternative for ﬁnal yield estimation at the sub-ﬁeld scale.


Introduction
Precision agriculture aims at optimizing input and output in field operations in order to achieve maximum economic profit while maintain environmental sustainability [1]. The sub-field level crop monitoring can provide finer spatial resolution (meter level) information and reveal the within-field crop variability. Information on the spatial variation of crop biomass and yield at the sub-field level is directly relevant to increasing farm profit by addressing the low-productivity areas within a field. Remote sensing has long been recognized as an effective means to provide multi-temporal information on crop growth over large areas in support of precision agriculture [2][3][4]. For example, moderate resolution imaging spectroradiometer (MODIS), Landsat, and RapidEye optical satellite data have been used to monitor crop growth status throughout the growing season using vegetation spectral indices and crop models [5][6][7][8]. Although the spatial and temporal resolution of satellite imagery has been improved over the years, it is still incapable of providing timely and detailed information of within-field variations for operational applications [9]. The recent advancement of the unmanned aerial vehicle (UAV) system has overcome the spatial and temporal limitation of satellite data for precision agriculture [10][11][12]. The high spatial and temporal UAV-based imagery can provide important information for monitoring the within-field variabilities of crop status during the growing area, and can automatically generate LAIe maps showing within field variation [47]. The LAIe is the result of the indirect approach retrieved LAI value using a non-destructive method. If the canopy satisfies the assumption of a random spatial distribution, the LAIe could be calculated from the gap fraction of the canopy [48]. The UAV-based photogrammetry point cloud data contains crop structural information. The SOPC derived LAIe maps from the UAV-based point cloud data can clearly indicate the winter wheat LAIe spatial variability without using ground based LAIe measurements. In addition, the SOPC method is not affected by shadow and illumination, which does not require a radiometric calibration for UAV-based imagery. Currently, Canada has no restriction on UAV system less than 25 kg in agricultural operating. This may encourage the end users to adopt the innovative technology in crop monitoring. The demand for crop monitoring will require a low-cost and accessible approach to achieve the estimation of crop DAM and yield for farmers. Due to the low-cost and easy operation of SOPC method, the SOPC derived LAIe has a great potential for final DAM and yield estimation.
The overall objective of this study aims at developing a simple and low-cost UAV-based approach for generating a high-resolution final DAM and yield maps without ground-based measurements. The SOPC derived UAV-based point cloud LAIe (SOPC-LAIe) were applied to the SAFY model, using winter wheat as an example, to generate the final DAM and yield map to represent the DAM and yield spatial variabilities. The study is designed: (1) to determine winter wheat CSPs from DAM in the SAFY calibration instead of LAI measurements, (2) to generate high spatial resolution multi-temporal LAI maps using the newly developed SOPC method on UAV-based point cloud data, and (3) to generate winter wheat final DAM and yield map using SAFY model calibrated with the UAV-based LAIe maps.

Study Area
The study site locates in the west of Melbourne, southwestern Ontario, Canada (42.787707 • N, 81.594801 • W) (Figure 1a). This region has productive soil and abundant water supply. Due to its long cold winter and shorter growing season (April to September), there is only one harvest per year for field crops. Winter wheat is one of the major crops grown in this region, which is typically sown in the previous fall and regrow in the following spring after snowmelt. The selected winter wheat field is 41  (b) shows the location of S1 and S2 in the wheat field. The blue points are the sampling locations in S1, and the red points are the sampling locations in S2.

Field Sampling Design and Field Data Collection
Since the UAV imagery was to be collected at very high spatial resolution, the corresponding destructive biomass collection will affect the subsequent crop parameter measurements. To circumvent this situation, two sub-fields, S1 and S2, were selected within this winter wheat field ( Figure 1b). S1 was used to collect LAI and DAM, and S2 was used to collect other crop parameters, including LAI, crop height, phenology, and final DAM. The size of each subfield is 100 m by 200 m. The blue points were the sampling locations in S1, and the red points were the sampling locations in S2.
Fieldwork was conducted multiple times from 11 May to 22 July in 2019. In S1, destructive aboveground biomass samples were collected at each of the 12 sampling locations on 8 May, 17 May, 21 May, 27 May, 3 June, 11 June, and 20 July. Winter wheat plant samples were collected from two 0.5 m by 0.5 sections within a 4 × 4 m area at each sampling location. The fresh plant samples were placed in plastic bags and transferred back to the lab directly. All samples were oven dried at 80 °C for at least 24 hours to obtain the DAM. In S2, 32 sampling points were used to collect other data, including LAI, soil moisture, crop height and phenology on 11 May, 21 May, 27 May, 3 June, and 11 June. At each sampling location, LAI was obtained using a Nikon D300s camera and a 10.5mm fisheye lens following the procedures described in Shang et al. (2014). Crop phenology was identified in the field using the Biologische Bundesanstalt, Bundessortenamt und CHemische Industrie (BBCH) scale. Details on data collection are listed in Table 1. (b) shows the location of S1 and S2 in the wheat field. The blue points are the sampling locations in S1, and the red points are the sampling locations in S2.

Field Sampling Design and Field Data Collection
Since the UAV imagery was to be collected at very high spatial resolution, the corresponding destructive biomass collection will affect the subsequent crop parameter measurements. To circumvent this situation, two sub-fields, S1 and S2, were selected within this winter wheat field (Figure 1b). S1 was used to collect LAI and DAM, and S2 was used to collect other crop parameters, including LAI, crop height, phenology, and final DAM. The size of each subfield is 100 m by 200 m. The blue points were the sampling locations in S1, and the red points were the sampling locations in S2.
Fieldwork was conducted multiple times from 11 May to 22 July in 2019. In S1, destructive aboveground biomass samples were collected at each of the 12 sampling locations on 8 May, 17 May, 21 May, 27 May, 3 June, 11 June, and 20 July. Winter wheat plant samples were collected from two 0.5 m by 0.5 sections within a 4 × 4 m area at each sampling location. The fresh plant samples were placed in plastic bags and transferred back to the lab directly. All samples were oven dried at 80 • C for at least 24 h to obtain the DAM. In S2, 32 sampling points were used to collect other data, including LAI, soil moisture, crop height and phenology on 11 May, 21 May, 27 May, 3 June, and 11 June. At each sampling location, LAI was obtained using a Nikon D300s camera and a 10.5mm fisheye lens following the procedures described in Shang et al. (2014). Crop phenology was identified in the field using the Biologische Bundesanstalt, Bundessortenamt und CHemische Industrie (BBCH) scale. Details on data collection are listed in Table 1.

Combine Harvester Yield Data Collection
The true spatially variable winter wheat yield data was collected by the producer using a 10-m swath John Deere Combine Harvester (John Deere, USA). This harvester equipped with a grain yield monitor and real-time kinematic (RTK)-global navigation satellite systems (GNSS) to record the dry grain weight every second and measured the harvested mass flow, moisture content, and geographic position, in addition to generating high spatial resolution yield map. The moving area of the combine harvester will be calculated by multiplying the harvester moving speed and time interval during crop weight recording. The yield map for S2 was composed of many points containing the yield values ( Figure 2). The yield data were then processed in ArcMap 10.7 (ESRI, USA) for duplicate-point removal, spatial resampling, and data extraction.

Combine Harvester Yield Data Collection
The true spatially variable winter wheat yield data was collected by the producer using a 10-m swath John Deere Combine Harvester (John Deere, USA). This harvester equipped with a grain yield monitor and real-time kinematic (RTK)-global navigation satellite systems (GNSS) to record the dry grain weight every second and measured the harvested mass flow, moisture content, and geographic position, in addition to generating high spatial resolution yield map. The moving area of the combine harvester will be calculated by multiplying the harvester moving speed and time interval during crop weight recording. The yield map for S2 was composed of many points containing the yield values ( Figure 2). The yield data were then processed in ArcMap 10.7 (ESRI, USA) for duplicate-point removal, spatial resampling, and data extraction.

UAV-Based Image Collection
Multi-temporal UAV-based RGB imagery was collected using a DJI phantom 4 RTK UAV system on 11 May, 21 May, and 27 May when crop phenology was at BBCH scale were 21, 31, and 39 of early leaf development stage to the end of the stem extension stage for S2. A 5K high-resolution digital camera was mounted on this system to collect information on red, green, and blue bands. An RTK base station was placed on the ground and combined with the RTK system on the UAV to achieve high precision location estimation for all imagery. The UAV flights were performed between 10 am and 2 pm and flown at the altitude of 30 m with the front and side overlapping of 90%. The total time of the operation is 55 min to acquire all images for S2. The Pix4Dmapper Pro v2.4 (Pix4D, Lausanne, Switzerland) software was used to process the UAV-based imagery to generate 3D point cloud data (PCD). The time of point cloud generation and LAIe calculation depends on the hardware of computer. It took about 30 h with a 12 core XEON processor and Quadro M4000 video card for this field in this study. The output of the 3D PCD has a similar format to LiDAR data, which contained the crop structure and optical RGB information. This type of PCD has a low cost on both sensor and UAV system.

Simulated Observation of Point Cloud Method
The SOPC method was designed to retrieve the spatial distribution of crop canopy and bare ground points in a simulated observation area from the UAV-based PCD and generate high spatial resolution crop LAIe map [47]. First, the SOPC method divided the study area into many observation areas based on the final resolution of LAI map. In each simulate area, the crop vegetation and bare ground points in the UAV-based PCD will be classified into two groups. The gap fraction will be calculated from the ratio of crop canopy and bare ground points at multi-view angles. Finally, the LAIe will be calculated from the gap fraction in the simulate area. The general principle of the SOPC method is shown in Figure 3. This method achieves the crop canopy LAIe estimation from the UAV-based PCD instead of the traditional optical information, which has limited effects by the shadow and view angles. In addition, the SOPC method retrieves the LAIe estimates without ground based LAI measurements and has a good agreement with downward-looking digital hemispherical photograph method derived LAIe. This method can successfully retrieve winter wheat LAIe at early growth stages from leaf development to the booting stage. The SOPC-LAIe maps are shown in Figure 4.
Lausanne, Switzerland) software was used to process the UAV-based imagery to generate 3D point cloud data (PCD). The time of point cloud generation and LAIe calculation depends on the hardware of computer. It took about 30 hours with a 12 core XEON processor and Quadro M4000 video card for this field in this study. The output of the 3D PCD has a similar format to LiDAR data, which contained the crop structure and optical RGB information. This type of PCD has a low cost on both sensor and UAV system.

Simulated Observation of Point Cloud Method
The SOPC method was designed to retrieve the spatial distribution of crop canopy and bare ground points in a simulated observation area from the UAV-based PCD and generate high spatial resolution crop LAIe map [47]. First, the SOPC method divided the study area into many observation areas based on the final resolution of LAI map. In each simulate area, the crop vegetation and bare ground points in the UAV-based PCD will be classified into two groups. The gap fraction will be calculated from the ratio of crop canopy and bare ground points at multi-view angles. Finally, the LAIe will be calculated from the gap fraction in the simulate area. The general principle of the SOPC method is shown in Figure 3. This method achieves the crop canopy LAIe estimation from the UAV-based PCD instead of the traditional optical information, which has limited effects by the shadow and view angles. In addition, the SOPC method retrieves the LAIe estimates without ground based LAI measurements and has a good agreement with downward-looking digital hemispherical photograph method derived LAIe. This method can successfully retrieve winter wheat LAIe at early growth stages from leaf development to the booting stage. The SOPC-LAIe maps are shown in Figure 4.

Weather Data
Weather data was retrieved from a nearby weather station located on the main campus of Western University, London, Ontario. This weather station has been in operation since 2016 and collects weather data every 30 min, including solar radiation (MJ/m 2 ), temperature (°C), rainfall (mm),

Weather Data
Weather data was retrieved from a nearby weather station located on the main campus of Western University, London, Ontario. This weather station has been in operation since 2016 and collects weather data every 30 min, including solar radiation (MJ/m 2 ), temperature ( • C), rainfall (mm), and wind speed (m/s). The weather data was used to represent the weather conditions in the region of the study area. The distance between the weather station and the study site is 35 km. The daily shortwave solar radiation from 1 October 2018 to 31 July 2019, was also extracted as the sum of daily solar radiation for the study site ( Figure 5a). The daily mean temperature was also calculated from the average of the daily maximum and minimum air temperature (Figure 5b).

Weather Data
Weather data was retrieved from a nearby weather station located on the main campus of Western University, London, Ontario. This weather station has been in operation since 2016 and collects weather data every 30 min, including solar radiation (MJ/m 2 ), temperature (°C), rainfall (mm), and wind speed (m/s). The weather data was used to represent the weather conditions in the region of the study area. The distance between the weather station and the study site is 35 km. The daily shortwave solar radiation from 1 October 2018 to 31 July 2019, was also extracted as the sum of daily solar radiation for the study site ( Figure 5a). The daily mean temperature was also calculated from the average of the daily maximum and minimum air temperature (Figure 5b).

SAFY Model Calibration
The SAFY model has been used to estimate winter wheat DAM [34]. This model determines the optimized biomass production in the crop growing season based on the LUE [41] and leaf partitioning function [42] theories. Firstly, the daily DAM (∆DAM) accumulation was calculated using the simple LUE theory with the equations shown below: where ELUE is the effective LUE, which is the LUE under environmental stress except temperature stress [36], Rg is the incoming shortwave solar radiation, and ε C is the climate coefficient, which is the ratio of photosynthetically active radiation (PAR) to the shortwave solar radiation. In this study, the fixed value was adopted, ε C . = 0.48 [33,49,50], f APAR is the fraction of absorbed photosynthetically active radiation; F T (T a ) is the temperature stress, and the light-interception coefficient k is 0.5 under the assumption of the leaf angle distribution is uniform and the leaf inclination is a spherical distribution [44,51].
Secondly, the daily increase of GLAI (∆GLAI) can be calculated from ∆DAM which is portioned to leaves (P L ) according to a given coefficient of specific leaf area (SLA). The equation given below: where the T a is the sum of optimal air temperature accumulated since plant emergence. P L is the fraction between leaf and dry aboveground biomass, which is determined by air temperature and another two parameters (partition to leaf function parameters P La and P Lb ) [42]. The equation can be written as follows: After the air temperature reached the threshold S TT , the daily GLAI can be calculated from the following equation, where Rs is the rate of senescence. After simulating the final DAM, the final crop yield can be calculated by multiplying the harvest index with the DAM. The harvest index (HI) was calculated from the ground biomass and final yield in S1, the average HI for all 12 sample points was 0.45. The yield can then be calculated using the equation below:

Winter Wheat Parameter Estimation from Ground-Based Biomass Measurement
The first step was attempted to determine the CSPs (P La , P Lb , S TT , Rs) of winter wheat and the range of ELUE based on the experimental DAM data collected in S1. The CSPs depend on the genetic characteristics of the type and cultivar of the winter wheat. Five parameters affect the biomass partitioning which include two parameters P La and P Lb in the partition to leaf function P L (Equation (4)), the sum of temperature to start senescence S TT ( • C), rate of senescence Rs ( • C/day), and ELUE is the ratio of photochemical energy produced as DAM from absorbed PAR (APAR).
Nine parameters identified in the literature, weather station measurements, and in-situ measurements ( Table 2) were used to calibrate the SAFY model and determine the crop CSPs [34]. The nine parameters include (1) climatic efficiency (ε C ), which is the ratio of PAR to the shortwave solar radiation; (2) minimum, optimal and maximum temperature (T min , T opt , and T max ) required for winter wheat growth; (3) specific leaf area (SLA), which is the unit weight of crop leave; (4) initial value of DAM at the day of plant emergence; (5) the light-interception coefficient (k ext ) in Bear's law which is related to the plant leaf area index and the fraction of APAR; (6) the day of plant emergence; (7) the day of senescence; (8) the daily shortwave solar radiation (Rg); and (9) daily mean air temperature (T a ). Detailed values used in the SAFY model are given in Table 2. The five CSPs are calibrated in the first SAFY calibration.
The winter wheat CSPs and the range of ELUE were calibrated against the DAM observation collected in S1 using the global optimization method, Shuffled Complex Evolution-University of Arizona (SCE-UA) algorithm [52], to determine the optimal value for P La , P Lb , S TT , Rs, and ELUE [5,34,36].
After running the optimization procedure using the SCE-UA algorithm (Duan, Sorooshian, & Gupta, 1994), the winter wheat CSPs and ELUE can be determined. The RMSE between the simulated DAM (DAM sim ) and the In-situ true DAM (DAM true ) was used as the cost function of the calibration: Remote Sens. 2020, 12, 2378 9 of 21 Model calibration was performed by minimizing the RMSE between DAM sim and DAM true . In the SAFY model, the optimization procedure was run 10,000 times for each sampling location to achieve the optimal parameters with the lowest RMSE. In addition, RRMSE was also used to indicate model accuracy.

Fisheye-Derived GLAI and Model-Simulated GLAI
As the first calibration generated the daily GLAI values for 12 sampling locations for the entire growing season using the SAFY model, these GLAI values can be used to calibrate the fisheye derived LAIe. The relationship between the fisheye LAIe and the GLAI was established using the data from S1. This relationship was then used to convert the fisheye derived LAIe in S2 and used as input for the SAFY model to simulate the DAM and SAFY-GLAY in S2. A second calibration of the SAFY model was performed for S2 using the winter wheat CSPs (P La , P Lb , S TT , Rs), and the GLAI which is converted from fisheye-derived LAIe measurements in S2 to simulate the SAFY-GLAI and final DAM. The median values of these CSPs were then adopted and the ELUE was kept as a variable in the second calibration. The RMSE between the simulated GLAI (GLAI sim ) and the converted fisheye LAI measurements (GLAI true ) were used as the cost function in the SCE-UA algorithm during the calibration.
The calibration procedure was also performed 10,000 times to ensure the optimal GLAI simulation with the lowest RMSE.

Final DAM and Yield Estimation Using UAV-Based LAIe in S2
The relationship between the UAV-based LAIe and simulated SAFY-GLAI derived from the second SAFY calibration was then established from the second calibration. After converting the UAV-based LAIe to GLAI, the third SAFY calibration was performed for S2 using the winter wheat CSPs derived from the first SAFY calibration and the UAV-based GLAI to simulate the DAM in S2. The total of UAV-based LAIe has 5977 measurements for S2 on each monitoring date. The final DAM was estimated by optimizing the lowest RMSE of GLAI value using Equation (9). The final yield was calculated using Equation (6). After resampling the UAV-based final yield to the same resolution of the harvester yield data, the mean, standard deviation, and coefficient of variation (CV) were used to evaluate the performance of yield estimation. The flow chart below illustrates the steps of winter wheat yield estimation using the UAV-based LAI ( Figure 6).
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 22 Figure 6. Flowchart shows the steps to perform UAV-based winter wheat yield estimation.
Step 1 determines the cultivar-specific parameters of winter wheat.
Step 1 and 2 established the relationship between simulated green leaf area index (GLAI) and UAV-based point cloud LAIe.
Step 3 applies the UAV-based GLAI maps and cultivar-specific parameters in simple algorithm for yield (SAFY) model to estimate final biomass and yield map.

Determination of Cultivar-Specific Parameters
After the initial SAFY model calibration using the In-situ DAM, the ranges of , , , and ELUE were determined. Table 3 shows the results of the five parameters, including maximum, minimum, median, mean, and standard deviation (STD). The median value of , , , and were adopted from the second and third SAFY model calibration [34]. The RMSE between the In-situ DAM measurements and simulated DAM in S1 was 81 g/m 2 , and the RRMSE is 13.89%. Step 1 determines the cultivar-specific parameters of winter wheat.
Step 1 and 2 established the relationship between simulated green leaf area index (GLAI) and UAV-based point cloud LAIe.
Step 3 applies the UAV-based GLAI maps and cultivar-specific parameters in simple algorithm for yield (SAFY) model to estimate final biomass and yield map.

Determination of Cultivar-Specific Parameters
After the initial SAFY model calibration using the In-situ DAM, the ranges of P La , P Lb , S TT , Rs and ELUE were determined. Table 3 shows the results of the five parameters, including maximum, minimum, median, mean, and standard deviation (STD). The median value of P La , P Lb , S TT , and Rs were adopted from the second and third SAFY model calibration [34]. The RMSE between the In-situ DAM measurements and simulated DAM in S1 was 81 g/m 2 , and the RRMSE is 13.89%. Table 3. The cultivar-specific parameters and effective light use efficiency (ELUE) derived from the initial SAFY calibration. P La is the parameter a of P L function; P Lb is the parameter b of P L function; S TT ( • C) is the sum of temperature for senescence; Rs ( • C day) is the rate of senescence; ELUE (g/MJ) is the effective light-use efficiency.

Relationship between Simulated SAFY-GLAI and Fisheye-Derived LAIe in S1 and S2
After the first SAFY calibration ( Figure 6, step 1), the simulated daily SAFY-GLAI were generated from DAM calibration for 12 sampling locations in S1, and the relationship between the simulated GLAI and fisheye-derived LAIe were established (Figure 7). The coefficient of determination (R 2 ) was 0.75 for all 60 measurements.

Relationship between Simulated SAFY-GLAI and Fisheye-Derived LAIe in S1 and S2
After the first SAFY calibration ( Figure 6, step 1), the simulated daily SAFY-GLAI were generated from DAM calibration for 12 sampling locations in S1, and the relationship between the simulated GLAI and fisheye-derived LAIe were established (Figure 7). The coefficient of determination (R 2 ) was 0.75 for all 60 measurements. Using this relationship between fisheye-derived LAIe and simulated GLAI in S1, the fisheyederived LAIe measurements in S2 were converted to simulated GLAI. The second SAFY model calibration ( Figure 6, step 2) was conducted using the simulated GLAI values of 32 sampling locations in S2, the final DAM and simulated daily SAFY-GLAI of these locations was determined in the second SAFY calibration. After estimating the final DAM for the 32 sampling locations in S2 from the second calibration, the In-situ measured final DAM and simulate DAM derived from SAFY model were compared (Figure 8). The R 2 between the estimation and measured grain yield was 0.48, the RMSE was 54 g/m 2 , and the RRMSE was 9.37%.  Using this relationship between fisheye-derived LAIe and simulated GLAI in S1, the fisheye-derived LAIe measurements in S2 were converted to simulated GLAI. The second SAFY model calibration ( Figure 6, step 2) was conducted using the simulated GLAI values of 32 sampling locations in S2, the final DAM and simulated daily SAFY-GLAI of these locations was determined in the second SAFY calibration. After estimating the final DAM for the 32 sampling locations in S2 from the second calibration, the In-situ measured final DAM and simulate DAM derived from SAFY model were compared (Figure 8). The R 2 between the estimation and measured grain yield was 0.48, the RMSE was 54 g/m 2 , and the RRMSE was 9.37%. Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 22 Figure 8. Relationship between the measured and estimated dry above-ground biomass (DAM) using SAFY model for S2.

DAM Estimation Using UAV-based LAIe Measurements
After the second calibration of SAFY model in S2, the simulated daily SAFY-GLAI were generated from the GLAI (converted from fisheye-LAIe) for 32 sampling locations. By comparing the results of simulated SAFY-GLAI and the SOPC derived UAV-based LAIe of 32 sampling location, the relationship between the simulated SAFY-GLAI and UAV-based LAIe was established for three monitoring days and shown in Figure 9. The R 2 was 0.82 for all 96 measurements. Based on this relationship, the multiple UAV-based LAIe maps will be converted to GLAI maps. These GLAI maps will be used in the third SAFY calibration ( Figure 6, step 3) to estimate the final DAM in S2.

DAM Estimation Using UAV-Based LAIe Measurements
After the second calibration of SAFY model in S2, the simulated daily SAFY-GLAI were generated from the GLAI (converted from fisheye-LAIe) for 32 sampling locations. By comparing the results of simulated SAFY-GLAI and the SOPC derived UAV-based LAIe of 32 sampling location, the relationship between the simulated SAFY-GLAI and UAV-based LAIe was established for three monitoring days and shown in Figure 9. The R 2 was 0.82 for all 96 measurements. Based on this relationship, the multiple UAV-based LAIe maps will be converted to GLAI maps. These GLAI maps will be used in the third SAFY calibration (Figure 6

DAM Estimation Using UAV-based LAIe Measurements
After the second calibration of SAFY model in S2, the simulated daily SAFY-GLAI were generated from the GLAI (converted from fisheye-LAIe) for 32 sampling locations. By comparing the results of simulated SAFY-GLAI and the SOPC derived UAV-based LAIe of 32 sampling location, the relationship between the simulated SAFY-GLAI and UAV-based LAIe was established for three monitoring days and shown in Figure 9. The R 2 was 0.82 for all 96 measurements. Based on this relationship, the multiple UAV-based LAIe maps will be converted to GLAI maps. These GLAI maps will be used in the third SAFY calibration (Figure 6, step 3) to estimate the final DAM in S2.     Figure 11 shows the final DAM maps at the resolution at 2 × 2 m.
Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 22 same spatial resolution as the UAV-based LAIe map. Figure 11 shows the final DAM maps at the resolution at 2 × 2 m.
(a) (b) (c)   Table 4 shows statistical information of the harvester and estimated yield maps. Figure 12 shows the maps of harvester measured grain yield (Figure 12a) and estimated yield (Figure 12b) for S2. The accuracy of the estimated yield was evaluated by comparing their RMSE, mean, standard deviation (STD), and CV. Figure 13 shows the absolute difference map between true grain yield and estimated yield for S2. same spatial resolution as the UAV-based LAIe map. Figure 11 shows the final DAM maps at the resolution at 2 × 2 m.

Comparison of True Grain Yield and Estimated Yield
(a) (b) ©   Table 4 shows statistical information of the harvester and estimated yield maps. Figure 12 shows the maps of harvester measured grain yield (Figure 12a) and estimated yield (Figure 12b) for S2. The accuracy of the estimated yield was evaluated by comparing their RMSE, mean, standard deviation (STD), and CV. Figure 13 shows the absolute difference map between true grain yield and estimated yield for S2.  Table 4 shows statistical information of the harvester and estimated yield maps. Figure 12 shows the maps of harvester measured grain yield (Figure 12a) and estimated yield (Figure 12b) for S2. The accuracy of the estimated yield was evaluated by comparing their RMSE, mean, standard deviation (STD), and CV. Figure 13 shows the absolute difference map between true grain yield and estimated yield for S2. Table 4. The mean grain yield, coefficient of variation (CV), and standard deviation (STD) of grain yield measured by harvester and estimated by SAFY model. The root mean square error and relative root mean square error between the harvester and estimated yield.

Cultivar-Specific Parameters Derived from the First SAFY Model Calibration
The determination of the CSPs (P La , P Lb , S TT , Rs and ELUE) can help the SAFY model in crop final biomass estimation. Many studies have provided the ranges of these parameters. For instance, Duchenmin et al. (2008) determined the range of P La and P Lb for winter wheat as 0.05-0.5 and 10 −5 -10 −2 respectively. The median values of P La = 0.1573 and P Lb = 0.00196 were adopted during the final DAM estimation. In this study, the ranges of P La and P Lb used were 0.2038-0.2686 and 0.0015-0.0021, respectively. The ranges of the two parameters used in our study are much smaller in comparison with that of in the literature. This is because the In-situ DAM measurements were used to calibrate the SAFY model and the day of plant emergence was observed through In-situ observation. The P La and P Lb combined with the accumulated temperature will affect the value of P L (Equation (4)). Figure 14 shows the relationship between P L and the accumulated temperature in this study. The P L decreases exponentially from the value of 1 to the value of 0 with the accumulated temperature, which is the plant emergence to the end of the leave production phase [42]. The P L is the ratio of the daily increase of GLAI and the daily increase of DAM, the value of 0 marks the stopping point of leave developments. As illustrated in Figure 14, P L is zero when the accumulated temperature is at 881 • C, and the date was 14 June. According to in-situ observation on 16 June, winter wheat was at the end of flowering stage at that time.
The determination of the CSPs ( , , , and ELUE) can help the SAFY model in crop final biomass estimation. Many studies have provided the ranges of these parameters. For instance, Duchenmin et al. (2008) determined the range of and for winter wheat as 0.05-0.5 and 10 −5 -10 −2 respectively. The median values of = 0.1573 and = 0.00196 were adopted during the final DAM estimation. In this study, the ranges of and used were 0.2038-0.2686 and 0.0015-0.0021, respectively. The ranges of the two parameters used in our study are much smaller in comparison with that of in the literature. This is because the In-situ DAM measurements were used to calibrate the SAFY model and the day of plant emergence was observed through In-situ observation. The and combined with the accumulated temperature will affect the value of (Equation (4)). Figure 14 shows the relationship between and the accumulated temperature in this study. The decreases exponentially from the value of 1 to the value of 0 with the accumulated temperature, which is the plant emergence to the end of the leave production phase [42]. The is the ratio of the daily increase of GLAI and the daily increase of DAM, the value of 0 marks the stopping point of leave developments. As illustrated in Figure 14, is zero when the accumulated temperature is at 881 °C, and the date was 14 June. According to in-situ observation on 16 June, winter wheat was at the end of flowering stage at that time. The simulated total accumulated temperature ( ) determined by SAFY was 954 °C when senescence started. The date for achieving this temperature was 17 June, which was one day later than winter wheat flowering stage. According to the relationship between phenology and total accumulated temperature, the winter wheat starts seed fill at this temperature [54]. Beyond this threshold temperature, the GLAI calculated using Equation (5) started to decrease. In Equation (5), the is related to the threshold temperature and total accumulated temperature for maturity. The total accumulated temperature can be used to determine the date of maturity with the aid of the weather data. Based on the value of and , the total accumulated temperature for maturity of all 12 samples were evaluated using the SAFY model for S1. The average total mature day is 282, which is five days earlier than the day of harvest (287). The values of the sum of temperature for senescence ( ) and the rate of senescence ( ) are consistent with the actual in-situ observations. The simulated total accumulated temperature (S TT ) determined by SAFY was 954 • C when senescence started. The date for achieving this temperature was 17 June, which was one day later than winter wheat flowering stage. According to the relationship between phenology and total accumulated temperature, the winter wheat starts seed fill at this temperature [54]. Beyond this threshold temperature, the GLAI calculated using Equation (5) started to decrease. In Equation (5), the Rs is related to the threshold temperature and total accumulated temperature for maturity. The total accumulated temperature can be used to determine the date of maturity with the aid of the weather data. Based on the value of S TT and Rs, the total accumulated temperature for maturity of all 12 samples were evaluated using the SAFY model for S1. The average total mature day is 282, which is five days earlier than the day of harvest (287). The values of the sum of temperature for senescence (S TT ) and the rate of senescence (Rs) are consistent with the actual in-situ observations.

ELUE
ELUE plays an important role in the SAFY model. ELUE is largely influenced by crop species, physiological factors, soil conditions, and weather conditions. Soil conditions such as soil nutrient, soil moisture, texture, organic matter, and pH vary across the field and can lead to variable ELUE values for the same crop. Conventionally, a constant ELUE value is used for the entire field in yield estimation. Dong et al. (2017) introduced the spatially variable crop maximum LUE for the first time and achieved significant improvements in biomass estimation accuracy for winter wheat and corn. Liao et al. (2019) also used variable ELUE values to estimate the yield of corn and soybean. Therefore, in this study the ELUE was adopted as a variable parameter in the SAFY model rather than a fixed value as the P La , P Lb , S TT , and Rs in the second and third calibrations. The range of ELUE in the first calibration was 2.93-3.18, which fits in the range of ELUE in the literature [5,34].

Uncertainties of the Estimated Crop Biomass and Yield
By comparing the yield maps in Figure 12, both the estimated and measured maps have similar patterns at the bottom corner, but the upper left area in the field has different yield estimations. Figure 13 shows the absolute difference between the true grain yield and estimated yield for S2. Since the disease and pests stress were not observed during the fieldwork in the winter wheat field. The absolute difference in the upper left area in the field may be affected by some external factors such as soil water content and nutrient. Because the SAFY model simplifies the crop growth model which considers the temperature stress only, the soil water and nutrient stress after the last LAI monitoring may lead to an inaccurate final DAM estimation. Except the upper left area, most of the test points have small differences which are less than 100 g/m 2 . Some test points with large yield differences are located at the end of rows which may cause by the basis of the harvester. Because the harvester needs to lift the head during turning. The uncertainties of the yield estimation on the upper left area might have been due to the limited number of UAV-based LAIe maps. The early LAIe measurements of the winter wheat were used in the SAFY model to estimate the final DAM and yield. The last UAV-based LAIe measurement was 50 days earlier than the actual harvest date. The LAIe information used only captured the crop growth condition up to the date of the last UAV flight which could not have taken into condition of the further biomass growth during the later season. The multiple data sources can be considered to increase the temporal resolution of the GLAI for the SAFY model.
The histograms of the UAV-based yield and the harvester measured yield map are shown in Figure 15. Overall, these two yield histograms exhibit a similar distribution, in which the range of harvester measured yield data was between 340 to 840 g/m 2 , and the range of the estimate yield data was between 420 to 780 g/m 2 . The UAV estimated yield has a lower standard deviation and CV than that of the harvester measured yield data in S2. This is likely due to the fixed P La and P Lb used in the SAFY model. In addition, the ground-based yield map, generated by the combine harvester, might have suffered yield loss in comparison with manual harvest. The destructive biomass is more accurate in estimating GLAI and yield during the SAFY model calibrations. This may be one of the reasons for the difference between the true and estimate yield maps. Due to the use of destructive biomass, the estimated yield has a slightly narrower histogram distribution in the yield map. In addition, the final yield might have been influenced by the fixed harvest index. The UAV-based yield was calculated from the average harvest index derived from 12 sampling points. However, the crop growth condition can influence the harvest index and lead to a different yield estimation [55].

Application and Contribution
This study, for the first time, applied the SOPC derived UAV-based point cloud LAIe maps to the SAFY model to generate the sub-field biomass and yield maps. It is also the first time to use the UAV-based data to the SAFY model. One of the potential applications of this framework is for yield

Application and Contribution
This study, for the first time, applied the SOPC derived UAV-based point cloud LAIe maps to the SAFY model to generate the sub-field biomass and yield maps. It is also the first time to use the UAV-based data to the SAFY model. One of the potential applications of this framework is for yield estimation. This can help reveal the spatial variability of yield potential for winter wheat at the sub-field scale. The LAI maps generated at various growth stages can also provide useful information to assist in making crop management decisions. The soft red winter wheat CSPs determined in this study can be used in future studies for the same region. Since this framework is designed for UAV application at farm level, the day of planting and emergence can be obtained from in-situ observation. With specific phenological date and the CSPs, yield prediction using the SAFY model and the UAV derived simulated GLAI only requires the solar radiation and temperature data. Therefore, it is possible to predict the winter wheat final yield using weather data collected at real-time and simulated for the rest of the growing season using the proposed framework.
The normalization of LAIe derived from different platforms can help the SAFY model application to final DAM estimation using different data sources. The normalization approach can be applied to other crops to determine the CSPs. In this study, the normalized UAV-based point cloud LAIe were used in the SAFY model calibration for final DAM and yield estimation. The UAV-based point cloud LAIe can provide a variable range of spatial resolution map from submeter to meter scale (50 cm to 5 m) which can clearly display the within-field final DAM and yield variability. Fu et al. (2020) [45] estimated the winter wheat biomass using the UAV-based multispectral imagery. The RRMSE of the DAM were 23.37%. In this study, the DAM values were estimated from only 30 individual sampling locations. In contrast, this study achieved a RRMSE of 15% for more than 1800 points. In addition, this study uses the UAV-based PCD derived from RGB imagery instead of multispectral vegetation indices. The SOPC method does not require any ground spectral calibration, which is more accessible and low cost. The final DAM and yield estimation could be conducted without ground measurements in this study. Furthermore, the method of this study could be applied on different geographical locations and other crops. The solar radiation and air temperature should be collected for the different locations. The SOPC-LAIe calculates the crop leaf area index based on the gap fraction theory under the assumption of the leaf angle distribution is uniform and the leaf inclination follows spherical distribution [44,51]. If the crop leaf area index can be calculated using the gap fraction theory with the same assumption. This method is practical for different locations and crops. For example, the gap fraction theory on leaf area index estimation and the SAFY model on yield estimation have been proved on corn and soybean [36]. This method may possible be applied on corn and soybean and to estimate the final yield.

Conclusions
In this study, destructive biomass was used to calibrate the SAFY model to derive CSPs for winter wheat instead of using the LAI measurements to calibrate the SAFY model for determining of CPSs. Through the normalization of the SOPC-LAIe maps to the simulated SAFY-GLAI maps, the UAV-based point cloud derived LAIe can be used as input to the SAFY model for winter wheat final DAM and yield estimation. This is the first time the SOPC derived LAIe was used for final DAM and yield estimation of winter wheat. It is also the first time to apply the UAV-based LAI to the SAFY model to estimate the final DAM. The results showed that the SOPC method derived UAV-based point-cloud LAIe and the SAFY model have a great potential in generating high-spatial resolution (2 × 2 m in this study) DAM and yield map. The final yield estimation achieved the RMSE of 88 g/m 2 and RRMSE of 15.22%. The UAV-based point cloud derived LAIe before the booting stages can be used to estimate the final DAM and yield of winter wheat. The CSPs determination approach developed in this study can be adopted in deriving crop CSPs for other crop types when using SAFY model to estimate and forecast final DAM and grain yield. After determining the CSPs of crops, the approach can be achieved without having to rely on ground measurements, which is a great advantage for operational near-real-time operations.
Given the success of the reported results, this stud still has the limitations. For future work, UAV imagery collected at the late growth stages can be incorporated into the analysis to achieve improved estimation accuracy. In addition, more detailed soil information, such as moisture content, nutrient level, and soil organic matter, should also be used as model input to improved model performance. Funding: This work is partially funded by NSERC discovering Grant awarded to Wang; a MITACS accelerate internship sponsored by MITACS and A&L Canada Laboratories Inc., and an NSERC CRD grant awarded to Leblon, Wang and Sabarinathan.