Phenotyping of Plant Biomass and Performance Traits Using Remote Sensing Techniques in Pea (Pisum sativum, L.)

Field pea cultivars are constantly improved through breeding programs to enhance biotic and abiotic stress tolerance and increase seed yield potential. In pea breeding, the Above Ground Biomass (AGBM) is assessed due to its influence on seed yield, canopy closure, and weed suppression. It is also the primary yield component for peas used as a cover crop and/or grazing. Measuring AGBM is destructive and labor-intensive process. Sensor-based phenotyping of such traits can greatly enhance crop breeding efficiency. In this research, high resolution RGB and multispectral images acquired with unmanned aerial systems were used to assess phenotypes in spring and winter pea breeding plots. The Green Red Vegetation Index (GRVI), Normalized Difference Vegetation Index (NDVI), Normalized Difference Red Edge Index (NDRE), plot volume, canopy height, and canopy coverage were extracted from RGB and multispectral information at five imaging times (between 365 to 1948 accumulated degree days/ADD after 1 May) in four winter field pea experiments and at three imaging times (between 1231 to 1648 ADD) in one spring field pea experiment. The image features were compared to ground-truth data including AGBM, lodging, leaf type, days to 50% flowering, days to physiological maturity, number of the first reproductive node, and seed yield. In two of the winter pea experiments, a strong correlation between image features and seed yield was observed at 1268 ADD (flowering). An increase in correlation between image features with the phenological traits such as days to 50% flowering and days to physiological maturity was observed at about 1725 ADD in these winter pea experiments. In the spring pea experiment, the plot volume estimated from images was highly correlated with ground truth canopy height (r = 0.83) at 1231 ADD. In two other winter pea experiments and the spring pea experiment, the GRVI and NDVI features were significantly correlated with AGBM at flowering. When selected image features were used to develop a least absolute shrinkage and selection operator model for AGBM estimation, the correlation coefficient between the actual and predicted AGBM was 0.60 and 0.84 in the winter and spring pea experiments, respectively. A SPOT-6 satellite image (1.5 m resolution) was also evaluated for its applicability to assess biomass and seed yield. The image features extracted from satellite imagery showed significant correlation with seed yield in two winter field pea experiments, however, the trend was not consistent. In summary, the study supports the potential of using unmanned aerial system-based imaging techniques to estimate biomass and crop performance in pea breeding programs.


UAS Data Collection
In the winter and spring pea experiments, a total of 10 GCPs were uniformly distributed over each experimental area, including the field edges to minimize the planimetry error [38]. A marker stake was placed at each GCP location and remained in place throughout the season. Prior to each flight, boards (0.8 m x 0.5 m) that could be seen in the resulting UAS images were placed at each GCP position. The coordinates of each point were recorded at the end of the experiment with a RTK system based on SPS850 Global Navigation Satellite System receivers from Trimble Inc. (California, USA), which integrates a 450-900 MHz transmitter/receiver radio and a 72-channel L1/L2/L2C/L5/GLONASS GPS receiver.
RGB data was collected with a DJI-Phantom 4 Pro (Shenzhen, China) using its original 20 MP resolution, 25.4 mm CMOS camera with lens characteristics of 84° field of view and 8.8 mm/24 mm (35 mm format equivalent). DJI-phantom 4 Pro is powered with 6000 mAh LiPo 2S battery and the speed The spring pea field was located in the Plant Materials Center of the USDA, Washington, USA (46 • 43 12.83 N; 117 • 8 33.88 W). The plant materials in this experiment were the USDA Pea Single Plant Derived Core Collection (PSP), a genome wide association mapping population that has been previously phenotyped and genotyped [37]. The 307 accessions were planted in a randomized complete block design with three replications. This experiment was planted on 14 May 2018, plots consisted of two, 1.2 m long rows (Figure 1b). Once the plants reached 50% flowering, CH GT , lodging, and leaf type were measured, and AGBM (as total dry weight) was assessed through destructive sampling of Sensors 2019, 19, 2031 4 of 23 the entire plot. Lodging was measured as the ratio of the height of the canopy divided by the total length of the plant, i.e., the closer the ratio is to 1.0, the more erect (less lodged) the plants are. Data collection occurred on 28 June (1231 ADD), 5 July (1424 ADD), and 12 July (1648 ADD). ADD was calculated from the planting date.

UAS Data Collection
In the winter and spring pea experiments, a total of 10 GCPs were uniformly distributed over each experimental area, including the field edges to minimize the planimetry error [38]. A marker stake was placed at each GCP location and remained in place throughout the season. Prior to each flight, boards (0.8 m × 0.5 m) that could be seen in the resulting UAS images were placed at each GCP position. The coordinates of each point were recorded at the end of the experiment with a RTK system based on SPS850 Global Navigation Satellite System receivers from Trimble Inc. (California, USA), which integrates a 450-900 MHz transmitter/receiver radio and a 72-channel L1/L2/L2C/L5/GLONASS GPS receiver.
RGB data was collected with a DJI-Phantom 4 Pro (Shenzhen, China) using its original 20 MP resolution, 25.4 mm CMOS camera with lens characteristics of 84 • field of view and 8.8 mm/24 mm (35 mm format equivalent). DJI-phantom 4 Pro is powered with 6000 mAh LiPo 2S battery and the speed during data acquisition was 2 m/s; it works with the Global Navigation Satellite System (GNSS: GPS and GLONASS constellations) with average horizontal and vertical accuracies of~0.5 m and 1.5 m, respectively. The high-density data were collected in a double grid pattern with 90% overlap (both directions) at 20 m AGL (0.005 m of ground sample distance/GSD) to generate high accuracy digital surface models. As high-density images were collected from different angles (more points of view for each object on the field), it was expected that the process would improve the quality of the 3D reconstruction. The multispectral information was captured using a Double 4K camera (Sentera LLC, Minneapolis, USA) of 59 × 40.9 × 44.5 mm dimensions with 12.3 MP (0.005 m GSD) resolution of five spectral bands. The central wavelength and full-width half maximum data for R, G, B, red edge (RE), and NIR spectral bands were 650 nm and 64 nm, 548 nm and 44 nm, 446 nm and 52 nm, 720 nm and 39 nm, and 839 and 20 nm, respectively. This sensor was mounted on an ATI-AgBOT TM (ATI LLC., Oregon, USA) quadcopter with 1012 400 kv motor and dual 6000 mAh batteries; its positioning system is 3DR uBlox GPS (UAV Systems International, Las Vegas, USA) that works with a 3 V lithium rechargeable battery at 5 Hz update rate and a low noise regulator of 3.3 V. The multispectral data were collected in a single grid pattern with 80% frontal overlap and 70% side overlap, also at 20 m AGL. A white reference panel (0.25 m × 0.25 m; Spectralon Reflectance Target, CSTM-SRT-99-120) (Spectra Vista Cooperation, New York, USA) was placed on the field for radiometric correction during image processing.

Satellite Data Acquisition
A multispectral 1.5 m-GSD SPOT-6 satellite image was acquired from AIRBUS Defense & Space (Leiden, The Netherlands). The image captured on 3 June 2018 (close to 784 ADD) is composed by four spectral bands with the following range: R (625-695 nm), G (530-590 nm), B (455-525 nm), and NIR (760-890 nm). The original image was atmospherically corrected, but not geo-referenced. The satellite information was not used for spring pea plot evaluation for three reasons: (1) at 1.50 m GSD, it was not possible to differentiate between plots (~1.20 × 0.30 m), (2) at the time of the data capture, the plants in the spring pea trials were in early growth stages and small, and finally (3) alternative satellite images matching the UAS data collection dates were unavailable.

UAS-Based Imagery Analysis
Pix4D TM software was used to create the mosaics and DSM from both sensors (RGB and multispectral) through the 3D map template. During the stitching process, each RTK-GCP was fixed by identifying its position with 10 to 15 checkpoints representing GCP location on individual images (both fields and all data points). For the winter pea experiments, 5 RGB, 5 multispectral and 5 DSM mosaics were generated; while 3 RGB, 3 multispectral and 3 DSM mosaics were generated for the spring pea experiment. The white reference panel (99% reflectance in RGB-RE-NIR spectral range) imaged during each data collection was used to correct the image pixels in each band. Following this, using the "Array" command in AutoCAD (version 2018), the polygons representing each winter pea plot were digitized in a *.dxf format and further translated into * .shp. As the spring pea plots did not present a uniform grid pattern, they were directly digitized in * .shp format using Quantum GIS (QGIS, version 2.18.22). Each plot was labeled with plot ID based on experimental details.
The green-red vegetation index, normalized difference vegetation index, and normalized difference red edge index were computed using the following equations.
where R, G, RE, and NIR represents the reflectance in the red, green, red edge, and near infrared bands. The DSM (in m above the mean sea level) was obtained from the stitched image data. To extract the CSM, with the canopy height (in m AGL) information, the DTM was created based on the interpolation of elevation data over bare soil points, and subtracted from the DSM (Equation (4)).
Using data from the winter pea field plots at 1268 ADD as reference, an assessment of the quality of the RTK geo-rectification was performed by estimating the vertical position error (VPE) and the horizontal position error (HPE) [39,40] of the rectified and non-rectified mosaic images from the two sensors (RGB and multispectral). The HPE was calculated using Equation (5). The VPE was estimated as the sum of the changes in elevation among adjacent points calculated with the non-rectified image (∆ZNR) subtracted from those calculated with the rectified image (∆ZR) (Equation (6)).
where HPE and VPE are horizontal and vertical position errors, EE and NE are East and North direction errors, ∆ZNR and ∆ZR are elevation differences from non-rectified and rectified images, and total number of samples (n) is 4 ( Figure 2). The ∆Z is the sum of absolute difference in the elevation between two contiguous points (∆Z 1 +∆Z 2 +∆Z 3 +∆Z 4 , Figure 2). From the CSM, the UAS-based CH (CH UAS ), Canopy Coverage (CC) and Plot Volume (PV) were estimated. The CSM was segmented into two categories where pixels above 0.15 m AGL were classified as "canopy", and pixels below 0.15 m AGL were classified as "non-target canopy" to eliminate weeds and other noises from the crop of interest. The 0.15 m was set as empirical threshold selected manually based on observations. The count of "canopy" pixels of a single plot was multiplied by the pixel area (e.g., 25 × 10 −6 m 2 ) to get the CC (m 2 ). The PV (m 3 ) was computed by multiplying the CH UAS with CC. The binary image (non-canopy and canopy) was also used as a soil mask image.  From the CSM, the UAS-based CH (CHUAS), Canopy Coverage (CC) and Plot Volume (PV) were estimated. The CSM was segmented into two categories where pixels above 0.15 m AGL were classified as "canopy", and pixels below 0.15 m AGL were classified as "non-target canopy" to eliminate weeds and other noises from the crop of interest. The 0.15 m was set as empirical threshold selected manually based on observations. The count of "canopy" pixels of a single plot was multiplied by the pixel area (e.g., 25 × 10 -6 m 2 ) to get the CC (m 2 ). The PV (m 3 ) was computed by multiplying the CHUAS with CC. The binary image (non-canopy and canopy) was also used as a soil mask image.
With the "Zonal Statistics" plugin in QGIS, the mean (as relative vigor) and sum (as absolute vigor) statistics of the three VIs, CHUAS, CC, and PV were extracted and recorded in the attribute table of the plot polygons, where each plot was differentiated based on its specific ID. In order to verify the consistency of the data across time, the three VIs and the CHUAS were plotted as a function of the ADD and compared with a reference dry matter curve [41].
In addition to the features specified above, lodging assessment was performed in spring pea. The changes in the CHUAS between first and second data points, and between the first and third data points were employed to calculate the lodging in spring pea. When a plot lodges, not only does the CH decrease, but the CC increases, due to an increase in surface area. For these reasons, both features were utilized during lodging estimation. For the lodging estimation between data points 1 and 3, the difference in absolute CC values was multiplied with the differences between CHUAS data (Equation (7)).
where 1 and 3 represent data collected at time points 1 and 3, 1231 ADD and 1648 ADD, respectively.
Green band (from RGB orthomosaic) frequencies were plotted for the two leaf types in the spring and winter peas. Additionally, the mean and the standard deviation of the green reflectance were also computed as indicators of greenness and its variability. This processing was carried with the multispectral mosaic collect at 1231 ADD in the spring pea plots and 1268 ADD in winter pea plots.

Satellite-Based Imagery Analysis
Using the "Georeferencer" tool in QGIS, the satellite image was rectified to the correct location. First, based on satellite archive Bing imagery (Bing aerial with layers) displayed with the "Open Layers" plugin, the original image was geo-located to its respective region with an error that would oscillate between 1-5 m. Second, the specific location of the winter pea experimental field was corrected to a sub-meter accuracy using the UAS RTK-mosaics as reference by matching the corner points of the field. In order to increase the resolution of the multispectral data from 6.0 m GSD to 1.5 m GSD, a pan-sharpening processing, based on a higher resolution panchromatic band, was performed in Erdas Imagine (version 14.1, Hexagon Geospatial) using the high pass filtering algorithm, which presented the clearest contrast between soil and vegetation pixels, compared with other methods like principal component analysis, hyperspectral color sharpening, and Brovey transform. GRVI and NDVI were computed with the satellite image following Equation 1 and With the "Zonal Statistics" plugin in QGIS, the mean (as relative vigor) and sum (as absolute vigor) statistics of the three VIs, CH UAS , CC, and PV were extracted and recorded in the attribute table of the plot polygons, where each plot was differentiated based on its specific ID. In order to verify the consistency of the data across time, the three VIs and the CH UAS were plotted as a function of the ADD and compared with a reference dry matter curve [41].
In addition to the features specified above, lodging assessment was performed in spring pea. The changes in the CH UAS between first and second data points, and between the first and third data points were employed to calculate the lodging in spring pea. When a plot lodges, not only does the CH decrease, but the CC increases, due to an increase in surface area. For these reasons, both features were utilized during lodging estimation. For the lodging estimation between data points 1 and 3, the difference in absolute CC values was multiplied with the differences between CH UAS data (Equation (7)).
where 1 and 3 represent data collected at time points 1 and 3, 1231 ADD and 1648 ADD, respectively. Green band (from RGB orthomosaic) frequencies were plotted for the two leaf types in the spring and winter peas. Additionally, the mean and the standard deviation of the green reflectance were also computed as indicators of greenness and its variability. This processing was carried with the multispectral mosaic collect at 1231 ADD in the spring pea plots and 1268 ADD in winter pea plots.

Satellite-Based Imagery Analysis
Using the "Georeferencer" tool in QGIS, the satellite image was rectified to the correct location. First, based on satellite archive Bing imagery (Bing aerial with layers) displayed with the "Open Layers" plugin, the original image was geo-located to its respective region with an error that would oscillate between 1-5 m. Second, the specific location of the winter pea experimental field was corrected to a sub-meter accuracy using the UAS RTK-mosaics as reference by matching the corner points of the field. In order to increase the resolution of the multispectral data from 6.0 m GSD to 1.5 m GSD, a pan-sharpening processing, based on a higher resolution panchromatic band, was performed in Erdas Imagine (version 14.1, Hexagon Geospatial) using the high pass filtering algorithm, which presented the clearest contrast between soil and vegetation pixels, compared with other methods like principal component analysis, hyperspectral color sharpening, and Brovey transform. GRVI and NDVI were computed with the satellite image following Equations (1) and (2). The mean and sum statistics were extracted from the plot polygons layer created for low altitude satellite imagery.

Statistical Analysis
Pearson's correlation matrix between the ground truth and UAS-based data, averaged by entry, was calculated in RStudio (Version 1.1.423). For spring peas, the correlations were calculated using plot-by-plot comparisons, since the replicates of the same entry were not always harvested on the same day. Using RStudio, the least absolute shrinkage and selection operator algorithm [25] was employed to predict AGBM using the well-correlated image features. To assess the Lasso prediction accuracy, using 85% of the original dataset, a cross-validation of the mean absolute error and the correlation coefficient (r) between the estimated and actual values were computed. The features were centered and scaled using the 'preProcess' function for comparable coefficient generation, where mean data was subtracted from each value of vegetation indices and divided with standard deviation. With the winter pea data, the resulting coefficients from Lasso were used to estimate the AGBM in experiments 1822 and 1823 where ground truth AGBM data were not available.

UAS-Based Position Data Accuracy
The horizontal and vertical dilutions of precision (HDOP and VDOP) representing the GCP coordinates reading accuracy are shown in Table 1. The horizontal and vertical position errors using RGB and multispectral mosaic images without RTK rectification were higher than those with RTK rectification ( Table 2). After image rectification, the correlation coefficient between CH GT and CH UAS increased for both RGB-and multispectral-CSMs (Table 3). The differences between CH GT and CH UAS can be attributed to human error during ground truth data collection, and some variances in the grid pattern and overlap percentages, since the resolution and the flight altitudes of sensors were similar. Despite observing higher accuracy and correlation between CH GT and CH UAS , the use of a double grid pattern and high overlap percentage may not be necessary to monitor research plots with simple geometry, such as in evaluated winter and spring pea breeding experiments. This will in-turn save battery life (thus increasing flight time and efficiency), data storage space, and image processing time. The use of the double grid pattern and higher overlap percentages with RTK-GPS rectification are necessary for monitoring and 3D mapping of more complex crop geometry, such as plant architecture with thin and narrow canopies (e.g., apple orchards and grape vineyards). Furthermore, the centimeter to sub-centimeter accuracy in the horizontal and vertical positions obtained with the RGB mosaic suggest its functionality to generate accurate plot length data, which is an important trait frequently monitored in breeding programs to estimate yield per unit area. Table 3. Correlation between CH GT and CH UAS for the CSMs obtained from RGB and multispectral imageries (with its respective flight parameters differences) with and without RTK rectification, based on 1268 ADD time point data. All correlation coefficients were significant at P < 0.001.

Mosaic
Flight

Winter Pea Growth and Development
The average vegetation index and plot volume data across different time periods in winter pea showed a similar pattern as the reference dry matter curve [41]. At the beginning of the season, the three VIs values were about 0.30 units. At 784 ADD, the GRVI and NDVI indices increased, while the CH and PV averaged approximately 0.24 m and 0.20 m 3 , respectively (Figures 3 and 4). The GRVI, CH UAS , and PV continued to increase until 1268 ADD; however, the NDVI changed marginally, which could be due to saturation. After 1268 ADD (flowering), the VIs decreased as the plants approached physiological maturity and senescence. Similarly, the CH UAS and PV also decreased at the end of the season because of maturity and crop lodging. While a similar pattern was observed with NDVI and GRVI data, the NDRE data were low, which could be due to less abiotic stress in winter pea experiments during the season.
The data validates the generic vegetative growth stages, where the crop canopy vigor increases and reaches maximum photosynthetic activity at flowering, resulting in a higher NIR reflection and absorption in visible wavelengths, which can be observed from the increase in VI values from 365 ADD to 1268 ADD. During the seed development and pod filling stage that represents the translocation of photo-assimilates to the seeds after flowering, there is a decrease in leaf biomass accumulation, which can be observed with a decrease of VI values.

Correlation between Image Features and Performance Traits
In couple winter pea experiments, a strong correlation between image features (GRVI, NDVI, NDRE, CHUAS, CC, PV) and seed yield was observed, especially at 1268 ADD. At 1268 ADD, most image features were also correlated with FN (Table 4). In experiment 1823, at 1268 ADD, high correlation coefficients between image features and F50, PM, and SY were observed. In experiment 1822, high correlations between image features and SY were observed starting earlier in the season. In general, imaging between 1268 ADD (flowering) and 1725 ADD (pod development) is recommended for capturing yield differences. There were no significant differences between sum and mean vegetation index values.

Correlation between Image Features and Performance Traits
In couple winter pea experiments, a strong correlation between image features (GRVI, NDVI, NDRE, CH UAS , CC, PV) and seed yield was observed, especially at 1268 ADD. At 1268 ADD, most image features were also correlated with FN (Table 4). In experiment 1823, at 1268 ADD, high correlation coefficients between image features and F50, PM, and SY were observed. In experiment 1822, high correlations between image features and SY were observed starting earlier in the season. In general, imaging between 1268 ADD (flowering) and 1725 ADD (pod development) is recommended for capturing yield differences. There were no significant differences between sum and mean vegetation index values.
The lower correlations found at early stages (365 ADD) could be attributed to the distortions caused by the brightness of bare soil within the plots [42]. To overcome this limitation, Badgley et al. [43] proposed the NIR vegetation reflectance indicator (NIRv), where 0.08 is subtracted from the product of the total NIR and NDVI, which represents the proportion of the digital number of a pixel attributed to the vegetation. In the present study, when NIRv was used, an increase in the correlations between NIRv with SY and AGBM at 365 ADD (Table 5) was observed, which did not affect relationships at 1268 ADD (high canopy cover). This suggests the importance of NIRv usage as an indicator for UAS-based ABGM predictions under high soil exposure environments or early growth stages.  According to the Reference [41], the high yielding genotypes have a larger leaf area index that leads to accumulation of more biomass during flowering. This is validated by the strong correlations found between leaf area index related image features such as CC and yield at 1268 ADD. Furthermore, an increase in the correlations with the phenological traits was detected in the winter pea experiments at 1725 ADD during pod development and maturity. These stronger correlations are attributed to the contrasting canopies characteristics among early and late F50 and PM plots that were easily captured with remote sensing data ( Figure 5). product of the total NIR and NDVI, which represents the proportion of the digital number of a pixel attributed to the vegetation. In the present study, when NIRv was used, an increase in the correlations between NIRv with SY and AGBM at 365 ADD (Table 5) was observed, which did not affect relationships at 1268 ADD (high canopy cover). This suggests the importance of NIRv usage as an indicator for UAS-based ABGM predictions under high soil exposure environments or early growth stages. According to the Reference [41], the high yielding genotypes have a larger leaf area index that leads to accumulation of more biomass during flowering. This is validated by the strong correlations found between leaf area index related image features such as CC and yield at 1268 ADD. Furthermore, an increase in the correlations with the phenological traits was detected in the winter pea experiments at 1725 ADD during pod development and maturity. These stronger correlations are attributed to the contrasting canopies characteristics among early and late F50 and PM plots that were easily captured with remote sensing data ( Figure 5).  The FN is also related to the flowering time in the sense that lower FN [44] can be associated with earlier F50 and senescence in winter pea. This can explain the high correlations between FN and GRVI-mean, CH UAS , and PV at 1725 ADD and 1948 ADD, when image features from lower FN entries with early F50 and senescence could be differentiated from higher FN entries. Additionally, the impact of FN on yield can also be captured with VIs in winter pea, where earlier flowering entries may have more reproductive nodes per plant and therefore higher seed yield ( Figure 6).
The FN is also related to the flowering time in the sense that lower FN [44] can be associated with earlier F50 and senescence in winter pea. This can explain the high correlations between FN and GRVI-mean, CHUAS, and PV at 1725 ADD and 1948 ADD, when image features from lower FN entries with early F50 and senescence could be differentiated from higher FN entries. Additionally, the impact of FN on yield can also be captured with VIs in winter pea, where earlier flowering entries may have more reproductive nodes per plant and therefore higher seed yield ( Figure 6). At 1231 ADD, in the spring pea experiment, PV was significantly correlated with CHGT ( Table  6). The CHUAS, GRVI, and NDVI were correlated with CHGT at 1231 and 1648 ADD. Estimations based on elevation data, such as CHUAS and PV demonstrated high correlation with the CHGT measurements in spring pea in most cases, which could be because these features were correlated with F50, FN, and SY as found in some winter pea experiments.

Correlation between Image Features and AGBM
In the winter pea experiments, strong correlation between GRVI-sum, NDVI-sum, and NDREsum with AGBM were found at 1268 ADD ( Table 7). The correlation was lower between AGBM with NDVI and NDRE mean values. The elevation-based features were not significantly correlated to AGBM. However, in the spring pea experiments, better correlations between image features and AGBM were found, especially at 1231 ADD. The variability in F50 among the accessions in the spring pea experiment could have contributed to higher correlation with remote sensing data ( Table 7). The At 1231 ADD, in the spring pea experiment, PV was significantly correlated with CH GT ( Table 6). The CH UAS , GRVI, and NDVI were correlated with CH GT at 1231 and 1648 ADD. Estimations based on elevation data, such as CH UAS and PV demonstrated high correlation with the CH GT measurements in spring pea in most cases, which could be because these features were correlated with F50, FN, and SY as found in some winter pea experiments.

Correlation between Image Features and AGBM
In the winter pea experiments, strong correlation between GRVI-sum, NDVI-sum, and NDRE-sum with AGBM were found at 1268 ADD ( Table 7). The correlation was lower between AGBM with NDVI and NDRE mean values. The elevation-based features were not significantly correlated to AGBM. However, in the spring pea experiments, better correlations between image features and AGBM were found, especially at 1231 ADD. The variability in F50 among the accessions in the spring pea experiment could have contributed to higher correlation with remote sensing data ( Table 7). The major finding from the spring pea dataset was that the PV extracted from images was consistently correlated with AGBM, which could be useful in breeding programs.

AGBM Prediction with Model Development
With the winter pea data, the Lasso method was implemented using highly correlated image feature data at 1268 ADD to predict AGBM. The results from this model showed a R 2 of 0.99 at F < 0.001 significance with four features. The resulting equation for AGBM estimation is defined as: where a 1 to c 1 represent the coefficients generated by Lasso for GRVI-sum  (8)) was used to estimate the AGBM in the complete data set from experiments 1821 and 1821cc at 1268 ADD. The correlation coefficient between estimated and actual AGBM was 0.60 (P < 0.001) (Figure 7), with a mean absolute error of 2.82 kg.
With the spring pea data, the Lasso method was implemented using the information from PV, GRVI (sum), NDRE (sum), and NDVI (sum and mean) at 1231 ADD as the best correlated scenarios. The results from the model showed a R 2 of 0.74 at F < 0.001 significance with five features. The equation for AGBM estimation is defined as (Equation (9)): where a 2 to d 2 represent the coefficients for GRVI-sum (a 2 = 7.80), NDVI-sum (b 2 = 2.46), NDVI-mean (c 2 = 4.98), NDRE-sum (d 2 = 4.98) and PV (e 2 = 10.73), and the intercept (f 2 =87.83) of the function. During validation with the complete data set at 1231 ADD, the correlation coefficient between the estimated and the actual AGBM was 0.84 (P < 0.001) (Figure 7) with a mean absolute error of 19.28 g. major finding from the spring pea dataset was that the PV extracted from images was consistently correlated with AGBM, which could be useful in breeding programs.

AGBM Prediction with Model Development
With the winter pea data, the Lasso method was implemented using highly correlated image feature data at 1268 ADD to predict AGBM. The results from this model showed a R 2 of 0.99 at F < 0.001 significance with four features. The resulting equation for AGBM estimation is defined as: * * * * /1.5E⁴ (8) where a1 to c1 represent the coefficients generated by Lasso for GRVI-sum (a1 = 1.56), NDVI-sum (b1 = −0.83), NDRE-sum (c1 = 0.75), CC (d1 = −0.01), and the intercept (e1 = 8.85) of the function. To cross validate, the equation (Equation (8)) was used to estimate the AGBM in the complete data set from experiments 1821 and 1821cc at 1268 ADD. The correlation coefficient between estimated and actual AGBM was 0.60 (P < 0.001) (Figure 7), with a mean absolute error of 2.82 kg.
With the spring pea data, the Lasso method was implemented using the information from PV, GRVI (sum), NDRE (sum), and NDVI (sum and mean) at 1231 ADD as the best correlated scenarios. The results from the model showed a R 2 of 0.74 at F < 0.001 significance with five features. The equation for AGBM estimation is defined as (Equation 9): * * * * * /1.0Eᶟ (9) where a₂ to d₂ represent the coefficients for GRVI-sum (a₂ = 7.80), NDVI-sum (b₂ = 2.46), NDVI-mean (c₂ = 4.98), NDRE-sum (d₂= 4.98) and PV (e₂ = 10.73), and the intercept (f₂ =87.83) of the function. During validation with the complete data set at 1231 ADD, the correlation coefficient between the estimated and the actual AGBM was 0.84 (P < 0.001) ( Figure 7) with a mean absolute error of 19.28 g. In both pea types, sum data of GRVI, NDVI, and NDRE were variables that made a strong contribution to the AGBM estimations; alongside PV and mean NDVI in spring pea and CC in winter pea, that had an impact on the predictions. Equation (8) was used to estimate the winter pea AGBM based on 1268 ADD data. The estimated AGBM values for the 20 entries are plotted in Figure 8. The entries with the lower AGBM estimations had low canopy cover and low vegetation index values. In experiment 1822, entries 6 to 9 and 20 had AGBM estimations above the average; while, entries 12 and 14 to 18 were clearly below the average. However, in experiment 1823, the AGBM estimated for the majority of the entries were close to or below average, except for entries 1, 3, 4, 12, 15, 17, and 18. The higher AGBM accumulation entries predicted in field pea experiments shared a lower FN, F50, and PM, and higher SY. The results presented in this study need to be integrated with multiple season data in order to create a larger data-pool to build a robust machine learning prediction method. In both pea types, sum data of GRVI, NDVI, and NDRE were variables that made a strong contribution to the AGBM estimations; alongside PV and mean NDVI in spring pea and CC in winter pea, that had an impact on the predictions. Equation 8 was used to estimate the winter pea AGBM based on 1268 ADD data. The estimated AGBM values for the 20 entries are plotted in Figure 8. The entries with the lower AGBM estimations had low canopy cover and low vegetation index values. In experiment 1822, entries 6 to 9 and 20 had AGBM estimations above the average; while, entries 12 and 14 to 18 were clearly below the average. However, in experiment 1823, the AGBM estimated for the majority of the entries were close to or below average, except for entries 1, 3, 4, 12, 15, 17, and 18. The higher AGBM accumulation entries predicted in field pea experiments shared a lower FN, F50, and PM, and higher SY. The results presented in this study need to be integrated with multiple season data in order to create a larger data-pool to build a robust machine learning prediction method.

Leaf Type Characterization
The GRVI, and average and standard deviation of green bands from spring and winter pea field plots are presented in Figure 9. The average green reflectance data in semi-leafless entries were higher in the winter pea plots than spring pea plots, while GRVI values showed an opposite pattern. In winter pea, the variation in greenness was clearly different between experiments, which could be resulting from different pigmentation (anthocyanin), but not between leaf types. Higher variability was detected in the normal leaf type. Based on green band average, the leaf type was classified with 90%, 73%, and 87% accuracy in winter pea experiments 1821, 1821cc, and 1822, respectively, and with 74% accuracy in spring pea.

Leaf Type Characterization
The GRVI, and average and standard deviation of green bands from spring and winter pea field plots are presented in Figure 9. The average green reflectance data in semi-leafless entries were higher in the winter pea plots than spring pea plots, while GRVI values showed an opposite pattern. In winter pea, the variation in greenness was clearly different between experiments, which could be resulting from different pigmentation (anthocyanin), but not between leaf types. Higher variability was detected in the normal leaf type. Based on green band average, the leaf type was classified with 90%, 73%, and 87% accuracy in winter pea experiments 1821, 1821cc, and 1822, respectively, and with 74% accuracy in spring pea. Further investigation of leaf type characterization as well as pigmentation (chlorophyll, anthocyanin) and differences in anatomical and morphological structures is needed. The leaf type influences the total leaf area (hence biomass) and also strongly affects lodging tolerance. When the correlation analysis between GRVI-sum with AGBM was calculated by leaf type, the correlation coefficients increased (Table 8). Thus, it may be important to integrate leaf type classification with remote sensing data analysis for more robust variety selection. Table 8. Correlation coefficients (r) between GRVI-sum with AGBM by leaf type (af and Af) in spring peas at 1231 ADD. All correlations were significant at P < 0.001.

Lodging Estimation
Despite the small size of spring pea plants and the plots, the use of elevation data was promising for monitoring lodging. It is hypothesized that the relationships will grow stronger with larger plot sizes, where the changes in CH and CC between time points can be captured with ease. Lodging is influenced by stem strength as well as the weight of the developing pods, and can be defined as a change in the vertical height of plants producing its inclination [45] and can be estimated based on decrease in CH and increase in CC ( Figure 10) across time points. Lodging assessments based on the detection of changes in CH between dates 1 to 2 and 1 to 3, correlated with ground truth lodging observations with r of 0.58 and 0.57, respectively. Furthermore, including the absolute CC value for the lodging estimation between dates 1 and 3 (Equation (7)), the correlation coefficient increased up to 0.70. Further investigation of leaf type characterization as well as pigmentation (chlorophyll, anthocyanin) and differences in anatomical and morphological structures is needed. The leaf type influences the total leaf area (hence biomass) and also strongly affects lodging tolerance. When the correlation analysis between GRVI-sum with AGBM was calculated by leaf type, the correlation coefficients increased (Table 8). Thus, it may be important to integrate leaf type classification with remote sensing data analysis for more robust variety selection. Table 8. Correlation coefficients (r) between GRVI-sum with AGBM by leaf type (af and Af ) in spring peas at 1231 ADD. All correlations were significant at P < 0.001.

Lodging Estimation
Despite the small size of spring pea plants and the plots, the use of elevation data was promising for monitoring lodging. It is hypothesized that the relationships will grow stronger with larger plot sizes, where the changes in CH and CC between time points can be captured with ease. Lodging is influenced by stem strength as well as the weight of the developing pods, and can be defined as a change in the vertical height of plants producing its inclination [45] and can be estimated based on decrease in CH and increase in CC ( Figure 10) across time points. Lodging assessments based on the detection of changes in CH between dates 1 to 2 and 1 to 3, correlated with ground truth lodging observations with r of 0.58 and 0.57, respectively. Furthermore, including the absolute CC value for the lodging estimation between dates 1 and 3 (Equation (7)), the correlation coefficient increased up to 0.70.

Comparison with Satellite Data
The GRVI data extracted from the satellite images were not significantly correlated with AGBM and seed yield in experiments 1821 and 1821cc (Table 9). In experiments 1822 and 1823, the correlations between GRVI and NDVI mean, and SY were significant but lower than those using UAS data. The lower spatial resolution of the satellite image and spectral mixing (canopy and soil) on the field edges resulted in a poor relationship between image features and ground-reference data.

Comparison with Satellite Data
The GRVI data extracted from the satellite images were not significantly correlated with AGBM and seed yield in experiments 1821 and 1821cc (Table 9). In experiments 1822 and 1823, the correlations between GRVI and NDVI mean, and SY were significant but lower than those using UAS data. The lower spatial resolution of the satellite image and spectral mixing (canopy and soil) on the field edges resulted in a poor relationship between image features and ground-reference data. Table 9. Correlation coefficients (r) between UAS-based (1268 ADD) and satellite-based (887 ADD) GRVI and NDVI (sum and mean) with AGBM and SY in the winter pea experiments. The decrease in the satellite image resolution was not directly proportional to the decrease on its relationship with ground reference data. As an example, decreasing the resolution of an UAS image, by resampling its pixel size with the nearest neighbor method, to the level of the satellite image [46] SPOT 6 (1.50 m), resulted in an image in which the field edge was un-recognizable. At the same resolution, the satellite image provided more details of the field (Figure 11), which could be the reason for the significant correlations between satellite image features and seed yield obtained with the 1.5 m satellite image resolution. The decrease in the satellite image resolution was not directly proportional to the decrease on its relationship with ground reference data. As an example, decreasing the resolution of an UAS image, by resampling its pixel size with the nearest neighbor method, to the level of the satellite image [46] SPOT 6 (1.50 m), resulted in an image in which the field edge was un-recognizable. At the same resolution, the satellite image provided more details of the field (Figure 11), which could be the reason for the significant correlations between satellite image features and seed yield obtained with the 1.5 m satellite image resolution. In the future, satellite data is anticipated to have higher spectral and spatial resolution. The actively sensed data from orbital Synthetic Aperture Radar sensors offers new opportunities for plant phenotyping [47], because of its feasibility for crop phenology monitoring [48,49], crop height [50], and lodging estimations [51]. Furthermore, the anticipated launch of the FLuorescence EXplorer satellite mission, will provide sun-induced crop fluorescence spectral data that will create new research opportunities [52,53].

Conclusions
In this study, the potential of UAS-based imaging techniques to estimate biomass and crop performance in pea breeding programs was evaluated. In winter pea experiments, a strong correlation between all the image features and seed yield was observed at flowering; while at pod development and maturity, an increase in the correlations with phenological traits was detected. Spectral data was also found to be useful in leaf type identification. Overall, elevation based remote sensing data was highly correlated with CHGT and was also suitable for lodging assessment in spring pea; furthermore, in some winter pea experiments, this type of information was correlated with F50, FN, and SY. These results were obtained regardless of the use of flight plans with double grid pattern and high overlap percentage.
AGBM was found to be highly correlated with image features at 1268 and 1231 ADD (flowering) in the winter and spring peas, respectively. The Lasso model developed with selected image features was able to estimate AGBM with a high level of accuracy. The proposed methods and feature extraction can be used for evaluating biomass in forage breeding trials as well. The satellite imagery needs to be further explored for phenotyping applications.

ADD
accumulated degree days; Af normal leaf type; af semi-leafless leaf type; AGBM above ground biomass; AGL above ground level; CC canopy cover; In the future, satellite data is anticipated to have higher spectral and spatial resolution. The actively sensed data from orbital Synthetic Aperture Radar sensors offers new opportunities for plant phenotyping [47], because of its feasibility for crop phenology monitoring [48,49], crop height [50], and lodging estimations [51]. Furthermore, the anticipated launch of the FLuorescence EXplorer satellite mission, will provide sun-induced crop fluorescence spectral data that will create new research opportunities [52,53].

Conclusions
In this study, the potential of UAS-based imaging techniques to estimate biomass and crop performance in pea breeding programs was evaluated. In winter pea experiments, a strong correlation between all the image features and seed yield was observed at flowering; while at pod development and maturity, an increase in the correlations with phenological traits was detected. Spectral data was also found to be useful in leaf type identification. Overall, elevation based remote sensing data was highly correlated with CH GT and was also suitable for lodging assessment in spring pea; furthermore, in some winter pea experiments, this type of information was correlated with F50, FN, and SY. These results were obtained regardless of the use of flight plans with double grid pattern and high overlap percentage.
AGBM was found to be highly correlated with image features at 1268 and 1231 ADD (flowering) in the winter and spring peas, respectively. The Lasso model developed with selected image features was able to estimate AGBM with a high level of accuracy. The proposed methods and feature extraction can be used for evaluating biomass in forage breeding trials as well. The satellite imagery needs to be further explored for phenotyping applications.