We demonstrate our approach using an agricultural field under active commercial cultivation of soybeans located in the Arkansas Delta. For this study, we acquired time-series datasets of UAV and geophysical data during the 2018 growing season. This provided the opportunity to investigate the relationship between plant growth and spatially heterogeneous soil physical properties, as well as to analyse the temporal dynamics of such a relationship. All the methodology and data collection comply with relevant institutional and national legislation.
Site description and crop management
The investigated study area (34° 24.458′ N, 91° 40.462′ W), located in Humphrey, Arkansas, is a crop field of about 20 hectars that is under active commercial cultivation of soybeans (Fig. 1).
Study area located in Humphrey, Arkansas: (a) map produced by cropScape (USDA) showing part of the Delta region with the main agricultural land covers; (b) map of the field under study reporting the locations of the ground data collection; (c) map showing the topographic elevation derived by an UAV-based DEM, acquired acquisition prior to planting; (d) gSSurgo map provided by USDA, showing the soil types. Black vertical lines separate the east, center, and west sections used in the statistical analysis. Maps made with QGIS (v 3.6, https://www.qgis.org)36.
The study area is part of the Lower Mississippi River Basin’s (LMRB) alluvial plain37. This region is characterized by a wide range of cropland, such as soybean, corn, rice, grass, and pasture, supported by a humid subtropical climate with hot, humid summers and mild, slightly drier winters.
The crop field is adjacent to the oxbow Glenwood Lake (shown in Fig. 1) and is characterized by an alternating of Hebert silt-loam and Rilla silt-loam, according to the description provided by the gSSURGO map38. The Hebert profile consists of a very deep, poorly drained, moderately slowly permeable soil that formed in silty alluvium, while the Rilla profile consists of a very deep, well drained, moderately permeable soils that formed in reddish silty and loamy alluvium.
High-quality soybean seed requires three appropriate conditions for germination—soil moisture, temperature and oxygen. Provided that the soil is not saturated, oxygen concentration is not a limitation, but germination will not occur in flooded soil dueto lack of oxygen. Research shows that soybeans are in general very susceptible to soil saturation and anoxic conditions, which would provoke damage to the root system and in some case the termination of the plant growth39.
Planting started on April 21st in rows (raised beds) with a spacing of 0.9 m. Plants reached their first vegetative stage of emergence (VE) on the first days of May and developed their first node with fully developed unifoliate leaves around mid-May (V1). The stage of initial bloom, which represents the starting of the reproductive stage (R1–R8) after the last vegetative stage, started in mid-June, reaching the full maturity (R8) about the end of September. Subsequent to R8, the leaves dry then fall from the plant at the time of harvest, which was performed on October 24th.
A common channeled (furrowed) surface irrigation approach was used as an irrigation system, which consists of a plastic pipe (polytube) located along one side of the field to deliver water into the channels between the rows40. At the site, the plastic pipe was located on the east side of the field with the rows and furrows running east to west following the main slope gradient (Fig. 1) from east to west, in order to facilitate the water flow. Smoothing of the field’s ground surface was performed prior to the planting to ensure a good surface drainage. Irrigation started on June 12th, after the rain period, and performed every 14 days.
Temperature and precipitation were measured by a weather station that we installed in the southest corner of the field (Fig. 1). Based on the data collected, the average daily temperature increased from 14.6 to 24.7 °C between planting and the first vegetative stage (April–May), reached 27 °C during the vegetative period through the initial reproductive period (June–July), and decreased from 27 to 24 °C during the final period of maturity (August–September). The 2018 growing season was characterized by sporadic heavy-rain events, especially during the vegetative period (from April to mid-June), with a total cumulative precipitation (period April–September) of about 650 mm (Supplementary Figure S1). As a result of rain events, the soil moisture in April ranged from 0.17 (m3/m3) to 0.3 (m3/m3) (Supplementary Figure S2).
Instrumentation, data acquisition, and preprocessing
UAV-data acquisition
Multitemporal image acquisitions were performed over the field using UAV-mounted sensors. The images were acquired on: May 8th, a few days after planting, May 28th, during the early vegetative development, and June 25th, July 23rd, September 19th, the beginning, middle, and end of the reproductive period.
We used the DJI Matrice 600 flying platform with the DJI FC350 on-board RGB (red, green, blue) camera as the primary image acquisition sensor. A flight plan to acquire high-resolution images at centimeter resolution was designed using this acquisition system. Flights were performed during the maximum Sun’s elevation and with clear sky to avoid possible changes in light conditions, at an altitude of about 60 m, with a 75% image overlap, reaching a ground resolution of 2 cm.
In order to minimize the impact of possible variations in illumination conditions during each acquisition, two normalized mosaics were computed by using the pictures of a reference panel taken before and after each flight. The final mosaic was obtained by averaging the two normalized mosaics.
A total of 13 Ground Control Points (GCP) were positioned along the field perimeter and within the field for the mosaic reconstruction (Fig. 1). These consisted of flat targets and were surveyed at each acquisition with a Topcon real‐time kinematic differential global positioning system (Hyper-V RTK-GPS). The system is composed of a rover and a fixed base station, ensuring high-resolution geolocation information at centimeter resolution. All the collected geolocation data were post-processed obtaining an average planar error of about 2 cm and a vertical error of about 3 cm.
We used the commercial PhotoScan (Agisoft) software for photogrammetry to construct the final mosaics. The software uses photogrammetry techniques based on structure from motion from multi-view stereo (SfM-MVS) images42. The processing produced mosaics with spatial resolution less than 10 cm, but the mosaics used in this work have a have a pixel size of 10 cm. Additionally, we used the software to compute digital elevation models (DEMs) from the RGB images with a final vertical resolution less than 5 cm.
Electromagnetic Induction (EMI) system
Soil ECa was measured with an electromagnetic induction (EMI) system (CMD Mini-Explorer, GF Instruments), which consists of a transmitter and three receiver coils, allowing an effective acquisition of soil apparent electrical conductivity (ECa) averaged across a depth range of 0.5 m, 1.0 m, and 1.8 m from the surface. The instrument provides soil ECa measurements in milliSiemes per meter (mS/m) and the in-phase in parts per thousand (not considered in this study). The system was mounted in a sled pulled by an all-terrain vehicle and supported by a Differential Global Positioning System (DGPS) in order to tag each measurement with high-resolution positioning information. Time-lapse soil ECa data collection was performed four times, including before and during the growing season: April 21st, May 12th , June 11th, and July 9th. EMI acquisitons were planned according to the irrigation plan and UAV acquisitions. Since the EMI instrument is sensitive to changes in soil moisture, we avoided the collection during or right after irrigation events. The data were acquired consistently every 10 rows with a distance between traverses of about 9 m.
In situ soil monitoring, ground imagery, and soil sampling
Soil sensors (5TE, METER) were used for continuous acquisition of volumetric water content (VWC) and soil temperature. The sensors were deployed in five areas (A, B …, E, see Fig. 1) based on two factors: a) spatial varibaility of geophysical properties, b) differences in yield during the 2017 growing season (historical data available for this study). Each area has one logger that accomodates 4 sensors within a maximum radius of 4 m. To capture both shallower and deeper soil dynamics, we positioned 2 sensors vertically aligned at the depths of 12 cm and 25 cm on both the south (side1) and north (side2) of the logger. The areas A and D were positioned in high and low soil EC zones, respectively, while B was in an area of transition. Sensors A, B, and D were positioned along the same rows to forms two transects, one on the south and one on the north. Sensor E was positioned within a a very low 2017 yield in the west side of the field but at a mid-range of soil EC values. The logger C was positioned across a soil EC boundarie, but it stopped recording data at early season and was discarded from this analysis. The collected temprature data are used to correct the EMI acquisitions, whereas soil moisture is used to track the irrigation events and possible saturations. As a support system for the sensors’ data transmission, we used EM60G (METER) data loggers and ZENTRA cloud (METER) for cloud data storage.
Ground RGB images were taken at several plots (Fig. 1) for ground-based plant spatial abundance assessment. We used a wooden-frame of 1 m by 1 m to delimit the plot’s area and took a picture from above (1.4 m from the ground). This assessment will be used to validate UAV-based plant spatial abundance estimates.
Soil samples were collected at twelve locations (Fig. 1). Physical laboratory testing was performed at the Fayetteville Agriculture Diagnostic Laboratory of the University fo Arkanas to characterize soil properties, such as clay density, porosity, dry bulk density, volumetric water content, and organic matter. Such information was used for statistical analysis and to provide an interpretation of the soil ECa measurements. No plant tissues or seeds were collected in this study.
Yield data from a combine harvester
Crop yield was harvested on October 24th. The combine harvester recorded yield information as cloud-point data in bushels per acre (bu/ac) and reported here in kilogram per hectare (kg/ha), and the relativre geolocation information. The combine harvester is able to collect 11 rows at once and records a point measurement every 2 s. Areas with artefacts were determined based on the swath parameter recorded during the harvesting. The wrong choice of the swath parameter on the combine harvester would produce unreliable yield values. Such areas were not considered in further analysis.
Methodology
We designed a data analysis framework to combine the UAV optical images and soil geophysical measurements for characterizing plant phenological and physiological properties and soil spatial heterogeneity, respectively. We then performed statistical analysis of their co-variability to quantify the influence of soil properties on plant development. This framework is constructed by the following three steps:
- 1.
UAV-based monitoring of plant dynamics: This step presents the data processing pipeline used for computing UAV products to characterize plant phenology and structure, such as plant spatial abundance, plant-specific vigor, and plant height.
- 2.
EMI-based monitoring of soil characterization: This step focuses on quantifying the soil spatial variability. We perform statistical analyses to identify the relationship between soil ECa signal and textural information derived by soil samples collected during the ground data acquisition.
- 3.
Quantifying soil–plant spatiotemporal co-variability: This step focuses on the statistical analysis of plant physiological and structural properties and near-surface properties, as well as identifying key environmental factors.
UAV-based plant monitoring
Figure 2 depicts the data processing pipeline developed to assess the spatial variability of plant characteristics during the growing season. The steps to obtain the plant spatial abundance and plant vigor maps are presented in detail in the following sections.
Pipeline developed for the UAV data processing.
Plant spatial abundance
Plant spatial abundance is estimated by the RGB mosaics. In the first step of the procedure, a vegetation map was obtained by performing a binary classification to discriminate the two classes of vegetation and soil. We used a supervised spectral-spatial image classification approach43 that combines both the spectral and the contextual (i.e., geometrical) information. The contextual information represents the pixel spatial arrangement and the spatial reletionship between adjacent pixels. By modelling such infromation, we can extract structures within the scene. For this step, we use mathematical operators, such as attibute profiles44,45,46 that belong to the image processing sub-field of mathematical morpohlogy. Such operators are 2D filters that act on regions of connected pixels based on the evaluation of an attribute (e.g., scale, defined as number of pixel composing the region). The advantage of using this operator is the ability to preserve the geometrical detail of the structures in the scene. In a filtered image (or feature), regions with similar properties (i.e., scale) are preserved or otherwise merged to their surroundings (i.e., filtered). A multiscale contextual characterization is obtained by applying a filter recursively and relaxing the scale parameter of the filter at each iteration, resulting in a stack of filtered images. An automatic procedure developed in45 is used to perform such a characterization. Such operation serves to extract homogeneous structures present in the scene at different spatial scales, allowing us to better delineate the soybean rows. The original RGB image, together with the filtered images, compose the feature space used as input to a support vector machine (SVM) algorithm with a radial basis function (RBF) kernel. The algorithm is based on the LIBSVM library47 developed for the MATLAB environment, using a one-against-one multiclass strategy. In a second step, we compute the spatial abundance of the detected plants as the percentage of pixels that belong to the vegetation class within a grid unit. In our study we chose a grid unit of 2 by 2 m, which would allow to obtain a good approximation of the local change, while capturing large scale patterns. The quality of the plant spatial abundance map is mainly dependent on the performance of the classification algorithm used to separate the plants from the soil in the background. To validate this product, we computed the classification performance based on ground-truth data and compared the UAV-based estimates of plant abundance to ground-based assessments. The ground assessment of the plant spatial abundance was performed by selecting 31 plots in conjunction of the UAV acquisition occurring in May 28th, 2018. The plots were located mainly along the east and west sides of the field within a distance that ranges between 12 and 25 m from the edge of the field towards the field center; and some plots were collected along the south and north sides within a distance of 50 m from the field edge. These plots were selected to captured representative spatial abundances, ranging from 0% (almost bare) to 60%. Ground images were systematically taken at a height of 1.4 m from the ground. From these images, we computed the GCC and identified a threshold to separate the plant covered area from the soil backgroun. The plant spatial abundance was then computed as the percentage of pixels identified as soybean over the entire area of the plot. All the ground-data were geolocated using the RTK-GPS system previously described.
Plant-specific vigor estimation
Remote-sensing-derived VIs have been extensively used as a proxy to estimate plant vigor and investigate plant physiology and response to possible ecosystem changes. Widely used VIs for plant characterization are in general derived by the near-infrared (NIR) channel (e.g., NDVI, SAVI, etc.), which is sensitive to changes in chlorophyll and leaf area. However, RGB-based indices have been developed and used in several phenology studies, showing their effectiveness in estimating plant characteristics compared to NIR-based VIs48. The GCC was chosen based on authors’ recent studies, showing that this RGB-based VI is less affected by saturation compared to NDVI31. GCC is formally defined as follows:
$${text{GCC}} = frac{{{text{G}}_{{{text{DN}}}} }}{{{text{B}}_{{{text{DN}}}} + {text{G}}_{{{text{DN}}}} + {text{R}}_{{{text{DN}}}} }},$$
(1)
where R, G, and B are the red, green, and blue channels, and DN refers to the pixel intensity values as digital numbers41. The plant-specific vigor was computed considering only those pixels identified as vegetation, the information for which is provided by the binary classification. In order to have products spatially comparable, the plant-specific vigor was transformed into a 2 by 2 m grid by computing the average within each grid unit. This procedure has the advantage to minimize the soil component in estimating plant vigor, avoiding the nonuniqueness of lower productivity versus increased soil coverage. Specifically, it allows a more accurate way of capturing the spatial variability of the plant’s vigor, in particular during the early vegetative stage, when the soil component is overabundant compared to the plant spatial abundance. This is a clear advantage over a satellite-derived VI with a coarser resolution that would provide an underestimation of plant vigor due to soil coverage.
Plant height estimation
Plant height was computed by subtracting the digital elevation model (DEM) of bare soil (i.e., reference DEM) from the digital surface model computed at each acquisition. The reference DEM was computed from acquisition occurred on May 8th, which was planned a week after planting so that it is not affected by surface disturbances by the planter. DEMs were computed using the PhotoScan software (Agisoft) during the mosaic reconstruction phase. We assume that the ground surface elevation stays the same during the growing season.
EMI-based soil characterization
We used the ECa data to evaluate the spatial heterogeneity of soil properties and the temporal variability in soil moisture. We primarily used the EMI measurements associated with the 1 m effective depth, which is the depth range expected to encompass the soybean roots. Temperature corrections were applied to the point-cloud data for correcting the measured soil ECa to the reference temperature of 25 °C. We employed the exponential model presented in Corwin and Lesch49 and developed by Sheets and Hendrickx50, and formally described as follows:
$${text{E}}C_{25} = {text{E}}C_{a} cdot left[ {0.4470 + 1.4034e^{{left( { – T/26.815} right)}} } right],$$
(2)
where ({text{E}}C_left{ {25} right}) represents the corrected soil ECa at the reference temperature of 25 °C, ({text{E}}C_a) represents the measured soil ECa, and ({text{T}}) represents the soil temperature at the time of the acquisition. The model was chosen based on the on the methodological comparison51, which shows that the exponential model provides very similar estimations compared to the widely used ratio model31,52 in the temperature range of 3–43 °C, but with a smaller total residual.
A spatial correction was applied to the data to resolve a spatial shift caused by the distance from the EMI instrument and the actual position of the GPS antenna, which were 4.5 m apart due to the system setup. We performed statistical analysis between EMI data and soil samples collected during the April data acquisition to identify the main control on the spatial variability of soil ECa.
Co-variability among plant and soil signatures
We performed a suite of statistical analyses, based on the high-resolution data layers, on soil and plant characteristics (i.e., plant spatial abundance, plant vigor, plant height, soil ECa, and crop yield) to investigate the impact of soil properties on plant development during the growing season. For the analysis, we considered the yield point-cloud data as geo-reference.
Since the EMI and crop-yield data had a lower vertical resolution (along the y axis) of about 10 m (distance between rows selected for the data collection), we computed the average of each metric within a window of 20 m by 20 m. The exploratory data analysis included run charts and scatter plots to investigate the temporal evolution of plant development during the growing season, and the co-variability between plant characteristics and soil ECa. We evaluated the correlation between variables using Pearson’s correlation, as well as their spatial association. To control for the spatial non-independence, we derived a corrected Perason’s correlation by using the hypothesis testing procedure propsed by Clifford, Richardson, and Hémon (1989)53. The method uses Moran’s I54 to compute the spatial autocorrelation in the spatial data sets and estimate the correct degrees of freedom, which are then used to assess the significance of the correlation. The method uses the Sturges’ formula to identify the number of distance classes55.
Source: Ecology - nature.com
