Regression Kriging for Improving Crop Height Models Fusing UltraSonic Sensing with UAV Imagery

A crop height model (CHM) can be an important element of the decision making process in agriculture, because it relates well with many agronomic parameters, e.g., crop height, plant biomass or crop yield. Today, CHMs can be inexpensively obtained from overlapping imagery captured from unmanned aerial vehicle (UAV) platforms or from proximal sensors attached to ground-based vehicles used for regular management. Both approaches have their limitations and combining them with a data fusion may overcome some of these limitations. Therefore, the objective of this study was to investigate if regression kriging, as a geostatistical data fusion approach, can be used to improve the interpolation of ground-based ultrasonic measurements with UAV imagery as covariate. Regression kriging might be suitable because we have a sparse data set (ultrasound) and an exhaustive data set (UAV) and both data sets have favorable properties for geostatistical analysis. To confirm this, we conducted four missions in two different fields in total, where we collected UAV imagery and ultrasonic data alongside. From the overlapping UAV images, surface models and ortho-images were generated with photogrammetric processing. The maps generated by regression kriging were of much higher detail than the smooth maps generated by ordinary kriging, because regression kriging ensures that for each prediction point information from the UAV, imagery is given. The relationship with crop height, fresh biomass and, to a lesser extent, with crop yield, was stronger using CHMs generated by regression kriging than by ordinary kriging. The use of UAV data from the prior mission was also of benefit and could improve map accuracy and quality. Thus, regression kriging is a flexible approach for the integration of UAV imagery with ground-based sensor data, with benefits for precision agriculture-oriented farmers and agricultural service providers.


Introduction
In modern agriculture, fertilization and plant protection can be highly optimized if spatial data about the crop canopy can be made available in a timely manner, and integrated into agronomic measures [1,2].The crop height of a canopy permits statements about the growth stage, crop health, plant biomass, and crop yield [3,4].Therefore, crop height models (CHM) estimated for entire fields can be an important element for the success of any site-specific management strategy.A number of sensors and principles have been tested and are already adopted by the farmers or agricultural services that would enable the construction of CHMs.One of the most cost-effective ways of determining CHMs for precision agricultural needs is through the use of consumer-grade cameras on unmanned aerial vehicles (UAV) [5,6] and by using distance sensors attached to vehicles that pass the fields for making agricultural measurements [7].
UAVs have recently become a very versatile sensor platform for obtaining crop information in fields for the following reasons [8,9]: UAV systems have been made consumer-friendly and cheap, often as an all-in-one solution combining a copter, camera, and flight assistance system.Today's photogrammetry software uses fast algorithms based on structure from motion techniques to obtain 3D information from large sets of overlapping 2D images to generate orthoimages and surface models in reasonable time [10,11].In contrast to satellite-based remote sensing, UAVs deliver images even under cloudy conditions, making its use available most of the time during the season.Due to the high customizability of the UAV system, e.g., flight altitude, camera lens, or mission planning, adaptations of many agronomic problems are possible [9].To construct crop height models from UAV images (UAV-CHM), the difference method can be used.This leads to better results than with estimating the crop height from one 3D point cloud with dense cloud classification directly, because there is a high uncertainty in classifying pure ground pixels [12].The difference method needs two surface models to be generated.These are subtracted from each other, one from the crop surface and a corresponding one from the bare ground, for example before sowing or after harvesting.The accuracy of this approach has been tested by Bendig et al. and Grenzdörffer during field trials [5,12].Subsequent studies showed the good relationship of CHMs in combination with the color information of the orthoimage with biomass [5,6,13,14].
Another way to obtain CHMs is by using proximal sensing.In precision agriculture, proximal sensing has been used to collect data about the plant canopy or soil variability over the entire field, with vehicle-based sensors very near to, or within the canopy or soil, in order to adapt fertilization or crop protection [15][16][17][18].Next to laser-based sensing (LiDAR), ultra-sonic (ULS) sensing is one of the direct ways to estimate the distance between crop surface and sensor.In principle, the sensor emits a sound, the sound wave hits an object and is reflected, and the returning echo is recorded at the sensor.As long as the speed of the emitted sound in the air medium is known, the distance to the object can be calculated by the time difference between sound emission and echo recording.When calibrated to the bare ground, the crop height is yielded from ULS sensing of the canopies' surface.The first study that investigated ULS sensing for determining crop height was conducted by Shibayama et al. [19].They showed that echo intensity mainly depended on plant components oriented perpendicular to the wave burst.Subsequently, it was shown that ULS sensing delivers reliable information for crop height and weed height within field studies [20,21].Other studies showed that ultrasonic sensing can also be used to estimate crop density and biomass [22,23].ULS sensing has also been effectively integrated on phenotyping platforms for plant breeders [24,25].However, despite the good relationship with several agronomic parameters that can be obtained by ultrasonic transducers, research is still underrepresented compared with studies evaluating laser scanning or optical sensors [26,27].While the UAV sensor platform allows free positioning over the canopy, ULS sensors are proximal sensors.Their signal attenuates fast and the distance measurement ability is mostly below 10 m [28].These are best attached in short proximity over the canopy at the farm machinery.Therefore, ULS crop height measurements are bound along the tracks the fertilizer spreader or field sprayer will take as it travels over the field.Due to the low altitude above the canopy, the measurement have a small field of view, so that large measurement gaps perpendicular to the driving direction must be interpolated to construct CHMs for the entire field.
UAV based and ULS based derivation of crop height has obvious pros and cons.This leads to the idea of combining them together to improve the CHM and its applicability for estimating agronomic parameters.In general, the fusion of data from different sensors was favorable and improved the estimation of agronomic parameters [29,30].Hybrid geostatistical interpolation regression kriging would be a suitable choice because it can be assumed that data from both approaches are correlated with each other and inherit a spatial structure.Regression kriging combines regression and kriging by first conducting an arbitrary regression with a covariate and then interpolating the residuals of the resulting model with kriging [31][32][33].Regression kriging was used intensively as a spatial modeling tool in soil science and environmental science when there was a sparse variable sampled alongside correlated exhaustively available variables to improve the spatial mapping of the sparse variable.Case studies involved the mapping of soil properties with remote sensing data, geological data or terrain attributes [34,35].Specifically, in precision agriculture, soil electrical conductivity, a parameter that can be easily obtained from the bare soil with proximal soil sensors, was often used with regression kriging for improving the spatial estimation of various soil properties within the field [36,37].The use of regression kriging with UAVs is rare.One study used universal kriging in combination with model predictive control to track and predict oil spills boundaries with UAV [38].Using regression kriging in combination with UAV imagery for determining crop height has not been investigated so far.
Can we enhance the performance of ULS sensing for estimating crop height models for wheat fields with data derived from UAV photogrammetry?This data fusion approach makes sense because ULS sensors are installed on vehicles driven on the ground, and measured data is therefore bound to locations along specific traffic lanes in the field, whereas UAV data comprises the whole field continuously.Our assumption was that if UAV data was integrated as a covariate within the framework of a hybrid interpolation model, namely regression kriging, as a kind of data fusion approach, the accuracy of the CHM based on ULS measurements can be increased and the relationship with agronomic parameters could be improved.It was further tested whether earlier acquired UAV data in the growing season might improve the CHM as well.

Study Sites
The study was conducted in two fields (A and B) located in Eastern Germany during the spring seasons of 2015 and 2016 (51 • 49 N, 12 • 42 E).The investigated field sizes were for field A and B, 9.6 ha and 7.1 ha.The soils in both fields are influenced by remnants of the recent floodings of the nearby Elbe River.The main soil type is a fluvic cambisol and the soils are spatially variable in soil texture, ranging from sand and loamy sand to sandy loam.The crop type in both fields was winter wheat (Triticum aestivum L.) with the variety 'Linus' in field A and 'Montana' in field B, with an average crop density between 440 and 480 ears m −2 and a seed row distance of 0.12 m.The crop rotation is generally cereal dominated.Ground water influence contributed to crop variation.

Sensor and Sensor Platform Specifications
As a UAV platform, the hexacopter P-Y6 (Hexapilots, Dresden, Germany) with the navigation control system DJI Wookong M (DJI Innovations, Shenzhen, China) was used with an integrated Global Navigation Satellite System GNSS receiver to enable user-defined waypoint travelling for flight mission planning.A Sony Nex 7 camera was mounted on a Gimbal underneath the copter, and positioned to receive nadir view.The specifications of the camera were: 24 megapixels, 23.5 × 15.4 mm sensor, and E16 mm F2.8 fixed lens (Sony Corporation, Tokyo, Japan).The images were acquired from a flying altitude of 50 m corresponding to a ground resolution of approx.0.012 meter per pixel.Images were recorded in raw and lossless converted to Tagged Image File Format TIFF.
As an ultrasonic sensor, the commercial sensor device mic-340/D/M (Microsonic GmbH, Dortmund, Germany) was used.The sensor is built within a stable metal design for harsh outdoor environments.Its sensor range lies between 350 and 5000 mm, and its transducer frequency is at 120 kHz.Due to internal temperature compensation, the sensor produces reliable measurements between −25 and 70 • C outdoor temperature.The analogue output voltage signal was converted to digital values with Bus-Powered M Series Multifunction Data Acquisition (DAQ) (National Instruments Corp., Austin, TX, USA) on a 16-bit basis.The ultrasonic sensor was mounted on a tripod arm connected to the side of a tractor 2 m distance from the bare ground, and 2 m distance from the tractor.This ensured undisturbed top (nadir) view from the sensor to the wheat canopy.The ULS measurement has a circular footprint on the ground with a radius of approximately 1 m.Conversion to height values was made possible by prior calibration of two different voltages received from echoes reflected from two different measurement heights relative to the ground (1 m and 2 m).Differential global positioning system (GPS) measurements were recorded using a Trimble AgGPS with ±1-2 m accuracy to locate the current position during field passage.DASYLab software (measX, Mönchengladbach, Germany) was used for direct signal processing, conversion to height values, and synchronization with GPS measurements.

Field Measurements
Two measurement campaigns (missions) were conducted in field A (2015) and in field B (2016) between mid-May and the end of June.UAV and ultrasonic measurements (ULS) were both acquired on the same day.Additionally, one UAV mission was conducted over bare soil after harvesting, to determine the ground surface model.The flight missions were planned with the DJI software in such a manner that image capturing was triggered with a frequency to ensure a 60% overlap along and parallel to the flying direction.For georeferencing, 20 reference panels were allocated as ground control points (GCP) clearly visible from above within the traffic lanes and along a rectangular grid spanned over the field.ULS measurements were taken with a tractor speed of 8 km•h −1 and with a recording time interval of 1 s, resulting in a measurement density along track with 2 m/point.Before each measurement campaign, the null level of the ultra-sonic sensor was re-calibrated to the bare ground.Measurements were calibrated with one plant height measurement taken before mission start.
For ground truthing, 20 measurement plots were arranged along two different traffic lanes of the ULS sensor, with a dimension of 1 × 1 m 2 .These measurement plots were visually located in such a manner as to represent overall crop height variability within the fields.From these measurement plots, crop height and biomass were determined.For crop height, the heights of 10 individual wheat plants were randomly measured with a meter ruler and averaged.For biomass, all wheat plants within the plot area were harvested by using a short reaping hook.Fresh biomass was immediately measured after cutting.Measurement plots and GCPs 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.
The wheat crops were harvested on 31 July 2015 and 21 July 2016 for fields A and B, respectively.The crop yield was recorded with a yield monitoring system installed on a combine harvester along specific tracks in each field (T550, John Deere, Moline, IL, USA).Only one harvester was used for harvesting these tracks to prevent interference and systematic errors.The software GreenStar 2 was used to extract single crop yield values along with their GPS position.

Generating the Crop Height and Orthoimage Data from the UAV Images
The UAV crop height models (CHM) and orthoimages were generated from photogrammetric processing of the aerial images that were captured during the UAV flight missions.A number of different pre-processing steps were employed on the UAV images prior to photogrammetry to improve the photogrammetric results.This included vignette correction and contrast enhancement [39,40].These steps are described in full detail in Schirrmann et al. [17].To generate the surface models, the pre-processed UAV images were used in the photogrammetry software Agisoft PhotoScan (v.1.2.4,Agisoft LLC, St. Petersburg, Russia).We used Agisoft Lens (v.0.4.2,Agisoft LLC, St. Petersburg, Russia) to estimate the calibration coefficients of the used camera sensor, instead of using the empirical, internal estimation of the camera coefficients from the images itself by PhotoScan.For the first alignment of the images, these calibration coefficients as well as the GPS coordinates from the images were fundamental.The alignment was done by detecting local, matching feature points among the overlapping images, and following the movement of those points throughout the scene using advanced structure-from-motion algorithms [41].The image alignment was then further optimized by using the GCPs.Afterwards, the dense point cloud was generated from the GCP optimized image alignment and a full multi-view 3D stereo-reconstruction was performed.Thus, the GCP were fully integrated in the structure using the motion process flow in PhotoScan.The CHMs were generated with a reprojection error between 0.8 and 1.3 pix for field A, and 2.0 and 3.0 pix for field B.
CHMs and orthoimages were resampled with bilinear interpolation on a 1 × 1 m grid.CHMs were normalized by their mean values.The orthoimage information was calculated based on the following ratio: where R, G and B were the red, green and blue image color information, and the constant a was set to 0.667 [42].The VEG ratio showed high effectiveness for separating plant information in wheat crops [43].It was used as a covariable in addition to CHM to improve crop height estimation and is termed "ortho" throughout this paper.

Kriging Interpolation
The ULS measurements were interpolated with ordinary block kriging on a 1 × 1 m grid spanning the entire fields with correspondence of the UAV-CHM and UAV-ortho grids [44].To apply ordinary kriging, experimental semivariograms γ(h) were first calculated from the ULS data given by: where n(h) is the number of point pair comparisons between the points z(x i ) and z(x i + h), which are separated by the distance (lag) vector h [45].On the basis of the experimental semivariograms, semivariogram models were fitted with a weighted least squares approximation and a specific semivariogram function, e.g., a spherical or exponential function.The semivariogram model was then used for ordinary kriging to estimate the ULS-CHM [44]: where λ denotes the kriging weights solved from the kriging system of normal equations using the semivariances, z(x i ) the measurement points, and ẑ(x) the new prediction point.
To include UAV data in the spatial estimation of the ULS-CHM, regression kriging was used.Regression kriging is a multivariate extension of the kriging interpolation that allows for including additional spatial information when it is related to the target variable.In regression kriging, a prediction is estimated by the sum of two components originated from the regression between the target variable and related independent variables [31].
The first sum in Equation ( 4) describes the fitted trend estimated with a generalized least-squares regression model (GLS) with βk being the estimated drift model coefficients and q p being a set of predictors.The second sum describes the kriging prediction of the residuals e(x i ) based on the kriging weights λ, which were solved from the semivariogram model of the residuals.The drift model was calculated by regressing the ULS measurements on the UAV-CHM and UAV-ortho data.From the residuals of this model, the experimental semivariograms were calculated given Equation ( 2), and the semivariogram models were estimated in order to use them for regression kriging.
Geostatistical analysis, e.g., semivariogram modeling, ordinary kriging and regression kriging were performed with the gstat package [46] using the statistical programming language R Project for Statistical Computing [47].

Data Analysis
For all fields and missions, spatial estimations of crop height using ULS measurements with ordinary kriging and regression kriging with UAV-CHM and UAV-ortho data were computed and compared with the ground truthing data.We also tested if the crop height acquired from a later ULS mission could be improved with UAV data acquired from an earlier mission.For this, ULS measurements from Mission 2 were interpolated with UAV data from Mission 1 as covariate with regression kriging.
To compare ULS measurements, UAV data, and ground truthing directly, the ULS measurement track went through the ground truthing plots so that we had to simulate the gap-like character of the ULS measurements before comparing regression kriging and ordinary kriging-generated ULS-CHMs with each other.Thus, spatial buffers around the ground truthing plots were generated with a radius of 6 m, and ULS measurements falling within these spatial buffers were deleted.Linear regression models were calculated using crop height and fresh biomass as response.In order to compare the performance of the data as well as the different interpolation methods, we reported the root mean square error (RMSE) of the residuals and the r 2 from the models: with y i being the data or observed values, ŷi being the model predictions, and y being the mean of the data or observed values.We furthermore investigated the relationship between crop yield and ULS-CHM generated by ordinary kriging and regression kriging with a Spearman rank correlation.Differences between correlation coefficients were investigated with Bonferroni-adjusted confidence intervals computed by bootstrapping (n = 1000).

ULS and UAV Measurements and its Relation with Crop Height
The ULS measurements are shown in Figures 1 and 2, which superimpose the relative UAV-CHM data derived from the UAV imagery.The ULS sensor was attached to a ground vehicle, so that the sensor measurements were confined to specific traffic lanes within the fields.When considering only the ULS measurements, the values were agreeable between the nearest positions of the neighboring tracks, and patterns over all tracks were recognizable in each of the four missions.Thus, between-track errors were assumed to be relatively small, which underlined the consistency of the ULS measurements.The relationships between the ULS measurements and manual plant height measurements at the ground truthing plots were linear and highly significant (p < 0.001).The error deviations from the regression lines were between 2.50 and 4.55 cm, with an r 2 ranging between 0.71 and 0.95 (Table 1).The ULS measurements were a further highly significant predictor for fresh biomass.However, the relationship was slightly less strong than crop height, except in Mission 1 of field A. The UAV-CHMs were generated from the photogrammetric processing of the overlapping image sets acquired during the missions.Thus, UAV measurements were available on a raster grid covering the whole fields and not confined to individual measurement tracks.However, due to some problems with photogrammetric processing with the data of Mission 1 in field B, there were some missing areas in the eastern part of the field because images could not be aligned properly here.In general, the UAV-CHM data showed the crop height with a lot of detail (Figures 1 and 2).Even the narrow traffic lanes within the fields were visible in the CHM data.The relationships between the UAV-CHM and manual plant height measurements at the ground truthing plots were linear and highly significant (p < 0.001), but the error deviations in terms of the RMSE were mostly higher compared to the ULS measurements, and ranged between 3.13 and 10.43 cm with an r 2 between 0.56 and 0.96 (Table 1).For field A, the UAV-CHM was more in line with the manual crop height measurements than in field B. Specifically, Mission 2 of field B had a very high RMSE.The UAV-ortho data in terms of the VEG ratio were also in relation with crop height.In fact, in some missions, the ortho data had an even stronger relationship with crop height than the UAV-CHM data.In Mission 1, they even outperformed the ULS measurements.The UAV data were also a good predictor for fresh biomass but it did not outperform the ULS measurements.Specifically, the UAV-ortho data were a better predictor than the UAV-CHMs for fresh biomass.As can be seen from the figures, the ULS and UAV data had a quite high resemblance to each other and, in most cases, sudden changes of the UAV-CHM could be reproduced in the ULS measurements.They were significantly and linearly correlated at similar positions with an r between 0.62 and 0.89 (p < 0.001) in field A and 0.49 and 0.84 (p < 0.001) in field B (Table S1).Even the form of the experimental semivariograms of the ULS and the UAV data had a high similarity, and all semivariograms were well fit with an exponential model (Figure S1).The nugget-to-sill ratio of all semivariogram models were below 35%, which indicated high spatial dependence between the measurements, and a low contribution of stochastic measurement error.Due to its strong spatial structuring, choosing geostatistical methods for interpolation should be reasonable and beneficial because the spatial structure inherent to the data sets can be well implemented with a semivariogram model in the interpolation process.

Regression Kriging of the ULS Measurements
One of the major concerns of measurements acquired with ground vehicles during the growing season is that they are confined to specific tracks within the field.Thus, these are data sets which have spatial gaps that have to be filled with a valid interpolation model.As we learned from the semivariograms in the previous section, ULS measurements had a strong spatial dependency with each other.Thus, interpolating the measurement points would be desirable.However, choosing a univariate approach, such as ordinary kriging, as the interpolation method, would lead to very smooth prediction surfaces lacking the detail of a natural environment.We also saw that the ULS measurements had a significant linear relationship with the UAV-CHMs.Since the latter are normally available for the entire field, it seemed sensible to apply regression kriging, which allows the inclusion of additional spatial information in the interpolation process.
To apply regression kriging, we first computed the semivariograms of the residuals which resulted from the linear model of the ULS measurements regressed on the UAV data.The residuals of the linear models had a smaller degree of spatial structuring than the ULS measurements themselves, as indicated by their larger NSR or the smaller sill that the model reaches (Table 2, Figures S2 and S3).This effect arose because some spatial autocorrelation has been already explained by the UAV data in the linear model.For example, the high NSR and the low sill reached by the model, and the short autocorrelation range of the UAV-CHM semivariogram model for Mission 2 in field A, were the direct result of the strong relationship between the UAV-CHM and ULS measurements (r = 0.89).
For comparing the performance of regression kriging with ordinary kriging, the predictions were regressed with the manual crop height measurements (Table 2, Figure 3).What we can glean from this table is that regression kriging outperforms ordinary kriging in all missions in this setting.Even for the weakest relationship between ULS and UAV data (mission 2 of field B), the scattering around the regression line was reduced, and a gain in prediction accuracy was obtained by regression kriging.This gain was obtained for both, the crop height and fresh biomass.The latter benefitted even more strongly from regression kriging.The degree of outperformance was highly linearly related with the degree of linear correlation between ULS and UAV data (Table S1).
When comparing the visual results of the ULS-CHMs generated either by ordinary kriging or by regression kriging, the ordinary kriging predictions were very smooth and lacking details, whereas the regression kriging predictions were rich in detail, with many sharp boundaries (Figures 4 and 5).The missing part in the UAV map of Mission 1 of field B was filled with spatial information after  Since crop height is one of the better predictors for estimating crop yield, we also investigated if its correlation with crop yield could be improved with regression kriging.In the scatter plots of the Figures 6 and 7, the predictions with ordinary kriging using only ULS data, and with regression kriging using UAV-CHM and UAV-ortho data to enhance the ULS-CHM, were shown in relation with crop yield values measured by a yield monitor.All predictions were significantly related with crop yield.Generally, the predictions of the second mission had a higher correlation with crop yield than the first mission due to lesser influences of growth effects, weather effects and other environmental factors.Regression kriging predictions were slightly better or equally well correlated with crop yield than using the ordinary kriging predictions.Specifically, higher crop yield (>70 dt/ha) was significantly stronger correlated with regression kriging predictions than for ordinary kriging prediction for Mission 2 and field A (from r = 0.35 to 0.61) and B (from r = 0.23 to 0.32) on p < 0.05, although the global correlation strength was only slightly different in field A (from r = 0.89 to 0.93) or even not significantly different in field B (from r = 0.60 to r = 0.66) with p < 0.05.

Discussion
Despite its plain and inexpensive design and implementation, the ULS sensor provided amazingly stable measurements over the large fields.Post-processing of the measurements, e.g., cleaning the data sets from erroneous measurements or tracks, was not necessary, and we were able to use the data sets "as is".In comparison with other studies that investigated proximal sensors for estimating crop height, the ULS measurements showed a very good performance with an RMSE below 6.0 cm.In the study of Scotford and Miller [20], the crop heights for three different varieties of winter wheat were measured with an ULS sensor and they reported an error between 4.6 and 7.2 cm.However, the slope of the models for regressing the crop height on the ULS measurements slightly differed between the missions and there was a substantial underestimation of the true plant height values.Shibayama et al. [19] showed that the relative frequency of detected echo peaks was not only dependent on plant height, but also on leaf angle, leaf area and leaf density.The integrated ULS signal was therefore also substantially influenced by the leaves within the canopy, which led to an underestimation of the true plant height.This may explain the underestimation and the slight differences of the model slopes but this may also explain the quite high correlation with fresh biomass with r 2 values between 0.76 and 0.90 (p < 0.001).Azis et al. [48] investigated ULS sensing of plant height on 18 plants under laboratory conditions and also yielded different coefficients of the linear regression models when models from different growth stages were compared with each other.In our study, we could not show that the use of UAV data within the regression kriging setting could decrease the effect of underestimation of the ULS CHM to the true plant height.However, the underestimation as well as the slope differences could be easily circumvented by using several plant height measurements as calibration values, because these measurements are non-destructive, inexpensive, and easily obtain from the field.
The accuracies of the UAV data for estimating crop height were in close proximity with other studies, e.g., an RMSE between 3 and 9 cm [30] or an r 2 between 0.60 and 0.72 [49].However, they were generally worse than the ULS measurements for directly estimating crop height.There were many factors that influenced the quality of the UAV data and may have led to a reduced accuracy.The copter was not a fixed sensor platform.It could be tilted and twisted through wind effects during the missions so that the camera lens may slightly diverge from nadir positioning.Substantial difference in camera position due to shaking of the copter may result in blurred images.However, in our case, the acquired images were rarely rejected from the photogrammetric software because of blurring.Also, the copter recorded the images at an altitude of 50 m above ground to collect the data for an entire field in an economically reasonable way.Thus, limitation of vertical resolution must be taken into account.Furthermore, the crop height values were generated from the difference of two photogrammetric processed height models.This indirect approach is not error free and cannot compete fully with the rather direct estimation of crop height using ULS sensors.For the 3D reconstruction of the height model it is necessary that the software aligns the overlapping image set properly.This process requires finding identical points in overlapping images within the image set.However, large areas in wheat fields can be quite uniform.Images acquired from those areas contain less specific information that can be extracted via key point extractors for identifying equal points in the overlapping images.Specifically, Mission 2 in field B was problematic.Although we increased image overlapping to 80% in Mission 2 and distributed additional signal plates over the field, alignment and 3D reconstruction results were worse than for the other missions.We concluded that the overall homogeneity that we observed in this field at this time negatively influenced the photogrammetric process flow in AgiSoft.
As we have seen, both approaches have their limitations.Thus, combining them with a hybrid method to improve the spatial mapping results may overcome some of them.Regression-kriging was able to fully exploit its potential here, because we had a sparse data set (ULS) and an exhaustive data set (UAV), and both data sets had favorable properties for geostatistical analysis, e.g., high number of measurement points, sufficient number of measurement points in small distance lags, and strong spatial structures as was visible within their semivariogram plots, while they had a significant relationship with each other over the entire fields.The strong relationship between ULS and UAV data was captured by regression kriging as a trend that could improve the local estimation of crop height, because at any estimation point, UAV data was available.This produced a quite strong visual difference from the rather smooth map computed by ordinary kriging, to a much more detailed map which resulted from regression kriging.It also produced maps with higher accuracy.This is in agreement with other studies comparing regression kriging with ordinary kriging if a strong relationship between main variable and covariate exists [34,36,50].We also corroborated that the improvement in accuracy and also improved modelling or correlation, with the agronomic parameters of biomass and, to a lesser extent, crop yield.The correlation between crop yield and ULS measurements were higher than correlations which were determined in soy bean fields in Bai et al. [25].However, although a good correlation between crop yield and ULS measurements could be stated, there was still a quite high uncertainty, especially for the more uniform field B. We could determine an increase of correlation with high crop yield.This might have been the result of the more detailed map regression kriging produces with UAV as covariates, which in turn takes extremes more strongly into account.In contrary, ordinary kriging leads to a smoothing that can destroy vital information important for precision agriculture.A better prognosis of subareas within the field where high crop yield is expected to be of high interest for the farmer because these are areas where special attention should be given for tolerance to plant disease or nutrition deficiencies, in order to minimize economic impacts.Furthermore, the estimated spatial maps with regression kriging did not show any artifacts that may arise due to an uneven physical relationship between covariates and the main variable [31].Regression kriging draws its power from seamlessly shifting between regression and kriging, from estimating the trend/drift to estimating residual information.If the relationship between the main variable and the covariate is weakening, the more intense the kriging part will become.For that reason, we could make use of the UAV data from prior missions.Despite different features were visible in the prior UAV maps, the features based on the ULS measurements were most prominent in the regression kriging maps.

Conclusions
The integration of UAV imagery with regression kriging for the spatial estimation of crop height models with ultra-sonic sensing successfully improved relationships with crop height, fresh biomass and, to a lesser extent, with crop yield, compared to univariate interpolation ordinary kriging.It was apparent that the maps generated by regression kriging were of much higher detail than the corresponding maps estimated with ordinary kriging, because with regression kriging at any estimation point, UAV data was available.The UAV imagery that was acquired earlier in the growing season still contributed to improve the spatial estimation of crop height using regression kriging.Thus, there is no need for immediate action for conducting the UAV mission as soon as the ULS measurements are acquired but this can be relaxed over a longer period during the growing season.Regression kriging turned out to be a very flexible data fusion approach for UAV imagery in combination with ground-based sensor data.Its use may help to improve decision-support in agriculture, specifically for farmers and agricultural services with a focus on precision agriculture applications.

Supplementary Materials:
The following are available online at www.mdpi.com/2072-4292/9/7/665/s1, Figure S1.Frequency distributions and semivariogram models for the ULS measurements, UAV-CHM and UAV-ortho data for field A and Mission 1 (a, b, c) and Mission 2 (d, e, f) as well as for field B and Mission 1 (g, h, i) and Mission 2 (j, k, l). Figure S2.Semivariogram models of the residuals from the ULS linear regression models for field A with UAV-CHM and UAV-ortho as independent variable for Mission 1 (b, c) and Mission 2 (e, f) in comparison with the ULS semivariograms (a, d).For g and h, the UAV-CHM and UAV-ortho of the prior mission 1 was used.Figure S3.Semivariogram models of the residuals from the ULS linear regression models for field B with UAV-CHM and UAV-ortho as independent variable for Mission 1 (b, c) and Mission 2 (e, f) in comparison with the ULS semivariograms (a, d).For g and h, the UAV-CHM and UAV-ortho of the prior mission 1 was used.Table S1.Correlation matrix of the sensor variables.

Figure 1 .
Figure 1.The unmanned aerial vehicle (UAV) crop height models superimposed with the ultrasonic (ULS) sensing measurements for Missions 1 and 2 (a,b) of field A.

Figure 2 .
Figure 2. The UAV crop height models superimposed with the ULS measurements for Missions 1 and 2 (a,b) of field B.

Figure 4 .
Figure 4. Ordinary kriging and regression kriging-interpolated ULS-CHMs for Mission 1 (a,b) and Mission 2 (c,d) using the UAV-CHM as covariate, as well as (e) the regression kriging interpolated ULS-CHM of Mission 2, using the UAV-CHM from Mission 1 as covariate for field A.

Figure 5 .
Figure 5.With ordinary kriging and regression kriging interpolated ULS-CHMs for Mission 1 (a,b) and Mission 2 (c,d) using the UAV-CHM as covariate as well as (e) the regression kriging interpolated ULS-CHM of Mission 2 with using the UAV-CHM from Mission 1 as covariate for field B.

Figure 6 .
Figure 6.Scatter plots and correlation (Spearman rank) between crop yield and ordinary kriging (a,d) and regression kriging with UAV-CHM (b,e) and UAV-ortho (c,f) as covariate for field A. For (g,h), the UAV-CHM and UAV-ortho of the prior to Mission 1 was used.

Figure 7 .
Figure 7. Scatter plots and correlation (Spearman rank) between crop yield and ordinary kriging (a,d) and regression kriging with UAV-CHM (b,e) and UAV-ortho (c,f) as covariate for field B. For (g,h), the UAV-CHM and UAV-ortho of the prior Mission 1 was used.

Table 1 .
Linear relationship and root mean square error (RMSE) of ULS measurements and UAV data (crop height model (CHM), orthoimage (ortho)) with crop height and fresh biomass (FM).Autocorrelation range and nugget-to-sill ratio (NSR) of the semivariograms for ULS measurements and UAV data (CHM, ortho).