Analyzing and overcoming the eﬀects of GNSS error on LiDAR based orchard parameters estimation

.


Introduction
Precision agriculture is benefited from the use of communication and information technologies, for enhancing the decision-making process and for automation of a wide range of agricultural tasks.In particular, the use of new sensors and sensing technologies improved the capabilities of the farmer to take preventive actions, such as in the case of weed monitoring and plague detection (Underwood et al. 2017).
Among the new sensing technologies, LiDAR (light detection and range) sensors are used mainly for geometric characterization of groves (Rilling et al. 2017;Trochta et al. 2017;Rosell and Sanz 2012;Li et al. 2018).From the geometric parameters of plants, the crown surface area and volume are of special interest since both combine the width, height, geometric shape and the structure of trees (Auat Cheein et al. 2015).Such parameters are commonly used for farmers for, e.g., stablish a biomass model for plants (Lin et al. 2017), herbicide management and pruning directives (Rosell et al. 2009).
Geometric parameters estimate from raw data are mainly affected by the uncertainties in the laser scanner, calibration of the setup and the georeferencing system.Different works have analyzed the error propagation in scanner laser systems (Mezian et al. 2016;Goulden and Hopkinson 2010;Hartzell et al. 2015).However, those analyses do not consider the GNSS positioning error.Similarly, in the analysis of tree volume sensitivity developed in (Palleja et al. 2010), the GNSS positioning error is not considered.Such work considers specific errors such as vehicle speed, the height of the LiDAR, distance of the measurement and the angular orientation of the LiDAR.
As shown in the references mentioned above, GNSS positioning error is not considered when building the LiDAR-based point cloud.This fact occurs although the magnitude of the navigation system positioning errors has a direct impact on the resultant point cloud, which is not dependent upon any of the other uncertainties (Glennie 2008).Instead, the point clouds merging relies on the ability of the programmer (Rosell and Sanz 2012;Rosell et al. 2009) or on commercial software (Sanz et al. 2018) where no considerations are taken into account regarding GNSS absolute errors.Therefore, error propagation, from the GNSS antenna to every point from the point cloud, is disregarded.
Commonly, LiDAR-based point cloud information (Gaulton and Malthus 2010;St-Onge et al. 2008) is used for mapping approaches.However, when the LiDAR is mounted on a vehicle, such information can also be used in scan matching techniques for vehicle localization (Grant et al. 2019;Malavazi et al. 2018).The goal of those approaches is to find the position of the sensor over time by finding the rigid transformation matrix which maximizes the overlap between point cloud scans (2D or 3D) obtained at different times.
The process of finding the rigid transformation matrix is known as registration.
There are several variations of registration algorithms, for example, the Normal Distribution Transform (Magnusson et al. 2009), which creates a voxel grid over the point clouds or the Gaussian Mixture Model (Boughorbel et al. 2004), which performs registration based on expectation maximization.
Other algorithms focuses on specific features from the scene to perform registration, such as corners or planes (Lamine Tazir et al. 2018;Peng et al. 2016).A well-known alternative is the iterative closest point registration, which has been widely used in the last years (Ren et al. 2019;Kim et al. 2018;Donoso et al. 2017); therefore, many variations have been proposed, such as EM-ICP (Granger and Pennec 2002) and Generalized-ICP (Segal et al. 2009).All these techniques can provide the pose of a vehicle when the LiDAR is mounted on it.Nevertheless, the position estimation process is constrained by the quality of the point cloud and the efficient matching among points (Zaganidis et al. 2018).
In this work, we first describe the error propagation attached to mobile terrestrial LiDAR measurements, and its effects on the estimation of orchard parameters when only GNSS measurements are considered as localization.
Then, we describe an approach which improves the vehicle localization, by fusing a scan matching algorithm with GNSS measurements.For the experimentations, we increment the GNSS uncertainty up to two meters to quantify its influence on orchard parameter estimation.We analyze the case for the following features: crown surface area, crown volume, and porosity.
To do so, we implement previously published techniques (Pfeiffer et al. 2018;Trochta et al. 2017;Li et al. 2018) and evaluate their sensitivity to changes in the error of the GNSS antenna, for two cases: (i) GNSS as a solo localization system and (ii) GNSS with scan matching as localization system.
The remaining of the paper is organized as follows: Section 2 describes in detail the experimental setup and the methodology followed in this work.
Section 3 presents the experiment developed and the results for the crown parameter estimation.Finally, Section 4 shows the conclusions and future work.

Materials and methods
This section describes the experimental setup and the methodology followed to estimate the orchard parameters.First, the methodology describes the error propagation of GNSS over the laser scan measurements; second, it focuses on describing the proposal to overcome such effect.

Experimental setup
To perform our analysis, we equipped a hydro-pneumatic sprayer with an RTK Leica GPS1200 and a LiDAR Velodyne PUCK; both are assembled in a rigid structure, as shown in Fig. 1.The LiDAR Velodyne PUCK can generate a 3D point cloud (x-y-z positions) of the scanned scene using up to 16 ray readings per frame, with a maximum range of 100 meters and an accuracy of ± 0.03 meters.The RTK-GNSS system GPS1200+(Leica Geosystems AG, Heerbrugg, Switzerland) provides absolute coordinates and UTC time (synchronized with the LiDAR) with a frequency up to 20 Hz and precision of approximately 20 mm.

Error propagation
When registering a point cloud using a LiDAR and a GNSS antenna, there are at least three coordinate frames: a global reference frame, < Global >; a coordinate frame attached to the GNSS receiver, < GN SS >; and a coordinate frame attached to the LiDAR sensor, < LiDAR >, as can be seen in The point P (in spherical coordinates) acquired in the < LiDAR > coordinate frame is first transformed to the < GN SS > coordinate frame and then to the < Global > reference frame.Twenty trees have been used for this experiment, where the distance between each tree in the same row is L=1 meter and the height of each tree is H=2 meters.
Each point registered from the environment has three coordinates.In our case, such three coordinates are spherical: one range ρ and two angles, θ and φ, as shown in Fig. 2 (it is worth to mention that 3D LiDARs might also register a point in cylindrical coordinates, depending on the manufacturer).
The point P acquired by the LiDAR in spherical coordinates is then converted into Cartesian coordinates using the transformation shown in Eq. 1.
It is to be noted that LiDAR readings are not extent of error.If we consider that the spherical coordinates associated with a single point are independent Gaussian random variables, then their covariance matrix will be of the form: where σ 2 ρ is the covariance associated with the LiDAR range reading, and σ 2 θ and σ 2 φ are the covariances of the associated angles; Σ LiDAR is the covariance matrix associated with a single point P <LiDAR> .
The point P <LiDAR> can be transformed into the global coordinate frame according to Eq. 3, considering that the GNSS antenna receives localization information and not orientation (Mezian et al. 2016;Hartzell et al. 2015).
To propagate the error, we followed the guidelines presented in (Auat Cheein and Carelli 2013), therefore, analytically we have that: where ∇ P is Jacobbian matrix of Eq. 1 and the Jacobbian matrix corresponding to the transformation from the GNSS coordinate frame to the global frame is the identity matrix.It is to be noted that all covariance matrices used here are positive semi-definite.
From the previous statements, we can see that if we have two GNSS positioning systems, each one of them with their corresponding covariance matrix associated with their errors, i.e., Σ 1 GN SS and Σ 2 GN SS , and if we previously know that, e.g., Σ 1 GN SS − Σ 2 GN SS 0, where stands for positive semi-definite, then, for the same point P : The later means that the higher the covariance associated with the GNSS antenna, the highest the error propagation.Since the above calculation is for a single point P and considering that a point cloud is a concatenation of single points, then the same reasoning can be applied to a covariance matrix associated with a point cloud.

Overcoming the GNSS positioning error
To overcome the error associated with the GNSS and without the use of external hardware, our work aims to improve the results of the parameters estimation based on the fusion of scan matching estimations and GNSS measurements.The methodology used in this work is depicted in Fig. 3.As can be seen from Fig. 3, the first step in our methodology is the scan matching approach.Although several scan matching approaches can be implemented, in this work, the registration process between the previous scan, Λ k , and the current scan, Ω k , (where k represents sampling time instant) was accomplished with the use of a 3D ICP algorithm (Besl and McKay 1992).

GNSS
Others can be used, but its analysis is out of the scope of this work.As a result of the registration, the corresponded transformation matrix, H Ω , was obtained.
The ICP method is an algorithm used to find the rigid transformation H Ω between a target point cloud Ω k and a reference point cloud Λ k so that the matching satisfies a metric criterion.Ω k and Λ k , needs to be found.First, the ICP compute the closest points between the two point clouds as: where Cp represents a closest point operator, therefore, Y k represents the closest points from Ω k to Λ k .To obtain the rigid transformation matrix H Ω between two sets of points, the ICP minimize the sum of square error between them, as shown below: here, t represents translation and R rotation.When applying the alignment to Ω k the resulting point cloud Ω k+1 can be obtained, as seen in Eq. 7.
Then, the distance d between Ω k+1 and Λ k is obtained by: This process is repeated until d is less than a previous establish value, τ , or when the maximum number of iterations, I max , has been reached.It is worth to mention that this process only warranties the convergence to a local minimum (Yang et al. 2016).
Following the guidelines published in (Manoj et al. 2015), we obtained a close-form covariance of the resultant transformation provided by the ICP algorithm.This close-form covariance does not make any assumption on the noise present in the sensor data and has no constraints on the estimated rigid transformation.Considering Ω k and Λ k , a cost J can be defined as: In Eq. 9, R is a rotation matrix and t is a translation vector.Then, the covariance was obtained as follows: where: Here, (x, y, z) represents the translations of H Ω in the x, y and z coordinates, respectively, and (a, b, c) represents the rotations of H Ω in yaw, pitch, and roll, respectively.Additionally, n is the number of correspondences between the point clouds Λ k i and Ω k i .On the other hand, cov( ) represents the noise present in the (x, y, z) components acquired by sensor readings, as shown in Eq. 11.
The above expressions were used for information fusion purposes for getting a more accurate position estimation.This fusion was done via a Kalman filter fusing the data associated with the GNSS and the scan matching.It is worth mentioning that this algorithm must be initialized, with prior knowledge of the position.For the objective of improving the orchard parameter estimation only with remote measurements, the process model of the Kalman filter was set as an identity matrix.More information about the approach implemented in this work can be found in (Sun and Deng 2004;Caron et al. 2006).The sensor fusion procedure for improving localization is shown in Algorithm 1.
In Algorithm 1, lines of code (1) to (3) use prior knowledge for initialization; lines of code ( 9) to ( 17) show the ICP scan matching approach with the information provided in lines of code ( 7) to ( 9); in lines of code ( 20) to ( 26), the procedure for obtaining ICP covariance is computed.The measurement matrix is defined in line ( 27) and the covariance matrix is defined in line ( 28).The estimated position, Xk , is obtained with a Kalman Filter (line 29) considering the created measurement and covariance matrices.Line (30) transform Ω k to global coordinate system according to the estimated position Xk .Finally, the transformed point cloud, S k , is stored in M (line 32) to create an environment point cloud.
Once the environment point cloud, M, is obtained, we performed the foliar parameter estimation.For this manuscript, we have selected the following orchard features: crown volume, crown surface area and porosity, since they are the most common features extracted from point clouds.In addition, we have implemented, as features extraction procedure, the following computational approaches: convex-hull, for crown volume and crown surface area; alpha-shape, also for crown volume and crown surface area; and voxelization, for crown volume, crown surface area and porosity, as presented in (Pfeiffer et al. 2018).
Algorithm 1: Sensor fusion using scan matching and GNSS readings.

Experimental results
The tests were carried out at an apple tree grove.The platform traversed approximately 20 meters back and forward through an apple tree corridor.
The sampling time of the RTK was set to 0.05 seconds, synchronized with the sampling time of the LiDAR; 150 trials were performed.Since the RTK has an absolute error of 0.01 meters (horizontally), it was used as ground truth.
Then, to test the sensitivity of the grove features (i.e., crown surface area, crown volume and porosity) to the GNSS antenna error, artificial Gaussian noise was added to the GNSS positioning until reaching a maximum of 2 meters of absolute error (which is consistent with low-cost GNSS antennas), for each trial, and repeated ten times.Figure 5 shows part of the experimental data used in this work, where the LiDAR was positioned in several parts from the environment to register a point cloud.Figure 5a shows a partial view of a tree (green points) reconstructed with a LiDAR and GNSS antenna with 0.01 meters of absolute error in axis x and y.The same point cloud is then reconstructed using a GNSS antenna with a positioning error of 2 meters, as shown in Fig. 5d.As can be seen, as the absolute error increases, the point cloud becomes more disperse; in fact, several stems appear duplicated.

Propagation of GNSS error in crown parameter estimation
Figure 6-8 shows the results obtained during the trials.The three methodologies: voxelization, convex-hull and alpha-shape, were used to estimate the crown volume and crown surface area.However, to estimate porosity, only voxelization was used (Béland et al. 2014).These figures show the mean in solid magenta, the standard deviations in dashed dark line and single estimations in dotted dark line.
From Figure 6-8, a number of lessons are learned: • Regarding crown volume estimation, as the positioning error in the GNSS antenna increases up to 0.8 meters approximately, the estimated volume also increases for the three approaches (see Figs. 6a-6c).For the convex-hull and alpha-shape case, the volume estimation experiences a saturation behavior after one meter of error.Thus, if we increase the error in the positioning system, the volume estimation does not increase.6c shows the results for voxelization, convex-hull and alpha-shape, respectively.In all figures, the solid magenta line represents the mean estimation, whereas the dashed line is the standard deviation, and the dotted dark line are single estimations for each trial.
• When using voxelization, the estimated volume presented in Fig. 6a shows an opposite behavior compared to the other two approaches: it decays after 0.8 meters in the positioning error.This reason of this outcome is that as the error in the positioning increases, the cloud point becomes more disperse, as shown in Fig. 5. Therefore, the size of the voxels is smaller and the total estimated volume decreases.
• Regarding crown surface area estimation, it can be seen from Figs. 7a-  7c show the results for voxelization, convex-hull and alpha-shape, respectively.In all figures, the solid magenta line represents the mean estimation, whereas the dashed line is the standard deviation, and the dotted dark line are single estimations for each trial.
7c that the three approaches behave in a similar fashion: the estimation increases until approximately one meter of positioning error and then it remains almost constant.The voxelization approach, however, shows a slight decrease of approximately 10% with respect to the maximum observed.The latter is also associated with the fact that as the error in the position increases, the point cloud becomes more disperse and voxels smaller.• Regarding porosity, from the three methods, only voxelization allows for its estimation.As it can be seen, as the positioning error increases, porosity drastically decreases, which is consistent with porosity definition presented in (Béland et al. 2014;Pfeiffer et al. 2018): as the points become more dispersed, less shadow is expected.

GNSS and scan matching as a localization system
The estimated position obtained from the Kalman filter was used to transform the point cloud data to global coordinates.Once all point clouds were transformed, the orchard parameter estimation was obtained following the same methodology described in Section 3.1.Then, the results of both the original approach and the proposal are compared using the mean square error.

Orchard parameter estimation
The results showed in surface area is to increase.Based on the results obtained, the maximum is reached when GNNS error is equal to 2 meters, quite the contrary to the results from Figure 6-7 where the maximum is reached when GNSS error is equal to one meter.Therefore, the results obtained in this case give a more accurate estimation of orchard parameters: for crown surface area and crown volume estimation, for the three methodologies (voxelization, convex-hull and alpha-shape).
Additionally, Fig. 11 shows the results after evaluating porosity.As can be seen, when compared to Fig. 8, the decay of the porosity estimation as the GNSS error increases is smoother.The latter implies that the point cloud is less sparse when using the localization fusion approach to overcome the error in the GNSS antenna.chard parameters estimation.As can be seen, there is a reduction in the error associated with the estimation of all the features: crown surface area, crown volume, and porosity.The convex-hull results for crown surface area and volume are shown in Figure 12a and Figure 12b, respectively.In both cases, there is a notorious reduction of the error percentage.For example, in one meter of GNSS error, the results show that considering only GNSS as a localization system, the error increases up to 158 %, and when fused with the scan matching, the error increases up to 79%.Thus, an improvement of almost 50% was reached.The alpha-shape results show an improvement in the error associated with the estimation, as can be seen in Figures 12c     and 12d.The estimation errors decrease in at least 20% for all cases.The voxelization approach presents the results for crown surface area, crown volume and porosity in Figures 13a, 13b and 13c, respectively.From these results, it can be noted that the estimation error of the crown surface area and porosity has the same tendency as the previous methods.However, in the volume case, the presented approach seems to be worse than the GNSS trial after 1.2 meters of error.This outcome was obtained because the number of voxels within the workspace begins to decrease when exceeding 0.8 meters of GNSS error.Although such value is closer to the real value, it does not represent a consistent measure since the reconstructed point cloud is extremely dispersed.Even though our approach reduces the estimation error, it is subjected to similar effects to those of only the GNSS approach.13b show the voxelization approach mean square error for crown surface area and volume, respectively.Figure 13c shows the voxelization approach mean square error for porosity.

Conclusion
This work has analyzed the sensitivity of orchard features obtained via LiDAR, to GNSS positioning error.In particular, the crown volume, the crown surface area, and the porosity were studied, using three different computational approaches: convex-hull, alpha-shape, and voxelization.The experimentation was carried out at an apple grove.The analysis started with a positioning error of 0.01 meters provided by an RTK and then it was arti-ficially increased, adding Gaussian random noise, reaching up to two meters of error, as in low-cost GNSS receivers.The results considering only GNSS measurements have shown that crown surface area and crown volume experienced an increment as the positioning error grew up to one meter -for convex-hull and alpha-shape cases-, and then they remained almost stationary.For the voxelization case, crown surface area estimate behaved as in the previous approaches, but crown volume estimate decreased due to the reduction in the voxel size.In the case of porosity estimation, it decreased until positioning error reached one meter, and then started to grow again.To improve the above results, based only only on the LiDAR data available -thus avoiding the need of extra hardware-a scan matching approach (specifically, the ICP) was implemented and fused with the GNSS measurements, in order to decrease the localization error.The results showed a notorious improvement against the original implementation.Errors in the crown surface area, crown volume and porosity estimation reduced up to 20%, by adding the ICP when the GNSS error was of 1.2 meters and up to 50% for smaller errors.These outcomes suggest that the proposed approach can improve the results of foliar parameters estimation, therefore, leading to better decisions in agricultural operations without an increment in costs of external hardware.

Figure 1 :
Figure 1: Hardware used in this work.The platform is equipped with a Velodyne PUCK LiDAR and an RTK Leica GPS1200.

Figure 3 :
Figure 3: Methodology used for improving the foliar paramater estimation results.The previous and the current laser scans enter in a scan matching algorithm and are fused with the GNSS data.Once all laser scans are transformed to global coordinate frame considering the output of the Kalman Filter, the crown foliar parameters are estimated.

Figure 4 :
Figure 4: The transformation matrix H Ω allows to align the point cloud Ω k with the point cloud Λ k .

Figure 5 :
Figure 5: The ground-truth and three point clouds with different positioning error associated with the GNSS antenna.The GNSS error was established in 0.01, 0.2, 1 and 2 meters for Figure 5a-Fig.5d, respectively.

Figure 6 :
Figure6: Experimental results for crown volume estimation.Figures6a-6cshows the results for voxelization, convex-hull and alpha-shape, respectively.In all figures, the solid magenta line represents the mean estimation, whereas the dashed line is the standard deviation, and the dotted dark line are single estimations for each trial.

Figure 7 :
Figure7: Experimental results for crown surface area estimation.Figures7a-7cshow the results for voxelization, convex-hull and alpha-shape, respectively.In all figures, the solid magenta line represents the mean estimation, whereas the dashed line is the standard deviation, and the dotted dark line are single estimations for each trial.

Figure 8 :
Figure 8: Porosity results using voxelization.The solid magenta line represents the mean estimation, whereas the dashed line is the standard deviation, and the dotted dark line are single estimations for each trial.

Figure 9 :
Figure 9: Experimental results.Figures 9a-9c show the crown volume estimated with the three different methods implemented in this work: voxelization, convex-hull and alphashape, respectively.

Figure 10 :
Figure 10: Figures 10a-10c show the results for crown surface area estimation.The solid magenta line represents the mean estimation, whereas the dashed line is the standard deviation, and the dotted dark line are single estimations for each trial.

Figure 12 -Figure 11 :
Figure12-13 shows a comparison of the mean square error between the results obtained with GNSS only and its fusion with scan matching for or-

Figure 12 :Figure 13 :
Figure12: Mean square error percentage for the original (GNSS only) and the proposed (GNSS with scan matching) approach.The convex-hull results for crown surface area Figures 12a and volume 12b; Alpha-shape surface area and volume in Figures12c and 12d.The blue and the red bar represent the error percentage in the volume estimation process when only GNSS is consider and when GNSS is fused with scan matching measurements, respectively.