Estimation of Winter Wheat Yield from UAV-Based Multi-Temporal Imagery Using Crop Allometric Relationship and SAFY Model

: Crop yield prediction and estimation play essential roles in the precision crop management system. The Simple Algorithm for Yield Estimation (SAFY) has been applied to Unmanned Aerial Vehicle (UAV)-based data to provide high spatial yield prediction and estimation for winter wheat. However, this crop model relies on the relationship between crop leaf weight and biomass, which only considers the contribution of leaves on the ﬁnal biomass and yield calculation. This study developed the modiﬁed SAFY-height model by incorporating an allometric relationship between ground-based measured crop height and biomass. A piecewise linear regression model is used to establish the relationship between crop height and biomass. The parameters of the modiﬁed SAFY-height model are calibrated using ground measurements. Then, the calibrated modiﬁed SAFY-height model is applied on the UAV-based photogrammetric point cloud derived crop height and effective leaf area index (LAIe) maps to predict winter wheat yield. The growing accumulated temperature turning points of an allometric relationship between crop height and biomass is 712 ◦ C. The modiﬁed SAFY-height model, relative to traditional SAFY, provided more accurate yield estimation for areas with LAI higher than 1.01 m 2 /m 2 . The RMSE and RRMSE are improved by 3.3% and 0.5%, respectively.


Introduction
Crop biomass is an important parameter for crop yield potential prediction. Environmental conditions such as solar energy, temperature, soil nutrients, water, pests, disease, weeds, and other stresses can affect dry aboveground biomass (DAM) and potential yield of a crop [1]. As a crop management system, precision farming emphasizes monitoring and analysis to help farmers achieve optimal profit [2]. Hence, precision farming focused on targeted and cost-effective inputs to enhance efficiency and sustainability outputs. Therefore, mapping field-scale crop DAM play an essential role in precision farming. Conventional crop yield estimation methods are labor intensive and have limited sampling numbers due to the accessibility of fieldwork. These methods provide rough estimations and advising in farming activities. Remote sensing technology allows monitoring of variability for large-scale fields without a large amount of ground data. For example, satellite remote sensing has been used for regional crop monitoring. However, satellite imagery is hampered by long revisit times and cloud presence. Additionally, the cost of high temporal and spatial resolution satellite imagery limits its use. Unmanned Aerial Vehicle (UAV) collected high spatial and temporal imagery can monitor within-field DAM and yield variability in real-time. However, a practical and accurate approach is required to achieve within-field crop biomass monitoring and estimation using UAV-based data.
Two approaches have been adopted using remote sensing technology to estimate crop biomass. The first category is the statistical model, which uses the relationship between remote sensing vegetation indices and in-situ crop biomass measurements to estimate crop biomass [3][4][5]. The UAV-based high spatial resolution remote sensing data can indicate the crop biomass variation within the field using the vegetation indices. However, this method is applicable only for specific crop growth stages and geographic location, and variations in environmental conditions and location will affect the accuracy [6]. The second category is the crop growth model, which uses crop and environment parameters to simulate the development of crop growth and estimate crop biomass. Currently, crop growth models such as AquaCrop [7], CERES [6,8], STICS [9], CropSyst [10], and WOFOST [11] have been well developed in crop production estimation. However, these models require a comprehensive set of parameters to simulate crop growth status, which is difficult to collect and use for within-field DAM and yield estimation. WOFOST is a labor and time-intensive model which requires about 40 parameters, and may not be an affordable approach for UAV-based imagery in estimating within-field crop yield.
An allometric method is an alternative approach for estimating plant biomass from other crop parameters. Plants absorb solar energy and transfer this energy to biomass during the growing season. The allocation pattern of the biomass distribution between stem, leaves, fruit, and roots follows a special allocation principle for specific species. This pattern of energy allocation plays an essential role in the plant's reproduction, which affects the final yield in cereal crops [12]. The plant's allocation principle was defined as the allometric relationship in the plant, which is a part of the plant that has a similar variation rate compared with other parts of the plant or whole plant [13,14]. The information derived from the crop allometric relationships can also help farmers to monitor, analyze, and determine crop growth situations. Many remote sensing studies have proved that the allometric method can be applied to different crop species to determine the crop parameters [15][16][17][18]. The allometric method initially requires extensive destructive samplings to build an allometric model. The model could then be used as a non-destructive method to estimate crop biomass. Bakhshandeh et al., (2012) [14] addressed the allometric relationship between winter wheat plant height and biomass based on a ground-based study. However, limited ground sampling points presents a challenge for applying the allometric relationship to the entire field and displaying the within-field biomass variability.
Since intensive crop height data is either difficult or expensive to obtain from the conventional remote sensing methods such as satellite or airborne optical data directly, not many studies estimate crop biomass from plant height using remote sensing technology at a subfield scale. The recent development of the UAV system and computer vision have made possible dense 3D reconstructions to produce orthomosaic aerial image, DEM, and 3D photogrammetric point clouds using the Structure from Motion (SfM) method from UAV data. SfM is a computer vision technique that incorporates multi-view stereo images to match features, derive 3D structure, and estimate camera position and orientation. The SfM approach has been proven to be able to estimate of crop structural information such as, canopy size, height, and leave shape [19][20][21][22]. For example, Bendig (2014) [21] used a crop surface model to estimate the height of barley, highly correlated with the dry aboveground biomass. Thus, the allometric relationship between crop height and biomass could be meaningful information in estimating crop biomass.
A Simple Algorithm for Yield Estimation (SAFY) model uses the allometric relationship between daily crop leaf area index (LAI) and biomass accumulation to simulate the crop growth stages. In addition, it uses light use efficiency theory to determine the optimal parameters and estimate the final dry aboveground biomass of crops [18,[23][24][25]. Song et al., (2020) [26] used the UAV-based Simulated Observation of Point Cloud (SOPC) derived effective leaf area index (LAIe) [27] and the SAFY model to generate the fine-resolution final DAM and yield maps of winter wheat. The SAFY model may be used to derive a more accurate final DAM and yield estimation by combining the allometric relationship between winter wheat canopy height and biomass. The objective of this study is to incorporate the allometric relationship between crop canopy height and biomass into the SAFY model and assess the accuracy of crop biomass estimation. A modified SAFY-height model will be developed and calibrated using ground-based measurements. The final DAM and yield will be estimated from the calibrated modified SAFY-height model using UAV-based canopy height and LAIe maps. The canopy height and LAIe maps of winter wheat will be generated using the moving cuboid filter [28] and SOPC method [27].

Study Area
The study site locates in the west of Melbourne, Southwest Ontario, Canada (42.787707 • , -81.594801 • ). Since the study region has a shorter growing season, which can support one harvest per year, the soft red winter wheat (Brevatn Branson, Corteva Agriscience, USA) was seeded in October 2018 and harvested in July 2019. Two sub-fields were selected to collect different parameters. The size of the sub-field is 100 m by 200 m. The sub-fields are shown in Figure 1. The blue points are the sampling points in sub-field 1 (S1), and the red points are the sampling points in sub-field 2 (S2). The study site in Canada; (b) the sub-field 1 (S1) and 2 (S2) within the winter wheat field. The blue points are the sampling points in S1; and the red points are the sampling points in S2.

Ground Data Collection
The ground crop parameters collections were divided into two sub-fields since the destructive DAM will affect the accuracy of the UAV-based LAIe and height data. These two sub-fields were selected from one winter wheat field with the same cultivar, phenology, and nutrition supplies. S1 was used to collect winter wheat DAM, canopy height, LAIe, and final yield in 12 sampling locations. The Normalized Difference Water Index (NDWI) map derived from the Sentinel 2 imagery on 22 April 2019 was used for the selection of sampling locations. The NDWI map represents the water distribution of our study site, which will lead to a significant variation of crop final yield. S2 was used to collect winter wheat canopy height, LAIe, and final yield in 32 sampling locations. S2 used systematic sampling method to evaluate model performance. In addition, the UAV-based imagery was collected in S2. Multiple days of field work were conducted from 8 May to 20 July in 2019. The detailed parameters and dates are shown in Table 1. The DAM was collected from two blocks of a 0.5 m by 0.5 m section within a 2 m by 2 m area of each sampling location. After collecting the fresh winter wheat sample, the sample was transferred to the lab and oven dried at 80 • C for 24 h. The final dry weight of crop biomass was measured in the lab. The LAIe was calculated from the digital hemispherical photograph using the CAN-EYE v6.1. The digital hemispherical photograph was collected by a Nikon D300s camera and a 10.5 mm fisheye lens following the procedures described in Shang et al. (2014) [29]. The average crop canopy height was averaged from 10 measurements around the sampling location within a 2-m radius.

Combine Harvester Yield Data Collection
The crop final yield data comes from a John Deere combine harvester at the end of the season. The combine harvester had a 10 m wide head to collect the plant. The built-in real-time kinematic (RTK)-global navigation satellite system (GNSS) recorded the distance of the harvester moved every second. The area of the harvester moved will be calculated using the width of the head and driven distance. The grain yield monitor on the harvester recorded the grain weight, mass flow, and moisture content. The spatial yield map was composed of many points, which were shown in Figure 2. ArcMap 10.7 (ESRI, Redlands, CA, USA) was used to process the final yield data, including outlier removal, resampling, and data extraction.

UAV-Based Image Collection
A DJI Phantom 4 RTK UAV system was used to collect RGB imagery for S2 on 11 May, 21 May, and 27 May. The crop phenology for these dates at the BBCH scale was 21, 31, and 39, which covered the leaf development stage to the booting stage. The Phantom 4 RTK system can provide a 20MPix high-resolution imagery, and a ground RTK base station can achieve very high accuracy positional information. All UAV flights were performed at the altitude of 30 m with front and side overlapping of 90%. The UAV images have been processed to generate 3D point cloud data using Pix4Dmapper Pro v2.4 (Pix4D, Lausanne, Switzerland). The total processing time of photogrammetric point cloud generation is 30 h using a Windows computer with a 12 core Xeon cup and Quadro m4000 graphic card. The final resolution of the aerial image is 0.8 cm, and the point density of point cloud data is more than 4000 points/m 2 .

UAV-Based Plant Height and LAIe Maps
The UAV-based 3D point cloud data were used to generate the plant height and LAIe maps for three UAV operation days. The recently developed moving cuboid filter was used to retrieve the winter wheat plant height using the UAV-based 3D photogrammetric point cloud data without in-situ measurements [28]. The UAV-based photogrammetric point cloud data has very high point density with more than 4000 points/m 2 . The final resolution of the products derived from photogrammetric point cloud data can be adjusted based on the research requirement [27,28]. The spatial resolution of canopy height and LAIe maps in this study was 2 m. In addition, the photogrammetric point cloud data can generate canopy height and LAIe maps of winter wheat with high accuracy in the early growth stages. The canopy height maps of winter wheat in S2 were shown in Figure 3. The root mean square error (RMSE) for each canopy height map was 2.61 cm, 3.65cm, and 4.17 cm for the three dates respectively. The overall RMSE for all ground measurements on three dates was 3.48 cm. The LAIe maps of winter wheat in S2 were shown in Figure 4. The overall RMSE for all ground measurements on three dates was 0.19. The recent developed SOPC method was designed to estimate winter wheat LAIe from UAV-based 3D point cloud data using the gap fraction theory without in-situ measurements [26]. This method can successfully generate winter wheat LAIe maps at early growth stages from the leaf development to the booting stage. The UAV-based SOPC derived LAIe (SOPC-LAIe) maps are shown in Figure 4.

Weather Data
Weather data comes from a weather station located on the main campus of Western University, London, Ontario, which is 35 km away from the study site. The weather station collected the shortwave solar radiation and air temperature every 30 min. The daily shortwave solar radiation and the average air temperature were calculated; the data chart is shown in Figure 5.

Allometric Relationship Establishment
A piecewise linear regression model has been used to determine the relationship between plant height and DAM on wheat [14]. The wheat growth can be separated into two phases. The first phase covers the tillering and stems elongation stages; the second phase covers the rest of the stages beyond the stem elongation stages. The equations are shown below: where DAM is the final dry aboveground biomass, H is the plant height, α 1 is the increasing rate of DAM to H in the first phase, b 1 is the intercept in phase 1, a 2 is the increasing rate of DAM to H in the second phase, b 2 is the intercept in phase 2, the ∑ T a is the sum of optimal air temperature accumulated since plant emergence, T turn is the turning point between the two phases which is the total accumulated temperature.

SAFY Model
The SAFY model has been used to simulate the winter wheat daily DAM and green leave area index (GLAI) during the growing season based on the crop LUE and leaf partitioning function theories [18,30,31]. The LUE theory is modelled below: where ELUE is the effective LUE; Rg is the incoming shortwave solar radiation; ε C is the climate coefficient, which ε C = 0.48 [32][33][34]; 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 [35,36]; and F T (T a ) is the temperature stress. The daily increase of GLAI (∆GLAI) can be calculated from the fraction of the daily DAM (∆DAM) which is portioned to leaves (P L ) according to a given coefficient of specific leaf area (SLA). The leaf partitioning function shows below: where P L is the fraction between leaf and dry aboveground biomass, which is determined by air temperature and another two parameters (P La and P Lb ) [31]. The equation can be written as follows: After the air temperature reached the threshold S TT , the 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 final estimated yield was compared with the harvester measured yield to evaluate the performance of this study. The harvest index (H I) was calculated from the ground biomass and final yield in field 1, the average for all 12 sample points was 0.45. The equation shows below:

Modified SAFY-Height Model
The allometric relationship between crop canopy height and DAM is added to SAFY model to develop the modified SAFY-height model. This model is able to simulate a relationship between daily GLAI, canopy height, and DAM. First, the modified SAFY-height model will be calibrated using the ground-based DAM and canopy height measurements to determine the optimal parameters of the canopy height and DAM allometric relationship in S1 using the global optimization method, Shuffled Complex Evolution-University of Arizona (SCE-UA) algorithm [37]. In step one, Song et al. used the LAIe map derived from the UAV-based photogrammetric point cloud data to replace the GLAI and determined the cultivar-specific parameters (P La , P Lb , S TT , Rs) [26]. In step two, the modified SAFYheight model will combine the UAV-based canopy height map, UAV-based LAIe map, the cultivar-specific parameters, and allometric relationship parameters derived from step one to simulate the final DAM. The general flow chart was shown Figure 6. Both LAIe and canopy height data will be used to calibrate the modified SAFY-height model and calculate cost function. The cost function combines both the RMSE of crop LAIe and canopy height. The equation is shown as (Equation (7)): According to the results of Song's study [26], the cultivar-specific parameters were determined. The details are shown in Table 2.

Allometric Relationship between Plant Height and DAM
The DAM ranged from 15 to 1393 g/m 2 , and the plant height ranged from 18 to 90 cm for S1 from 8 May to 20 July. The piecewise linear regression model was well established from the plant height and DAM in S1. The relationship between canopy height and DAM was divide into two components, which is represented by two linear models. The turning point was the point where the data completely changes its behavior. Figure 7 shows the relationship. Figure 7. The relationship between winter wheat canopy height and dry aboveground biomass. Two phases of the relationship were determined from 8 May to 20 July in S1. The blue dots represent phase 1, which is the measurements from 8 May to 3 June; the orange dots represent phase 2, which is the measurements from 3 June to 20 July.
After the calibration of the modified SAFY-height model using the ground measured crop height and DAM in S1, the parameters a 1 , a 2 , b 1 , b 2 , and T turn , in Equation (1) was determined. The maximum, minimum, mean, median, and standard deviation (STD) list is shown in Table 3. According to the parameters of a 1 , a 2 , b 1 , and b 2 , question 1 can be written as

Estimated Plant DAM Using SAFY Model in S2
The median value of parameters of a 1 , a 2 , b 1 , and b 2 were adopted in the modified SAFY-height model to fix the relationship between winter wheat crop height and DAM. The final DAM in S2 was derived from the modified SAFY-height model using the cultivarspecific parameters of soft red winter wheat [26], height-DAM parameters, the UAV-based crop height maps, and UAV-based LAIe maps. Figure 8 shows final DAM maps using the modified SAFY-height model. Yield in the final DAM map ranged from 538 to 1808 g/m 2 . Figure 9 shows the final crop yield map derived from the SAFY model and the modified SAFY-height model which has same resolution as the harvester yield map. The yield ranged from 303 to 904 g/m 2 , and the STD and mean were 81.94 g/m 2 and 549 g/m 2 . Figure 10 shows the absolute difference yield maps for both SAFY and modified SAFY-height models compared with the harvester yield map. The comparison of mean, STD, coefficient of variation (CV), RMSE, and relative root mean square error (RRMSE) are shown in Table 4.      Figure 7 shows the piecewise linear relationship between the winter wheat canopy height and DAM in S1. It has a similar agreement to the results of Bakhshandeh et al. (2012) [14] where the piecewise linear regression model can model the relationship between the plant height and DAM. In this study, the winter wheat has different relationships between plant height and DAM during the tillering stages to the booting stages (phase one) and during booting stage to the heading stages (phase two). The rate of DAM changes with the plant height exhibits a significant decrease after the booting stage because the winter wheat started to produce the fruit after the booting stages, and plant height changes become minimal.

The Allometric Relationship between Winter Wheat Plant Height and DAM
The modified SAFY-height model using the ground measured winter wheat plant height and DAM in S1 was determined equation 1. After adopting the median values of the estimated parameters, Equation (1) can be written as Equation (8). The parameters a 1 , a 2 , b 1 , b 2 have similar values as the coefficients of the piecewise linear regression model for S1. The turning points of the total accumulated temperature was determined from the modified SAFY-height calibration, which was 712 • C on 3 June. The phenology of the winter wheat on 3 June was the beginning of the booting stages. This turning point of accumulated temperature indicated that the relationship between winter wheat plant height and DAM changes after the booting stages. Based on Equation (8), the RMSE between estimated DAM and ground measured DAM was estimated in S1. The first phase had a lower RMSE, 64 g/m 2 when the total accumulated temperature was less than T turn . The lower RMSE reveals that the plant height contributed more to the estimation of accumulated DAM in the early growth stages of winter wheat. The addition of plant height of winter wheat improved DAM estimation before the booting stage.
In contrast, the second phase of the relationship between plant height and DAM exhibited a larger RMSE, 197 g/m 2 . The larger RMSE revealed that the DAM derived from plant height has a lower accuracy after the booting stages. This is because the plant height has a lower variation in the later growth stages of winter wheat. Still, the accumulated DAM has a more significant variation due to the LAI and soil stress differences. The most contribution of the DAM comes from other factors rather than the plant height. Therefore, the plant height may not be an excellent factor to indicate the DAM after the booting stages. Although the piecewise linear regression model represents the allometric relationship between wither wheat canopy height and DAMwell (Bakhshandeh et al., 2010) [14], this model is incapable of demonstrating the yield variability of winter wheat after the booting stage due to the larger RMSE.

Uncertainty of Final Yield Estimation Using the Modified SAFY-Height Model
After calibrating the modified SAFY-height model, the final winter wheat DAM and yield map were generated in S2 (Figure 9) using the UAV-based photogrammetric LAI and height maps. A similar pattern of the final yield map was observed in this study by comparing the final yield map derived from the modified SAFY-height model with the harvester measured and SAFY model final yield map. The absolute difference maps were calculated for both the SAFY and modified SAFY-height models ( Figure 10). It shows that the upper left area in the field (blue rectangle area) has significant differences between these two absolute yield maps. The modified SAFY-height model significantly eliminated the estimation difference in the absolute yield map. However, the modified SAFY-height model still has a larger difference than the SAFY model in other parts. The total number of points with a difference (>150 g/m 2 ) is 206 out of 1821 points. Among these 206 points, 149 points have underestimated yield than the harvester measured yield. Most underestimated yield zones were located at the area with a lower LAI value, which is less than 1.01 m 2 /m 2 , on the 27 May LAI map ( Figure 11). The estimation error of the canopy height derived from the UAV-based photogrammetric point cloud should be considered since the height was an input for the modified SAFY-model. The RMSE of all three canopy maps increased as the crop grew, leading to a significant variation in later canopy height estimation. In addition, Song and Wang, (2019) [27] pointed out that the moving cuboid filter has an error if the threshold of the filter exceeds a specific range. Therefore, introducing the allometric relationship to the SAFY model also brings new errors for final products.
After removing the pixels with LAI value less than 1.01 on the 27 May LAIe map, the RMSE and RRMSE were calculated for both SAFY and the modified SAFY model. The modified SAFY-height model has lower RMSE and RRMSE, 91 g/m 2 and 16.41%, respectively, which is improved 3.3% and 0.5% from the SAFY model. In addition, the modified SAFY-height model has a similar standard deviation and CV with the harvested data. The details shown in Figure 12. It shows that the modified SAFY-height model has a minor improvement in estimating the final yield in higher LAI areas.

Future Study
Although plant height and DAM have a well-estimated allometric relationship for the winter wheat growing season, specifically in the early growth stages, the allometric relationship does not improve the final yield estimation when using the modified SAFYheight model. The lack of the relationship between plant height and LAI leads to a larger RMSE and an underestimation in the final yield estimation for the entire field. Since the new model improves the accuracy of yield estimation in high-LAI areas, the relationship between LAI and crop height can be considered in future studies to improve yield estimation in low-LAI areas. In addition, future research should consider external factors, such as soil moisture and rainfall, to improve the yield estimation accuracy. Due to the high correlation between plant height and DAM, the plant height can be used in many other machine learning methods such as random forest and Support Vector Machine, which interacts with different parameters and improves the yield estimation accuracy. Therefore, the machine learning approaches could be a potential research direction on final crop DAM and yield estimation.

Conclusions
This study developed the modified SAFY-height model by adding the relationship between plant height and DAM. The modified SAFY-height model uses both plant height and LAI information to calibrate the model and simulate the daily GLAI, plant height, and DAM. The allometric relationship between plant height and DAM was determined using the modified SAFY-height model in S1. In addition, the turning point of the accumulated temperature was determined using the modified SAFY-height model. The allometric relationship analysis reveals that the plant height is possible to estimate DAM before the turning temperature. Still, the plant height is not a good factor in DAM estimation after this turning point. The UAV-based plant height data before the booting stage of winter wheat, for the first time, was applied to the modified SAFY-height model to generate the final DAM and yield maps. The modified SAFY-height model has a minor accuracy improvement in estimation of the final yield in a high LAI area in the field.