Monitoring Agronomic Parameters of Winter Wheat Crops with Low-Cost UAV Imagery

Abstract

Monitoring the dynamics in wheat crops requires near-term observations with high spatial resolution due to the complex factors influencing wheat growth variability. We studied the prospects for monitoring the biophysical parameters and nitrogen status in wheat crops with low-cost imagery acquired from unmanned aerial vehicles (UAV) over an 11 ha field. Flight missions were conducted at approximately 50 m in altitude with a commercial copter and camera system—three missions were performed between booting and maturing of the wheat plants and one mission after tillage. Ultra-high resolution orthoimages of 1.2 cm·px⁻¹ and surface models were generated for each mission from the standard red, green and blue (RGB) aerial images. The image variables were extracted from image tone and surface models, e.g., RGB ratios, crop coverage and plant height. During each mission, 20 plots within the wheat canopy with 1 × 1 m² sample support were selected in the field, and the leaf area index, plant height, fresh and dry biomass and nitrogen concentrations were measured. From the generated UAV imagery, we were able to follow the changes in early senescence at the individual plant level in the wheat crops. Changes in the pattern of the wheat canopy varied drastically from one mission to the next, which supported the need for instantaneous observations, as delivered by UAV imagery. The correlations between the biophysical parameters and image variables were highly significant during each mission, and the regression models calculated with the principal components of the image variables yielded R² values between 0.70 and 0.97. In contrast, the models of the nitrogen concentrations yielded low R² values with the best model obtained at flowering (R² = 0.65). The nitrogen nutrition index was calculated with an accuracy of 0.10 to 0.11 NNI for each mission. For all models, information about the surface models and image tone was important. We conclude that low-cost RGB UAV imagery will strongly aid farmers in observing biophysical characteristics, but it is limited for observing the nitrogen status within wheat crops.

1. Introduction

Recently, unmanned aerial vehicles (UAVs) have been introduced into agricultural research to monitoring crops [1]. In contrast to satellite imagery and aircraft-based remote sensing, UAVs can be used frequently during the entire growth period. The main benefits are simple mission planning, instantaneous operation with low man power and imaging below cloud cover [2]. By carrying low-cost commercial camera systems, UAVs provide ultra-high resolution images of the crop canopy due to the low flight altitude. Current developments in photogrammetric algorithms are specifically adapted to the needs of UAV imagery. Thousands of images can be nearly automatically processed via ready-to-use software to produce orthoimages or surface models. Because wheat is the crop grown with the highest acreage [3] there is a strong interest in obtaining spatial and temporal information about the wheat canopy in high resolution, e.g., to adapt nitrogen and pesticide application site-specifically to improve production efficiency [4,5,6].

The periodic monitoring of the canopy biophysical parameters, such as biomass, leaf area index (LAI) and plant height (PHT), is essential to understanding crop development, variations in canopy reflectance and net ecosystem exchange or nitrogen and pesticide demand during the growth season [7,8,9]. These parameters and derivatives are important for precision agriculture, remote sensing, crop modeling, ecosystem modeling and climate modeling. Crop growth models use a wide range of biophysical parameters to estimate future yield as input and validation [10,11]. Biophysical parameters may further deliver vital information about the specific infection situation with fungal diseases to make field-specific decisions on plant protection [12]. The availability of this information at the field scale would enable a new generation of decision support systems that can optimize fungicide application in winter wheat, e.g., the prototype “proPlant expert” precisely recommends maximum application rates for up to three management zones within a field according to the yield expectation [7].

The monitoring of nitrogen (N) is an essential tool for investigating many metabolic and structural processes in maturing wheat plants, such as yield formation and health status. Because N is not immobilized in soils and an abundant reservoir of plant-available N is not present, it is important for optimal crop production to supply N by applying fertilizer throughout plant growth. However, excessive use of fertilizer eventually leads to unwanted N-leaching into groundwater or N-contamination of surface run-off water, contributing to environmental pollution, which is to be avoided [13]. The N status cannot be estimated from the leaf nitrogen content alone but biophysical parameters of the wheat canopy should also be taken into account. For example, the nitrogen nutrition index (NNI), i.e., the ratio between the actual N concentration (N_t) observed in the plants and the critical N concentration related to dry biomass (N_ct), has been used to reliably characterize the N status of wheat crops during the vegetative period [14]. NNI values >1 indicate excessive N supply (over-fertilization), whereas suboptimal N supply is indicated by values <1. Especially during flowering and the grain-filling growth stage, NNI is a good indicator of N deficiency affecting grain yield and grain protein [14,15].

Due to the complex interrelationship between many environmental factors, such as soil heterogeneity, cultivation and land surface, the parameters described above show high spatial and temporal variability so that a high measurement density would be needed to reflect their spatial patterns within the field [16]. It has been shown that site-specific management strategies in the context of precision farming increase management efficiency [4]. However, accurate measurements of these parameters are time consuming or destructive. For example, the calculation of the NNI involves the cutting of fresh biomass, subsequent drying and the determination of N_t using the Dumas or Kjeldahl method in the laboratory [17].

To become more efficient, indirect methods using sensors for estimating those parameters have been proposed and implemented. At the plot scale, sensor principles are available that enable on-spot measurements without destructing the canopy, mostly with direct contact to wheat plants. For example, leaf chlorophyll meters relate the light transmittance or fluorescence on leaf parts to chlorophyll content and, by using certain indices, to leaf N_t [18]. The LAI can be modeled by devices measuring sunlight interception in the wheat canopy using radiative transfer models [19]. Most of these principles involve time-demanding measurements or sophisticated measurement protocols that can only be performed manually in a stop-and-go mode. Their use for online or high-throughput measurements was discussed only very recently in research, e.g., mobile LAI determination by sunlight interception [20].

For the real-time determination of biophysical parameters and N status, many active sensing approaches have been studied, such as LiDAR, optical sensors and ultra-sound sensors [21,22,23,24], in addition to more refined solutions, e.g., the pendulum operating measurement principle Crop-Meter [25,26]. These sensors are mainly installed within or up to two meters above the canopy and are bound to a ground-driven vehicle such as a tractor. They have the advantage that they deliver information about the wheat canopy immediately and in high resolution. Optical sensors mostly use specific spectral vegetation indices (SVI) in the visual and infrared part of the spectrum. The GreenSeeker^® (Trimble, Sunnyvale, CA, USA) detect the reflection between 656 nm (VIS) and 774 nm (NIR), whereas the CropCircle^® (Holland Scientific Inc., Lincoln, NE, USA) allows for more freedom in the choice of wavelengths by using different filters [23]. To map larger areas, spectral vegetation indices (SVI) from satellite-based and aircraft-based remote sensing imagery have been related to those parameters [27,28]. In contrast to proximal sensing, these data have lower spatial and temporal resolution. UAV imagery might be able to close the gap between the plot scale, covered by manual sensors and proximal sensors, and the regional scale, covered by traditional remote sensing, because of their high spatial resolution and almost instantaneous availability for practitioners and experts in agriculture. Most studies have investigated the relationship between biophysical parameters and UAV imagery, e.g., biomass [29,30,31,32,33], LAI [34,35,36,37,38], plant height [30] and grain yield [39], whereas fewer studies have shown the relationship between nitrogen and UAV imagery [31,35,40,41] in wheat crops, so far. We summarized existing research studies relating UAV imagery with some agronomic parameters of wheat crops in Table S1. Only Pölönen, et al. [31] investigated dry biomass and total nitrogen content in combination with UAV hyperspectral imagery, which allowed them to calculate the under or over nourishment of the wheat plants during crop growth using the NNI. To the best of our knowledge, low cost UAV imagery has not been tried for estimating NNI in cereal crops. Moreover, many studies examined only fields or parts of fields with sizes smaller than 3 ha.

The objective of our work was therefore to study the prospects of low cost UAV imagery for monitoring biophysical plant parameters, N-content and NNI of wheat crops from a more practical viewpoint closer to agricultural routines. Three UAV missions were carried out between booting and maturity over an 11 ha field under soil-induced water deficit conditions and we acquired a large set of aerial images with a consumer level camera. Based on that, ultra-high resoluted orthoimages and surface models were computed through photogrammetric processing. We demonstrated the applicability of UAV imagery for recognizing the spatio-temporal patterns of the wheat canopy development, as observed in this field. Furthermore, we investigated the relationship of biophysical parameters, i.e., plant height, LAI and biomass, as well as, nitrogen status, with the image variables that were derived from the UAV imagery at specific dates. Regarding this relationship, we determined its strength to model the spatial variability of those parameters with linear regression in order to fulfil precision agricultural needs. This is important in order to establish a monitoring system for crop variability that allows fast production of maps for end users such as farmers or agricultural services adopting precision agriculture. This may improve typical applications of precision agriculture such as variable rate fertilization or precision plant protection.

2. Materials and Methods

2.1. Test Site

The study was conducted within a field in Eastern Germany during the spring season in 2015 (51°49′N, 12°42′E). The soil development in the field was influenced by recent flood plain deposits of the Elbe River. The main soil type is a fluvic cambisol and the soil texture varies between sand, loamy sand and sandy loam. The crop grown was winter wheat (Triticum aestivum L. var. ‘Linus’). The seed row distance was 0.12 m, and the average crop density was between 440 and 480 ears·m⁻².

2.2. UAV System and Flight Missions

We used a hexacopter (P-Y6, Hexapilots, Dresden, Germany) to carry a commercial camera system to acquire the aerial images for this study (Figure 1). The hexacopter consisted of six propellers aligned in a three arm mounting with two propellers attached to each arm using the push/pull principle for failure safety. The navigation control system was a DJI Wookong M (DJI Innovations, Shenzhen, China) with an integrated GNSS receiver, which enabled user-defined waypoint tracking. Power supply was ensured with lithium polymer batteries (10,000 mAh/5 s). With this setup, including the camera, approximately 15 min flight duration was achieved.

For image acquisition, we used a Sony Nex 7 point and shoot camera with the following specifications: 24 megapixels, 23.5 × 15.4 mm sensor and E16 mm F2.8 fixed lens (Sony Corporation, Tokyo, Japan). The camera was mounted onto a gimbal underneath the copter. The DJI software triggered image capturing, and the position of the copter was recorded by a GPS device.

The flight missions were planned to cover approximately 13 ha of ground area. The orthoimages and surface models were later clipped to 11 ha to correspond to the study area. Parallel tracks were flown with a between-track distance of 11 m at an altitude of 50 m which is a side-lap of 60%. The frequency of image capture was set to correspond to 60% image overlap. All images were taken from near nadir position with a ground resolution of approximately 0.012 m·px⁻¹.

2.3. Ground Truthing

Three flight missions (M1-3) were conducted over the wheat canopy during crop growth from booting of the wheat plants until mature growth stages. An additional flight mission (M4) was performed after tillage to yield a UAV image of the ground surface. Mission dates, weather conditions and objectives of the flight missions are enlisted in

Table 1. Flight mission objectives and environmental conditions during overflight.

Mission Objective	Processed Item	ID	Date	Growth Stage	Cloud Cover [okta]	Wind Speed [ms−1]
Crop canopy	Orthoimage/Surface model	M1	05-18-2015	BBCH * 41–47 (Booting)	8	2–3
Crop canopy	Orthoimage/Surface model	M2	06-04-2015	BBCH 61–71 (Flowering)	5	1–2
Crop canopy	Orthoimage/Surface model	M3	06-16-2015	BBCH 73–83 (Grain filling)	7	4
Ground model	Surface model	M4	07-31-2015	After tillage	8	4

Note: * Biologische Bundesanstalt, Bundessortenamt und Chemische Industrie (BBCH).

.

For each mission, 20 locations were chosen with respect to the wheat variability observed at that date. These locations were selected along the tractor lanes by eye. Locations were at least 2 m distance from the tractor lane in order to be non-influenced by the clearing. At each plot, the ground area spanned by 1 × 1 m² was taken as the sample size. For georeferencing the UAV images, 20 panels were laid out along a regular grid over the field within the tractor lanes with good visibility from above before the flight missions. Panels and plots were located using the differential GPS HiPer Pro system (Topcon Positioning Systems, Inc., Livermore, CA, USA) having a relative horizontal and vertical accuracy of 3 mm and 5 mm.

At each plot, the plant height (PHT), LAI, fresh biomass (FBM), dry biomass (DBM) and N_t were determined as agronomic reference measurements. The descriptive statistics and the abbreviations used in this paper are summarized in

Table 2. Descriptive statistics of the crop parameters measured at the sample plots for M1, M2 and M3.

Variable	Abbreviation	Unit	Min	Mean	Max	Standard Deviation	Median
Mission 1
Fresh biomass	FBM	kg·m−2	1.64	3.07	4.72	0.959	2.66
Dry biomass	DBM	kg·m−2	0.63	0.87	1.07	0.128	0.88
Leaf area index	LAI		1.77	3.08	4.81	0.956	2.59
Nitrogen	Nt	%	1.44	1.95	2.75	0.298	1.94
Plant height	PHT	m	0.44	0.58	0.7	0.076	0.58
Mission 2
Fresh biomass	FBM	kg·m−2	1.56	2.88	5.02	1.099	2.43
Dry biomass	DBM	kg·m−2	0.66	1.02	1.46	0.051	0.96
Leaf area index	LAI		2.15	3.66	6.29	1.140	3.52
Nitrogen	Nt	%	1.47	1.77	2.06	0.205	1.78
Plant height	PHT	m	0.45	0.61	0.82	0.111	0.58
Mission 3
Fresh biomass	FBM	kg·m−2	0.90	3.04	5.64	1.340	2.99
Dry biomass	DBM	kg·m−2	0.47	1.23	1.76	0.375	1.27
Leaf area index	LAI		1.07	2.77	5.42	1.069	2.57
Nitrogen	Nt	%	1.02	1.36	1.70	0.192	1.36
Plant height	PHT	m	0.30	0.62	0.80	0.154	0.67

. Wheat plants were classified according to the BBCH growth stages code of Lancashire et al. [42]. Measurements of PHT were integrated by the measured heights of the wheat plants using a folding yardstick at 10 random locations within the plot. LAI was measured using the SunScan Canopy Analysis System type SS1 (Delta-T Devices, Cambridge, UK). The LAI probe was positioned at 45° to the direction of the seed rows under the canopy, and the reference sensor was in immediate proximity above the canopy without disturbing the incoming sunlight. The LAI measurement was taken as an average of 10 individual measurements with the probe repositioned each time within the plot. All wheat plants within the plot area were then cut directly above the ground using a short reaping hook. FBM was determined immediately after cutting. DBM was measured by drying biomass samples at 60 °C for 48 h in a compartment drier. The dry plant material was milled and analyzed for N_t following the dry combustion method of Dumas [43] using an elemental analyzer (Vario MAX CN Elemental Analyser, Elementar, Hanau, Germany).

2.4. Image Pre-Processing

All aerial images were recorded in RAW data format and lossless converted to the tagged image file format (TIFF) in standard RGB color space. The full camera calibration matrix, including non-linear distortion coefficients, was estimated for the camera sensor using the Agisoft Lens software (v. 0.4.2, Agisoft LLC, St. Petersburg, Russia) by repetitively taking photographs of a checkerboard pattern. The Agisoft Lens software uses a pinhole camera model for lens calibration, and the distortion correction is modeled by Brown’s distortion model. We used a number of radiometric pre-processing steps to improve the final UAV image mosaicking and the surface model generation results, as shown below. These algorithms were programmed in Matlab 2015b. To reduce colorimetric alterations, the following corrections were conducted only on the value channel of the HSV transformed aerial images. The corrected HSV images were then transformed back to RGB space.

To diminish the effect of brightness reduction from the image center to the borders due to the sensor optics, all images were empirically corrected using the vignette correction algorithm proposed by Zheng et al. [44]. This processing involved calculating a mean image from a selection of homogeneous aerial images (n = 80). In our case, images were selected that depicted only the wheat canopy without prominent features such as trees. The further pre-processing steps were performed on the vignette-corrected images.

For better photogrammetric processing, contrast-enhanced images were produced from each aerial image. This was achieved by contrast limited adaptive histogram equalization (CLAHE). Each image was divided into eight tiles. Within each tile, the contrast was enhanced by histogram equalization following the Gaussian distribution. To avoid artificial boundaries, the neighboring tiles were then joined using bilinear interpolation. The images corrected by CLAHE were used during the mosaicking process and surface model generation.

Different incidence and viewing angles may result in unwanted radiometric variations in the aerial images related to the bidirectional reflectance distribution function (BRDF). This effect was reduced following the empirical BRDF correction algorithm suggested by Lelong et al. [35]. In short, each image was block-wise averaged to a smaller representation of the image by bilinear interpolation, then smoothed by Gaussian filtering and finally over-sampled onto the original image with bi-cubic interpolation. Then, the inverted, smoothed values were subtracted from the original values. Specifically, all images were subdivided into 400 × 400 pixel blocks, and the Gaussian filter size used was a 3 × 3 pixel window. The BRDF set of images was used for the final texturing of the orthoimage.

2.5. Photogrammetric Processing

The vignette-corrected and CLAHE pre-processed images were used in the semi-automated processing flow of the photogrammetry software Agisoft PhotoScan (v. 1.2.4, Agisoft LLC, St. Petersburg, Russia) to generate orthoimages and surface models for M1-4. The software implements Structure from Motion to estimate a 3D point cloud from the overlapping aerial images [45]. Using feature detection and description algorithms, key points between overlapping images that bear geometrical similarities in their immediate surroundings were selected. In the first step, a sparse 3D point cloud was generated to align the images and estimate the exact camera positions. We used the options ‘high’ and ‘generic’. In the second step, a dense, multi-view stereo reconstruction on the aligned image set operating on the pixel level was applied that generated a dense 3D point cloud. For this, we used the options ‘medium’ and ‘mild’ depth filtering. We used this setting because higher quality takes exponentially more time and demands more computational resources. Prior tests with the software showed “medium” as a good compromise with adequate processing time (4.6 cm/pix resolution error). The mesh was generated using the option ‘height field’ to produce the orthoimages and the surface models. To texture the orthoimages, the CLAHE pre-processed images were exchanged with the BRDF-corrected images. The texturing was performed with the software-implemented color correction. The ground resolution of the orthoimages was 0.025 m·px⁻¹.

2.6. Extraction of Image Variables

We calculated a set of image variables from the orthoimages and the surface models at each measurement plot to relate them to the agronomic reference measurements. The image variables were computed from averaging the pixels within the plot areas. Furthermore, crop pixels were identified to calculate crop coverage (CVR). This was achieved by first converting the respective orthoimages into the LAB color space because we found the best discrimination between soil and vegetation in the a-vector corresponding to the red-green variation of the images. Secondly, within the whole orthoimage, we randomly selected regions of interest (ROI) containing only soil-related pixels. These ROIs were combined into a simple vector and the value given at the 90% percentile was used as a threshold value to differentiate between soil and vegetation pixels in the measurement plots. CVR was computed as the percentage of all vegetation pixels identified within a plot. All image variables are described in

Table 3. Image variables description and calculation basis.

Image Variable	Description/Computation
CVR	Percentage of crop pixels within plot area
EXG 1	Related to the green channel (EXG = 2 green channel − red channel − blue channel)
RED 1	Red channel
BG 1	Ratio of the blue and green channel
RG 1	Ratio of the red and green channel
RB 1	Ratio of the red and blue channel
PHTUAV 2	Plant height from UAV images. Computed by subtracting the M4 surface model (ground surface model) from the M1-3 surface models.

Note: ¹ Before computation, the red, green and blue channels were normalized by the sum of all channels.

.

2.7. Statistical Analysis

The relationships between the image variables and reference measurements were investigated and tested for significant correlation using a Pearson correlation matrix for M1-3. Principal component analysis (PCA) was conducted to explore the structural variation and dimensionality within all image variables for M1-3. PCA was calculated using the correlation matrix to eliminate the influence of different standard deviations among the image variables. Loading plots were used to summarize the interrelationship between the image variables. Linear regression analysis was conducted using the agronomic reference measurements as dependent variables. As independent variables, the scores from the principal components (PC) that contributed sufficiently to the overall variance of the image variables were used to prevent multi-collinearity in the regression models. The best models were chosen by backward selection. Variables with p > 0.1 were deleted from the set of independent variables [46]. For model comparison, we reported the adjusted R² values according to Wherry [47].

The NNI was calculated using the approach of Lemaire and Salette [48]:

$$N_{ct}=aDB{M}^{-b},$$

(1)

$$NNI=\frac{N_t}{N_{ct}}$$

(2)

where N_ct is the critical N concentration related to a specific dry biomass (DBM). Justes et al. [15] specified the dilution curve for N_ct statistically on wheat crops by proposing values for a = 5.35 and b = 0.442, which were used for the NNI determination in this study.

All models were validated using an independent model and validation data set by splitting the field equally in two halves orthogonal to the tractor lane direction. As validation criteria, the R² of validation (R²val) expressed as the squared Pearson correlation coefficient, the root mean squared error (RMSE) and the mean error (ME) were calculated:

$$R^{2}val=\frac{cov\left(x,y\right)}{var\left(x\right)var\left(y\right)},$$

(3)

$$RMSE=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{x}_i{y}_i}{n}},$$

(4)

$$ME=\overline{x}-\overline{y}$$

(5)

where x and y denotes the reference and predicted values and $\overline{x}$ and $\overline{y}$ the means of the reference and predicted values at the validation locations, n the number of validation points and cov and var the covariance and variance, respectively.

3. Results and Discussion

3.1. Qualitative Assessment of the Acquired UAV Imagery

Three UAV missions were flown during the growing season on 5/18/15 (M1), 6/4/15 (M2) and 6/16/15 (M3). The phenological stages of the wheat plants ranged from booting with flag leaf sheath extending (M1) to flowering (M2) to maturity (M3) (Figure 2). At the specific dates, the growth stages varied to some extent spatially across the fields due to soil heterogeneity. Generally, at locations with coarser soil texture, where water availability was lower, the growth stages tended to be more mature, whereas plant biomass was underdeveloped. This differentiation gradually increased at later dates. Plant diseases were not detected throughout all mission dates.

The distribution of the RGB colors in the orthoimages was structured according to the spatial variation of the plant coverage, growth stage, leaf vitality and degree of senescence (Figure 3). Darker greenish colors indicate areas with denser crop canopy and higher plant vitality, whereas lighter yellow-brownish colors are associated with sparser crop canopy and lower vitality. Soil-induced water deficit stress reduced leaf area development and plant height, causing spatial variation in the crop canopy throughout the field. It also invoked an earlier senescence of the leaves [49]. For M1, the lighter-colored areas corresponded to average plant water content: 64%, LAI: 2.27 and dry biomass: 0.80 kg·m⁻², whereas within the greenish areas, the same parameters were higher, with plant water content: 74%, LAI: 3.61 and dry biomass: 0.92 kg·m⁻². The effect of the soil-induced water deficit stress increased from M1 to M3, leading to stronger structuring in the orthoimages with an advancing image tone to more yellow and brown colors. This structuring and pattern indicates near-surface sedimentological changes. They were established by the remnants of the former fluvial geomorphology in the Pleistocene shaped by the floodplain of the river Elbe. At that time, differences in the river flow rates caused a braided river system with lots of small river channels broken by temporary islands. This resulted in highly local differences in soil texture, which is typical for the fields in this area and invoked the patterning of the wheat canopy. As observed from the zoomed representations depicted in Figure 4a–c, the wheat canopy completely changed within only a few days from M2 to M3. That change cannot be observed from M1 and can only be vaguely outlined at M2. This sudden change emphasizes the need of high temporal availability of aerial image data for cultivated land—an advantage of UAV imagery over remote sensing imagery by satellite and airplanes due to its independence of cloud conditions, its easy mission planning and low costs. Moreover, the high spatial resolution of the UAV images allows the differentiation of these changes down to the individual plant level (Figure 4d). Using these images in a geographical information system, farmers would be able to locate and monitor changes in their crops with high accuracy over the entire field. Depending on the aim of the qualitative assessment of the farmer this zoom level may become useful, e.g., plant diseases (yellow rusts) spread from small nests.

The same strong patterning as in the UAV imagery occurred in the surface models (M1-3) of the wheat canopy. As observed in Figure S1, the patterning of M1 was not the same as the ground surface model (M4). Thus, changes in the surface model of M1 may show differences in plant growth.

3.2. Relationship between Agronomic Parameters and Image Variables

The biophysical parameters FBM, DBM, LAI and PHT were highly correlated with each other at M1-3 (

Table 4. Correlation matrix (Pearson) of the crop parameters (FBM—fresh biomass, DBM—dry biomass, LAI—leaf area index, PHT—plant height, N_t—total nitrogen content) and image variables (CVR—coverage, EXG—green SVI, RED—red channel, BG—blue green ratio, RG—red green ratio, RB—red blue ratio, PHT_UAV—plant height calculated from UAV). The grey color indicates the matrix area with correlations of reference measurements vs. image variables.

Mission	Variable	FBM	DBM	LAI	PHT	Nt	CVR	EXG	RED	BG	RG	RB
M1	DBM	0.87
	LAI	0.98	0.82
	PHT	0.94	0.84	0.90
	Nt	0.44 ns	0.23 ns	0.45	0.27 ns
	CVR	0.71	0.63	0.69	0.77	−0.01 ns
	EXG	0.69	0.52	0.73	0.73	−0.05 ns	0.83
	RED	0.54	0.56	0.51	0.41 ns	0.53	−0.02 ns	−0.16ns
	BG	−0.95	−0.82	−0.97	−0.93	−0.26 ns	−0.78	−0.82	−0.42 ns
	RG	0.02 ns	0.11 ns	−0.01 ns	−0.09 ns	0.43 ns	−0.47	−0.68	0.82	0.14 ns
	RB	0.88	0.80	0.87	0.78	0.52	0.43 ns	0.35 ns	0.85	−0.82	0.44 ns
	PHTUAV	0.83	0.68	0.83	0.88	0.17 ns	0.76	0.74	0.26 ns	−0.85	−0.18 ns	0.65
M2	DBM	0.97
	LAI	0.94	0.93
	PHT	0.94	0.90	0.90
	Nt	−0.06 ns	−0.04 ns	−0.16 ns	−0.19 ns
	CVR	0.81	0.78	0.82	0.86	−0.53
	EXG	0.74	0.65	0.75	0.81	−0.48	0.89
	RED	0.78	0.82	0.70	0.71	0.29 ns	0.47	0.26 ns
	BG	−0.92	−0.86	−0.89	−0.95	0.30 ns	−0.92	−0.93	−0.59
	RG	0.08 ns	0.20 ns	0.01 ns	−0.04 ns	0.68	−0.31 ns	−0.56	0.63	0.23 ns
	RB	0.97	0.95	0.91	0.93	0.05 ns	0.75	0.70	0.86	−0.90	0.19 ns
	PHTUAV	0.86	0.81	0.84	0.92	−0.38 ns	0.94	0.88	0.55	−0.94	−0.24 ns	0.82
M3	DBM	0.98
	LAI	0.96	0.92
	PHT	0.92	0.94	0.85
	Nt	−0.29 ns	−0.32 ns	−0.19 ns	−0.49
	CVR	0.88	0.92	0.77	0.96	−0.50
	EXG	0.94	0.90	0.94	0.83	−0.35 ns	0.80
	RED	0.87	0.80	0.85	0.74	−0.03 ns	0.71	0.78
	BG	−0.97	−0.92	−0.95	−0.86	0.28 ns	−0.84	−0.98	−0.88
	RG	0.41 ns	0.36 ns	0.38 ns	0.32 ns	0.31 ns	0.31 ns	0.19 ns	0.74	−0.39 ns
	RB	0.94	0.87	0.94	0.80	−0.14 ns	0.75	0.90	0.96	−0.96	0.58
	PHTUAV	0.92	0.95	0.85	0.98	−0.56	0.94	0.84	0.68	−0.85	0.23 ns	0.77

^ns not significant at the 95% probability level.

). Specifically, the correlation between LAI and FBM was high because plant water is mainly stored in plant leaves, and plant water varied significantly in this field.

The correlations between the biophysical parameters and the image variables were significant in most cases. Overall, the average correlation strength between all image variables and the biophysical parameters increased from M1 to M3 (r = 0.64 < 0.73 < 0.82). This increase can be explained by the higher variation and by the better differentiation of the image tone that was observed at the later missions in the field. The image variables BG and RB were the most strongly correlated. Moreover, PHT_UAV, calculated from the surface models, had a valuable contribution to explaining the biophysical parameters of the wheat canopy. From M1 to M3, the correlation of PHT_UAV with FBM, DBM, LAI and PHT was highly significant (r > 0.80). The only exception occurred at M1, in which PHT_UAV was somewhat more weakly correlated with DBM (r = 0.68). The weaker correlation can be explained by dry matter being mainly stored in the wheat leaves and not in the wheat stems at booting stage. The dry matter in the leaves does not significantly contribute to the actual plant height [50]. This result can also be anticipated from the lower correlation of PHT with DBM at M1 compared to M2 and M3.

In contrast, N_t was only weakly or not significantly correlated with the biophysical parameters throughout all missions. We expected a significant correlation between RED and N_t because a cross-correlation between the chlorophyll content in plant leaves and the nitrogen content exists. Because chlorophyll absorbs light in the red wavelengths, it may cause variation in the RED image variable or red-related SVIs, which could be related to N_t [51,52]. The correlation of N_t with the image variables was generally higher than correlation with the biophysical parameters but it was still relatively weak. On M1 and M2, significant correlations between RED and N_t were found. However, at M3, no significant correlation existed, which can be explained by the large appearance of senescent plant leaves with no or little chlorophyll content left. The linear relationship between the chlorophyll content and N_t breaks down when the vegetation becomes senescent [28]. As the wheat reaches maturity, wheat leaves lose chlorophyll and N is used for grain development [14]. The best relationship between the chlorophyll content and N_t was found near the flowering growth state for wheat [14,18]. RED and RG had the highest correlation with N_t at M2, whereas, at M3, the biophysically related image variables PHT_UAV and CVR were negatively correlated with N_t and RED was not correlated.

In Figure 5, the scree plots and the loading plots for the first few components were computed by PCA using the correlation matrix of the image variables for each of the missions as extracted at the measurement plots. The scree plots show the component variances against the components and the loading plots show the degree of the relationship to a specific PC. Loading vectors that point in the same or opposite direction indicate positively or negatively correlated variables whereas orthogonal loading vectors highlight dissimilar variables.

Most of the variance explained by the image variables was projected onto the first three principal components of the PCA as observed from the eigenvalues depicted in the scree plot (Figure 5). The loading plots show a specific multivariate grouping of the image variables. Two main groups can be identified: the image variables influenced by the red signal (RED, RG, RB) and the image variables related to the blue signal and biophysical parameters (BG, EXG, CVR, PHT_UAV). On M1 and M2, the first group was mainly related to PC 2 whereas the biophysical group was related to PC 1 and PC 3. At M3, the groups became similar.

3.3. Modeling the Biophysical Wheat Parameters

We performed linear regression on PC 1-3 as independent variables. The relevant PCs were chosen by backward selection with p < 0.1 (

Table 5. Linear regression results with validation for the biophysical crop parameters. The principal components of the image variables (PC) were used as independent variables.

Variable	Mission	R2	Significance	PC	R2val (n = 10)	RMSE (n = 10)	ME (n = 10)
FBM	M1	0.92	***	1***, 2***	0.93	0.46	−0.38
FBM	M2	0.93	***	1***, 2***, 3	0.87	0.53	−0.33
FBM	M3	0.97	***	1***, 2**	0.99	0.48	−0.40
DBM	M1	0.70	***	1***, 2	0.73	0.10	−0.06
DBM	M2	0.89	***	1***, 2***	0.83	0.14	−0.08
DBM	M3	0.94	***	1***, 2**, 3**	0.94	0.21	−0.15
LAI	M1	0.94	***	1***, 2***, 3*	0.96	0.46	−0.41
LAI	M2	0.83	***	1***, 2	0.90	0.48	−0.15
LAI	M3	0.90	***	1***, 3*	0.93	0.48	−0.36
PHT	M1	0.87	***	1***	0.87	0.06	−0.05
PHT	M2	0.93	***	1***, 2	0.90	0.06	−0.05
PHT	M3	0.96	***	1***, 2***, 3***	0.88	0.15	−0.11

Significant at the * 0.05; ** 0.01; or *** 0.001 probability level.

). All models were significant at p < 0.001. The best overall model fits were obtained at the later missions (M2 and M3) for FBM and PHT (R² > 0.92), and the lowest model fit was achieved by DBM at M1 (R² = 0.70). Thus, model quality increased towards the later dates. This gradient was most evident for DBM, due to the reduction of plant water within the wheat canopy caused by senescence of the plant leaves. The image tone increasingly represented the variation of dry biomass in the orthoimages of the later missions. Additionally, the contribution of PC 3 is more important at M2 and M3 because the plant height has a larger influence on DBM at flowering and maturity than at booting. Overall, the contribution of the PCs varied between the models, and in almost all cases, more than one PC was significant, which emphasizes the need to use different image variables to estimate the biophysical parameters. For all models, PC 1 was a highly significant predictor. Since CVR, PHT_UAV and BG load highest on PC 1, it shows that these variables are relevant for determining the biophysical parameters from UAV RGB imagery. However, the red-related image tone variables, i.e., RED, RG and RB, which loaded strongly on PC 2, had an influence on most models. Interestingly, the influence of PC 2 was relatively high for determining PHT, which means that the image tone may also support the estimation of plant height derived from the UAV crop surface models. Moreover, the inclusion of plant height information helps to linearize the other biophysical parameters in wheat crops when observed over several growth stages [53]. Therefore, it is advisable to use both the image tone and the surface models to estimate the biophysical parameters due to their strong interrelationship within the wheat canopy.

3.4. Modeling the Nitrogen Status

The models for N_t were weaker than the models for the biophysical parameters (

Table 6. Linear regression results with validation for the nitrogen content. The principal components of the image variables (PC) were used as independent variables.

Variable	Mission	R2	Significance	PC	R2val (n = 10)	RMSE (n = 10)	ME (n = 10)
Nt	M1	0.22	*	2*	0.43	0.20	−0.13
Nt	M2	0.65	***	1*, 2***, 3**	0.64	0.14	−0.08
Nt	M3	0.40	**	2**, 3*	0.24	0.17	0.05

Significant at the * 0.05; ** 0.01; or *** 0.001 probability level.

). The importance of PC 2 specifically at M1 and M2 was related to the influence of the red-related image tone variables because of their dependence on the chlorophyll content. A highly significant model was only obtained for M2 (p < 0.001) during flowering due to the well-established physical relationship between the chlorophyll content and N_t that can be observed during this growth stage. However, PC 3 also had a significant influence at M2 and M3 because of the increasing influence of PHT_UAV and the declining influence of the red image tone variables for determining N_t as was shown earlier in the correlation matrix.

To calculate the NNI, we combined the DBM and N_t estimations using Equations (1) and (2). The validation results are given in

Table 7. Nitrogen nutrition index (NNI) validation results based on the linear regression results of N_t and DBM calculated by Equations (1) and (2).

Variable	Mission	R²val (n = 10)	RMSE (n = 10)	ME (n = 10)
NNI	M1	0.73	0.11	−0.09
NNI	M2	0.58	0.11	−0.08
NNI	M3	0.37	0.10	−0.03

. For each mission, we obtained the same error of 0.11 NNI. The NNI reference values indicated a slight undernourishment of the wheat crops, with most NNI values smaller than 1.0, and the variability of the NNI was relatively small for the heterogeneous field (Figure S3). The scatter plots showed a linear relationship between validations and predictions but were influenced by some extremely erroneous predictions, specifically at M3.

4. Discussion

Compared to studies investigating proximal sensing to estimate biophysical crop parameters, the UAV image variables were competitive. Several studies tested LiDAR systems from a tractor as an online proximal sensing system for adapting application on the go. They related either the vegetation volume or the reflection height (distance between crop surface and laser scanning unit) derived by triangulation or time-of-flight of the reflected laser beam with biomass and LAI measurements. The studies reported R² values for FBM between 0.77 and 0.99, for DBM between 0.72 and 0.99 and for LAI between 0.70 and 0.96 [22,54,55]. Sensing LAI was further investigated using a combined radiometer and ultra-sound sensing system yielding an R² value of 0.84 [24] and a mobile LAI system based on canopy light transmittance yielding R² values ranging from 0.73 to 0.86 [20] between the sensor and LAI reference measurements.

These sensing approaches were not developed to map entire fields but to calculate application rates from sensor measurements online. Similarly, UAV systems might be used this way by having a cable-tethered UAV system flying in front of the application unit. This would enable greater influence on the field of view because the UAV system’s flight altitude would not be fixed, in contrast to the online sensors presented above [56]. Nevertheless, offline approaches, which are more related to UAV, can represent an entire field as a parameter map. Specifically satellite remote sensing products have been related to various biophysical parameters due to their large spatial coverage and their possibility for regional and global upscaling. Studies that related common SVIs used in satellite-based remote sensing such as normalized difference vegetation index (NDVI) or enhanced vegetation index (EVI) with biomass and LAI reported relationships in a wide range from non-significant to highly significant. These results depended on the specific SVI, the spatial and spectral resolution of the satellite system and the growth stage [27,28,57,58]. The newer commercial satellites make it now possible to depict the earth’s surface in sub-meter resolution. WorldView-3 (DigitalGlobe Inc., Longmont, CO, USA.), for example, delivers imagery with 0.30 m (panchromatic) and 1.1 m (RGB) ground resolution [58] with a revisiting time of less than five days. These properties are quite competitive with airborne or even UAV platforms. However, cloud coverage may limit the use of satellites for monitoring during crop growth. In addition, high scene prices and/or special order rules such as selling only large areas (10 × 10 km) may hinder the adoption of these data for precision agriculture.

UAV systems can include information about the canopy volume of the wheat crops estimated from point clouds calculated by the overlapped UAV imagery in addition to the spectral information contained in the image tone which might be well suited to estimate biophysical crop parameters. Bendig et al. [29] showed highly linear relationships between crop height measurements and plant heights calculated from UAV surface models with an R² of 0.92 over multiple growth stages in a controlled test site with different cultivars and N-treatments. Grenzdörffer and Zacharias [59] reported an R² between 0.60 and 0.72 within different cultivars. For comparison, when we pooled all the plant height measurements from the missions into one data set, we obtained a relationship with an R² of 0.85 (Figure S2). Hunt et al. [60] related the NIR, green and blue signal from a commercial camera system recorded from an airborne platform with DBM and LAI measurements and reported significant relationships with an R² of 0.65 and 0.85, respectively. LeLong et al. [35] modeled LAI with NDVI on the basis of UAV imagery. They reported an RMSE between 0.5 and 0.6 LAI shortly before and after flowering of durum and bread wheat.

In contrast to the biophysical UAV models, estimating the N status from UAV imagery was only moderately competitive with studies investigating proximal or hyperspectral sensing. Erdle et al. [23] investigated for its accuracy of N_t and NNI calculation a number of different proximal sensing systems used for online application such as the GreenSeeker^® and the CropCircle. They reported R² values between 0.41 and 0.83 for N_t or 0.52 and 0.91 for NNI shortly before flowering. Pölönen et al. [31] investigated a hyperspectral UAV system and found an R² value of 0.72 for N_t models. The study of Chen et al. [61] investigated the relationship between different SVIs and nitrogen status calculated from hyperspectral field measurements. They reported an error between 0.13 and 0.37% for N_t and an error between 0.13 and 0.17 for NNI shortly before booting. The R² ranged from 0.28 to 0.92, with an error of 0.82 to 0.90 for NNI. However, they observed a larger range of N_t and NNI values compared to the range in this study. They concluded that a mechanistic model that combines SVIs with a clear theoretical basis, i.e., SVIs related with N_t and DBM, should be more conclusive for the estimation of the NNI rather than the arbitrary use of SVIs, e.g., the NDVI. In the same way, UAV imagery can be considered a large improvement over airborne and satellite imagery because the combination of crop surface models with image tone variables provides a more direct approach to understand and model the NNI spatial distribution within fields. In addition, the estimation of N_t is also easier because of UAV imagery. Because the senescence of plant leaves is increasingly propagating within the wheat canopy during the later growth stages, biophysical parameters become increasingly important to describe N_t in the wheat canopy, and the use of PHT_UAV derived from crop surface models becomes more sensible for estimating N_t.

Online sensing delivers instant information about the crop variability. This is not possible with UAVs today. However, online sensing has the strong limitation that it is ground and vehicle based. The field of view of the measurement is quite narrow and linear along the driving direction. In contrast, UAV images deliver spatial information about crop variation over the entire field without interpolation as is with high spatial and temporal resolution. The offline approach of UAVs would also allow for more insight rather than the “black box” online sensing approach.

According to Colomina and Molina [62], UAV missions should use an 80%–90% image overlap to compensate for the instability of the platform. However, we did not experience any problems with image alignment and dense cloud computing in Agisoft using the 60% overlap setting. Therefore, we decided to use this more economical approach. Using the 80% setting would have increased mission time about 50% and doubled up lithium polymer batteries (10,000 mAh/5 s). In addition, memory storage for images and processing time for photogrammetry would have been increased drastically. Having in mind that many fields in Eastern Germany are even larger than 50 ha, the 80% overlapping approach is nearly unsuitably for farmers and even for agricultural services.

For adopting UAV to precision agriculture, the processing and analyzing of UAV images should be as easy as possible. However, in this study they were a number of processing steps included that might seem overstated. According to Lelong et al. [35] and Rasmussen et al. [63], UAV imagery should be preprocessed to account for vignetting or BRDF effects. In this study, we used algorithms that were openly described in the UAV literature and relatively simple to implement. Applying these pre-processing algorithms on the images before photogrammetric processing had a visible effect on the final UAV scene (Figure S4). Therefore, we decided to use the pre-processed UAV scenes for modeling. Furthermore, we extracted the plant pixel from the sample plots to relate only those with the crop parameters at first. Of course, plant separation as well as computation of the coverage would need some additional GIS analysis in order to work properly for mapping the entire image (field). Therefore it would be preferably to skip this process altogether in order to oblige to a more simplistic image processing. To prove this we extracted all image variables without segmentation on the basis of all pixels within the sample area. The average correlation strength between all image variables and the crop parameters did not change drastically. Therefore we decided to skip the plant separation step for the image variables and used it only for calculating the CVR variable (coverage).

It has to be acknowledged that the image variables were strongly interrelated with each other (Table 4). This might influence the linear regression due to the influence of multicollinearity. To avoid this problem, one can transform the variables into an alternative space in which the transformed variables become uncorrelated with each other. This can be avoided by using PCA. It is also a useful tool to explore relationships between variables in a more holistic way since most of the variance of the image variables is summarized within a few principal components (PC). To some extent, PCA is affected by strong non-linear behavior of the input variables. In our case, we do not see severe non-linearities among the image variables itself that would justify a more sophisticated non-linear PCA analysis to reduce the dimensionality such as kernel PCA. In order to back up the linear dimension reduction approach, we used the Kayser-Meyer-Olkin criterion and calculated the measure of sampling adequacy (MSA). The MSA value for the complete correlation matrix was 0.71. This is classified as middling (quite good). No MSA value was lower than 0.6 [64]. Of course, image variables and crop variables tend to have non-linear relationships in later growths stages. To some extent, this was reduced by summarizing the image variables via PCA into components (Figure S5).

5. Conclusions

With low-cost UAV imagery based on an RGB consumer-level camera, we were able to compute ultra-high resolution orthoimages and surface models from a cultivated 11 ha wheat field even by using only a 60% image overlap. We observed the development of spatial patterns in the wheat canopy from booting to grain filling with three flight missions and were able to qualitatively assess the changes in the wheat canopy down to the individual plant level.

Various biophysical parameters, i.e., the leaf area index, fresh biomass, dry biomass and plant height, were highly correlated with the blue-channel-related ratios, the plant height calculated from the surface models and the plant coverage calculated from thresholding the UAV imagery. Both the image tone and the surface models derived from the UAV imagery were important to describe the spatial variability of the biophysical characteristics observed in the wheat canopies. Linear regression models for these parameters with principal components calculated from the image variables yielded R² values between 0.70 and 0.97 for the entire data and 0.73 and 0.99 for the validation. Our study revealed that even under water-deficit conditions, UAV imagery can be used to estimate biophysical parameters properly. For example, deriving biomass or LAI maps from UAV imagery would help to establish irrigation systems more properly or improve plant disease forecasting.

Modeling the N content yielded only low R² values, with the best model obtained at flowering growth stage, with an R² of 0.65. The red image tone variables were the most important predictors because the red signal is influenced by chlorophyll absorption and chlorophyll can be related to the N content. During the later missions, the plant height and coverage derived from the UAV images became important because with further propagation of senescence, chlorophyll content decreases in the wheat leaves. When calculating the N status with the NNI, we obtained a constant error over all missions, with 0.10 to 0.11 NNI. However, the errors obtained here were comparable to other studies conducted with hyperspectral field measurements. The R² values might be underrepresented due to the low N variability because of the soil-induced water deficit stress conditions. Therefore, the prospects for low-cost RGB UAV imagery might be limited for wheat crops under these conditions to observe the spatial variability of N and NNI because the RGB image tone was only broadly related to changes in N during the flowering growth stage. However, the combination of morphological information, calculated from surface models, and the image tone is sensible for observing the NNI because it is dependent on the biophysical and biochemical characteristics of the wheat crops. Future studies should be focused on the relationship between UAV imagery and NNI under more variable N conditions using low cost RGB cameras and hyperspectral UAV sensing should investigate wheat crops under water-deficit conditions.