Estimating Canopy Parameters Based on the Stem Position in Apple Trees Using a 2D LiDAR
利用2D LiDART根据苹果树茎位置估计冠层参数

Abstract  摘要

Data of canopy morphology are crucial for cultivation tasks within orchards. In this study, a 2D light detection and range (LiDAR) laser scanner system was mounted on a tractor, tested on a box with known dimensions (1.81 m × 0.6 m × 0.6 m), and applied in an apple orchard to obtain the 3D structural parameters of the trees (n = 224). The analysis of a metal box which considered the height of four sides resulted in a mean absolute error (MAE) of 8.18 mm with a bias (MBE) of 2.75 mm, representing a root mean square error (RMSE) of 1.63% due to gaps in the point cloud and increased incident angle with enhanced distance between laser aperture and the object. A methodology based on a bivariate point density histogram is proposed to estimate the stem position of each tree. The cylindrical boundary was projected around the estimated stem positions to segment each individual tree. Subsequently, height, stem diameter, and volume of the segmented tree point clouds were estimated and compared with manual measurements. The estimated stem position of each tree was defined using a real time kinematic global navigation satellite system, (RTK-GNSS) resulting in an MAE and MBE of 33.7 mm and 36.5 mm, respectively. The coefficient of determination (R²) considering manual measurements and estimated data from the segmented point clouds appeared high with, respectively, R² and RMSE of 0.87 and 5.71% for height, 0.88 and 2.23% for stem diameter, as well as 0.77 and 4.64% for canopy volume. Since a certain error for the height and volume measured manually can be assumed, the LiDAR approach provides an alternative to manual readings with the advantage of getting tree individual data of the entire orchard.
冠层形态数据是果园栽培任务的关键。在这项研究中，一个二维光探测和距离(LiDAR)激光扫描系统安装在拖拉机上，在一个已知尺寸(1.81m×0.6m×0.6m)的盒子上进行测试，并应用于一个苹果园(n=224)，获得树木的三维结构参数。对考虑四边高度的金属盒进行分析，得到的平均绝对误差(MAE)为8.18 mm，偏差(MBE)为2.75 mm，均方根误差(RMSE)为1.63%，这是由于点云中的缝隙和入射角随着激光孔径与物体距离的增加而增加。提出了一种基于二维点密度直方图的方法来估计每棵树的树干位置。圆柱形边界被投影在估计的树干位置周围，以分割每一棵树。随后，对分割的树木点云的高度、树干直径和体积进行了估计，并与人工测量进行了比较。使用实时动态全球导航卫星系统(RTK-GNSS)定义了每棵树的估计树干位置，MAE和MBE分别为33.7 mm和36.5 mm。考虑人工测量和分段点云估计数据的决定系数(R ² )较高，高度的R ²

1. Introduction  1.引言

Information on geometrical and structural characteristics of fruit trees provides decisive knowledge for management within the orchard. The geometrical (height, width, and volume) and the structural (the leaf area, leaf density, trunk, and angle of branches) variables affect economically relevant parameters, such as the tree water needs [1,2,3], fruit quality [4,5], and the yield [6]. The manual measurement of tree variables is usually inaccurate and time-consuming. Thus, the performance of measurements in orchards for management purposes lack feasibility at the present state of the art.
果树的几何和结构特征信息为果园的管理提供了决定性的知识。几何变量(高度、宽度和体积)和结构变量(叶面积、叶密度、树干和树枝角度)影响相关的经济参数，如树木需水[1、2、3]、果实质量[4、5]和产量[6]。人工测量树木变量通常是不准确和耗时的。因此，出于管理目的在果园中进行测量在当前技术状态下缺乏可行性。

Among other instrumental methods for this purpose, the light detection and range (LiDAR) laser scanning seems to be the most adequate for the nondestructive monitoring of geometrical and structural parameters in fruit trees [7], because it (i) operates with an active light source that avoids perturbing effects of varying lighting conditions, (ii) enables higher frequencies than those provided by manual measurements, and (iii) enables us to measure from a distance in multiple directions using one single sensor. A laser diode emits photons with low dispersion, which are reflected back to the photodiode of the scanner when hiting an object. The time between the emission and the reflection of the beam is proportional to the distance among the scanner and the objects, considering the constant of light velocity and sensor response functions. A motor sweeps the laser beam in a plane, describing the shape of objects in two dimensions (2D).
在其他用于此目的的仪器方法中，光检测和距离(LiDAR)激光扫描似乎是最适合于无损监测果树几何和结构参数的方法[7]，因为它(I)使用主动光源，避免了不同光照条件的干扰影响，(Ii)实现了比人工测量更高的频率，并且(Iii)使我们能够使用单个传感器从一个距离的多个方向进行测量。激光二极管发射低色散的光子，当击中物体时，这些光子被反射回扫描仪的光电二极管。在考虑光速不变和传感器响应函数不变的情况下，光束从发射到反射的时间与扫描器与物体之间的距离成正比。马达在平面上扫描激光束，描述物体在二维(2D)上的形状。

However, in agriculture, the 2D LiDAR has been frequently mounted on a tractor together with a high precision geographical position system to obtain three dimensional (3D) data of trees [8,9,10]. Several studies used this configuration to estimate the structural parameters of fruit trees by segmenting the row point cloud into slices of a constant distance [4,9,10,11]. The canopy volume of the tree row revealed a high correlation coefficient (r) with the manual reference measurement (r = 0.906, p < 0.0001) in olive orchards using the aforementioned method [10]. In another study, Sanz and co-workers [12] used the volume of a tree row to obtain the leaf area, describing a high coefficient of determination (R²) in different chronological stages with a logarithmic model in apple trees (R² = 0.85), pear trees (R² = 0.84), and vines (R² = 0.86). Moreover, several methodologies have been developed for volume extraction such as the convex-hull clustering or voxel grid in vines, olives, and apples [2,13,14], while the alpha shape clustering was applied in citrus [15]. Jimenez-Berni and co-authors [16] described the analysis of plant height in wheat with a strong correlation (R² = 0.99), when using a 2D LiDAR. Furthermore, the leaf area index can be estimated by utilizing the cross-sectional area in vines and apple trees [5,13,17]. Similarly, correlations of pruning weight and vine wood volume resulted in high coefficients of determination with LiDAR measurements (R² = 0.91; R² = 0.76, respectively) [18].
然而，在农业中，2D LiDAR经常与高精度地理定位系统一起安装在拖拉机上，以获得树木的三维(3D)数据[8，9，10]。一些研究使用这种构型通过将行点云分割成恒定距离的切片来估计果树的结构参数[4，9，10，11]。用上述方法测得的橄榄园树行冠层体积与人工参考值具有较高的相关系数(r=0.906，p<0.0001)。在另一项研究中，Sanz和他的同事[12]使用树行的体积来获得叶面积，描述了不同时间阶段的高决定系数(R ² )，在苹果树(R ² =0.85)、梨树(R ² =0.84)和藤本植物(R ² =0.86)中采用对数模型。此外，已经开发了几种方法来提取体积，例如在藤本植物、橄榄和苹果上应用凸壳聚类或体素网格[2，13，14]，而在柑橘上应用阿尔法形状聚类[15]。Jimenez-Berni和合著者[16]描述了当使用2D LiDAR时，小麦株高的分析具有很强的相关性(R ² =0.99)。此外，叶面积指数可以通过利用葡萄藤和苹果树的横截面积来估计[5，13，17]。同样，修剪重量和藤材体积的相关性导致了与LiDAR测量的高决定系数(分别为R ² =0.91；R ² =0.76)[18]。

The point cloud clustering for separating single plants was already performed with fixed distances between the trees [13], the convex-hull approach [14], and hidden semi-Markov model [19,20]. Reiser et al. [20] determined the stem position of maize plants using an initial plant and the fixed plant spacing to search iteratively for the best clusters. The height profile of the clustered maize plants estimated, revealed a 33 mm deviation from the manual reference readings. Garrido et al. [21] used a LiDAR and a light curtain to detect stem positions of almond trees, with a detection accuracy of 99.5%. The position of the detected stems was obtained with the help of an optical wheel odometer that separated plants with no overlapping leaves. Zhang and Grift [22] proposed a LiDAR height measurement system for Miscanthus giganteus based on stem density. In static mode, the sensor kept stationary and an average error of 5.08% was achieved, while in the dynamic mode, the sensor was driven past a field edge, which resulted in an average error of 3.8%. In another study, the 2D histogram was used to obtain the height profile in 3D point clouds of cotton plants, resulting in a strong correlation (R² = 0.98) with manual measurements [23]. Moreover, the stem position in maize plants was determined based on bivariate point density histograms in order to find the regional maxima [24]. Thus, single plant segmentation was carried out by projecting a spatial cylindrical boundary around the estimated stem positions of plant point clouds. The maize stem positions were estimated with mean error and standard deviation of 24 mm and 14 mm difference from the actual stem position, respectively. For the overall plant height profile, an error of 8.7 mm was found while considering the manually obtained height.
用于分离单个植物的点云聚类已经使用树之间的固定距离[13]、凸壳方法[14]和隐式半马尔可夫模型[19，20]来执行。Reiser等人。[20]利用初始植株和固定的植株间距来确定玉米植株的茎位置，以迭代地搜索最佳聚类群。丛生玉米植株的高度轮廓估计，显示出与手动参考读数的33毫米偏差。Garrido等人。[21]利用激光雷达和光幕检测杏树主茎位置，检测准确率达99.5%。检测到的茎的位置是在光学轮式里程计的帮助下获得的，该里程计将没有重叠叶片的植物分开。张和格里夫特[22]提出了一种基于茎密度的大芒激光雷达高度测量系统。在静态模式下，传感器处于静止状态，平均误差为5.08%；在动态模式下，传感器被驱动通过磁场边缘，平均误差为3.8%。在另一项研究中，使用2D直方图来获得棉花植株3D点云中的高度轮廓，结果与人工测量结果有很强的相关性(R ² =0.98)[23]。此外，基于双变量点密度直方图确定了玉米植株的茎位置，以求出区域极大值[24]。因此，单个植物的分割是通过在估计的植物点云的茎位置周围投影一个空间圆柱边界来实现的。估计的玉米茎位与实际茎位的平均误差和标准差分别为24 mm和14 mm。对于整个植株高度轮廓，误差为8。考虑到手动获得的高度，发现7 mm。

The overall aim of this research was to propose a methodology for determining the stem position of apple trees by means of a 2D LiDAR system as a basis for separation of trees in commercial orchards. The specific objectives were to test the measuring uncertainty of the LiDAR system, when analyzing a box with known dimensions in the field and real-world trees. The analysis of tree point clouds was evaluated by estimating the height, the stem and the canopy volume of singularized trees. The innovation of this study is the application of the proposed methodology in fruit trees, providing an accurate and precise estimation of geometric variables of apple trees within the orchard.
这项研究的总体目标是提出一种通过2D LiDART系统确定苹果树茎位置的方法，作为商业果园中树木分离的基础。具体目标是在分析现场已知尺寸的盒子和现实世界的树木时测试LiDART系统的测量不确定性。通过估计奇异树木的高度、树干和树冠体积来评估树木点云的分析。本研究的创新之处在于将所提出的方法学应用于水果树，为果园内苹果树的几何变量提供了准确且精确的估计。

2. Material and Methods
2.材料和方法

2.1. Field Experiments  2.1.田间试验

The experiment was conducted in the Leibniz Institute of Agricultural Engineering and Bioeconomy (ATB) experimental station located in Marquardt, Germany (Latitude: 52.466274° N, Longitude: 12.57291° E) in early July of 2018. The orchard was located on 8% slope with southeast orientation. The orchard was planted with trees of Malus × domestica Borkh. ‘Gala’ and ‘JonaPrince’, and pollinator trees ‘Red Idared’ each on M9 rootstock (n = 224) with 0.95 m distance between trees, trained as slender spindle, which form the majority of apple trees in world-wide production, with an average tree height of 2.5 m. No pollinator trees were excluded from the study, since they appeared to have a similar tree architecture to the aforementioned tree types, and were also trained as slender spindles. Trees were supported by means of horizontally parallel wires. Data acquisition was performed at the end of fruit cell division, after a native fruit drop in June. The tractor that was carrying the sensor-frame was driven along each side of the two rows measured, with an average speed of 0.13 m s⁻¹.
该实验于2018年7月初在位于德国马夸特(北纬52.466274°，东经12.57291°)的莱布尼茨农业工程与生物经济研究所实验站进行。果园位于8%的坡度上，朝东南方向。果园内栽种了苹果树。‘Gala’和‘JonaPrince’和传粉树‘Red Idred’各在M9砧木(n=224)上，树间距0.95m，训练成细长纺锤，构成了世界范围内生产的大多数苹果树，平均树高2.5m。没有传粉树被排除在研究之外，因为它们似乎具有类似于上述树型的树体结构，也被训练成细长纺锤形。树是用水平平行的钢丝支撑的。数据采集是在果细胞分裂结束时进行的，这是在6月份本地果实脱落后进行的。承载传感器框架的拖拉机沿被测两排的两侧行驶，平均车速为0.13m S ⁻¹ 。

The stem position of each tree (n = 224) was defined with a real time kinematic global navigation satellite system (RTK-GNSS), while in parallel a wooden meter was used to estimate the tree height manually (H_manual). The average width of each tree (W_manual) was specified at three different heights across the x and y axis. Furthermore the stem diameter (S_manual) was measured at 0.3 m height, which was always above the grafting area and above the ground, using a tape scale. The value of the W_manual in x and y axes and H_manual of the tree was multiplied to determine the volume of each tree (V_manual) [4].
每棵树（n = 224）的树干位置使用实时运动学全球导航卫星系统（RTK-GNSS）定义，同时使用木表手动估计树高（H）。每棵树的平均宽度（W）在x和y轴的三个不同高度指定。此外，使用卷尺在0.3 m高度（始终在嫁接区域上方和地面上方）测量茎直径（S）。将树的x和y轴上的W值与H相乘，以确定每棵树的体积（V）[ 4]。

2.2. Instrumentation  2.2.仪表

A rigid aluminum sensor-frame carrying the sensors (Figure 1) was mounted on the front three-point-hitch of a tractor. A mobile terrestrial light detection and ranging (LiDAR) (LMS-511, Sick AG, Waldkirch, Germany) laser scanner was mounted at 1.6 m above the ground level. Τwo digital filters were activated to optimize the measured distance values: a fog filter that could not be avoided according to the manufacturer’s setting; and a N-pulse-to-1-pulse filter, which filtered out the first reflected pulse in cases where two pulses were reflected by two objects during a measurement (Table 1).
承载传感器的刚性铝制传感器框架（图1）安装在拖拉机的前部三点式挂钩上。移动式地面光检测和测量（LiDART）（MPS-511，Sick AG，Waldkirch，德国）激光扫描仪安装在地平面上方1.6 m处。激活了两个数字过滤器以优化测量的距离值：根据制造商的设置无法避免的雾过滤器;和N脉冲到1脉冲过滤器，在测量期间两个脉冲被两个物体反射的情况下过滤掉第一个反射脉冲（表1）。

Figure 1. (a) Representation of the sensor-frame system showing the coordinate system of LiDAR and a real time kinematic global navigation satellite system (RTK-GNSS). (b) Metal box of known distances.
图1. (a)传感器框架系统的表示，显示LiDART和实时动态全球导航卫星系统（RTK-GNSS）的坐标系。(b)已知距离的金属盒。

Table 1. LMS 511 specification data [25].

Functional Data	General Data
Operating range: up to 10 m	Light detection and ranging (LiDAR) Class: 1 (IEC 60825-1)
Scanning angle: 190°	Enclosure rating: IP 67
Scanning frequency: 25 Hz	Temperature range: −40 °C to 60 °C
Systematic error: ± 25 mm	Light source: 905 nm (near infrared)
Statistical error: ± 6 mm	Total weight: 3.7 kg
Angular resolution: 0.1667°	Laser beam diameter at the front screen: 13.6 mm

The three-dimensional tilt of the tractor was monitored by an inertial measurement unit (IMU) MTi-G-710 (XSENS, Enschede, Netherlands), which was placed 0.24 m aside from the LiDAR. The accuracy of the sensor was 1.0° root mean square error (RMSE) for the heading at static and dynamic mode, and 0.25° RMSE and 1.0° RMSE for both pitch and roll in static and dynamic mode, respectively [26]. The measured data was georeferenced by an AgGPS 542 RTK-GNSS (Trimble, Sunnyvale, CA, USA) mounted at 1.74 m height. The horizontal accuracy of the RTK-GNSS was ±25 mm +2 ppm and the vertical ±37 mm +2 ppm.
拖拉机的三维倾斜由惯性测量单元（IMU）MTi-G-710（XSENS，Enschede，荷兰）监测，该单元放置在距离LiDART 0.24 m处。静态和动态模式下航向传感器的准确度分别为1.0° RSSE，静态和动态模式下俯仰和横滚的准确度分别为0.25° RSSE和1.0° RSSE [ 26]。测量数据通过安装在1.74 m高度的AgGPS 542 RTK-GNSS（Trimble，Sunnyvale，CA，USA）进行地理参考。RTK-GNSS的水平精度为±25 mm +2 ppm，垂直精度为±37 mm +2 ppm。

Software with multi-thread architecture was developed in Visual Studio (version 16.1, Microsoft, Redmond, WA, US) in order to acquire the sensor data. The parallel data acquisition was achieved by creating a dynamic thread for every connected sensor. Consequently, when the raw sensor data were read by each thread, the corresponding measurement and the time stamp from real-time clock (1 ms resolution) were assembled into a text string. The string was stored in a text file for post-processing. The multi-threaded data acquisition of the sensors resulted in non-concurrent measurements. Therefore, the cubic spline interpolation method was used before synchronizing time instances, based on their individual timestamps.
具有多线程架构的软件是在Visual Studio（版本16.1，Microsoft，Redmond，WA，US）中开发的，以获取传感器数据。通过为每个连接的传感器创建动态线程来实现并行数据采集。因此，当每个线程读取原始传感器数据时，相应的测量结果和实时时钟（1 ms分辨率）的时间戳被组装到文本字符串中。该字符串存储在文本文件中以进行后处理。传感器的多线程数据采集导致非并发测量。因此，在同步时间实例之前，根据其单独的时间戳使用三次样条插值方法。

2.3. Sensor-Frame Analysis
2.3.传感器框架分析

A metal box with dimensions 1.81 m × 0.60 m × 0.60 m was constructed to validate the measuring uncertainty of the sensor-frame system (Figure 1b). Five bars with 0.3 m distance from each other were placed horizontally on the metal box. The sensor-frame mounted on a tractor was applied to scan each side of the box separately with 0.13 m s⁻¹. The measurements were carried out next to the apple trees in order to simulate the real-world operating conditions. The mean absolute error (MAE), the bias (MBE) and the RMSE considering the actual and estimated dimensions of the metal box were calculated (Equations (1)–(3))
建造了一个尺寸为1.81m×0.60m×0.60m的金属盒子，以验证传感器-框架系统的测量不确定度(图1b)。金属盒子上水平放置了五根相距0.3m的棒子。安装在拖拉机上的传感器框架以0.13m的S ⁻¹ 分别扫描盒子的两侧。这些测量是在苹果树旁边进行的，以模拟真实的操作条件。计算了考虑金属盒实际尺寸和估计尺寸的平均绝对误差(MAE)、偏差(MBE)和均方根误差(方程式(1)-(3))。

MBE (mm) = ∑ni=1(Dr,i−DLiDAR, i)n$$\mathrm{MBE}(\mathrm{mm})=\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{D}}_{\mathrm{r},\mathrm{i}}-{\mathrm{D}}_{\mathrm{LiDAR},\mathrm{i}}\right)}{\mathrm{n}}$$

(1)

MAE (mm) = ∑ni=1|Dr,i−DLiDAR, i|n$$\mathrm{MAE}(\mathrm{mm})=\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\left|{\mathrm{D}}_{\mathrm{r},\mathrm{i}}-{\mathrm{D}}_{\mathrm{LiDAR},\mathrm{i}}\right|}{\mathrm{n}}$$

(2)

RMSE (%) = ∑ni=1(Dr,i−DLiDAR,i)n−−−−−−−−−−−−−−−−−√$$\mathrm{RMSE}(\%)=\sqrt{\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{D}}_{\mathrm{r},\mathrm{i}}{-\mathrm{D}}_{\mathrm{LiDAR},\mathrm{i}}\right)}{\mathrm{n}}}$$

(3)

where D_r,i is the actual reference dimension of the metal box profile, D_LiDAR,i is the dimension of the metal box estimated by the LiDAR for time instance i, and n is the number of the measured points for each profile. Furthermore, the above equations were used to estimate the measuring uncertainty considering the manual and the LiDAR measurements of the stem position, tree height, stem diameter, and canopy volume.
其中，D^b0b0是金属盒轮廓的实际参考尺寸，D^b11是由LiDAR对时间实例i估计的金属盒的尺寸，n是每个轮廓的测量点的数目。此外，考虑到树干位置、树高、树干直径和冠层体积的手动测量和激光雷达测量，上述方程被用来估计测量不确定度。

2.4. Data Pre-Processing  2.4.数据预处理

The field data were processed using mainly the Computer Vision System Toolbox^TM of Matlab (2017b, MathWorks, Natick, MA, USA), the function of the statistical outlier filter (SOR), and cloud to cloud distance (C2C) of CloudCompare (EDF R&D) for point cloud processing.
现场数据的处理主要使用MatLab(2007b，MathWorks，Natick，MA，USA)的计算机视觉系统工具箱 ^TM 、统计离群滤波功能(SOR)和点云处理的CloudCompare(EDF R&D)的云到云距离(C2C)。

The LiDAR data set was filtered for considering only those detections with distance measurements between 0.05 m and 4.00 m. The resulting detections were transformed from polar to Cartesian coordinates (x_LiDAR, y_LiDAR, z_LiDAR). The 3D dynamic tilt orientation (yaw (ψ), pitch (θ), and roll (φ)) of the sensor-frame during the tractor movement was acquired by the IMU sensor in east north up (ENU) local frame. Thus, any point of interest on the sensor-frame could be calculated using coordinate transformation. The position data set, latitude (φ), longitude (λ), and altitude (h) (LLA) obtained by the RTK-GNSS receiver, was expressed in the geodetic frame defined by WGS84.
对LiDAR数据集进行过滤，只考虑距离测量在0.05m和4.00m之间的那些探测。结果探测从极坐标转换到笛卡尔坐标(x^{b0>，y _LiDAR ，z _LiDAR )。在拖拉机运动过程中，传感器框架的三维动态倾斜方向(偏航(ψ)、俯仰(θ)和横摇(φ))是由东北上(ENU)局部坐标系中的惯性测量单元传感器获取的。因此，可以使用坐标变换来计算传感器框架上的任何兴趣点。RTK-GNSS接收机获得的位置数据集，即纬度(φ)、经度(λ)和高度(H)(Lla)，以WGS84定义的大地坐标系表示。}

The RTK-GNSS and the IMU output data had to be converted in the same coordinate system in order to be fused. The fusion was carried out in ENU, since the IMU output data was already in this local coordinate system. However, the RTK-GNSS output data was converted twice, firstly the LLA to earth-centered earth fixed (ECEF) frame and secondly to the ENU frame.
RTK-GNSS和IMU的输出数据必须在同一坐标系中进行转换，以便进行融合。融合是在ENU进行的，因为IMU的输出数据已经在这个地方坐标系中。然而，RTK-GNSS的输出数据被转换了两次，第一次是LLA到地球中心固定(ECEF)框架，第二次是ENU框架。

The conversion from the geodetic coordinates to the ECEF coordinates was performed using the following equations [27]
使用以下方程执行从大地坐标到ECEF坐标的转换[ 27]

xe = (ν(φ) + h)cos(φ)cos(λ),$${\mathrm{x}}_{\mathrm{e}}=\left(\mathsf{\nu }\right(\mathsf{\phi })+\mathrm{h})\cos \left(\mathsf{\phi }\right)\cos \left(\mathsf{\lambda }\right),$$

(4)

ye = (ν(φ) + h)cos(φ)sin(λ),$${\mathrm{y}}_{\mathrm{e}}=\left(\mathsf{\nu }\right(\mathsf{\phi })+\mathrm{h})\cos \left(\mathsf{\phi }\right)\sin \left(\mathsf{\lambda }\right),$$

(5)

ze = [ν(φ)(1 − e2)+h]sin(φ),$${\mathrm{z}}_{\mathrm{e}}{=\left[\mathsf{\nu }\right(\mathsf{\phi }\left)\right(1-\mathrm{e}}^2)+\mathrm{h}]\sin \left(\mathsf{\phi }\right),$$

(6)

where e = 0.0818 is the first numerical eccentricity of the earth ellipsoid, Π_e = [xe,ye ,ze${\mathrm{x}}_{\mathrm{e}}{,\mathrm{y}}_{\mathrm{e}}{,\mathrm{z}}_{\mathrm{e}}$]^T is the position in ECEF coordinates and ν(φ) is the distance from earth’s surface to the z-axis along the ellipsoid normal given by
其中e = 0.0818是地球椭圆体的第一个数字同心度，ð=[]是ECEF坐标中的位置，v（ph）是地球表面沿着椭圆体法向z轴的距离

ν(φ) = α1−e2sin2φ−−−−−−−−−√$$\mathsf{\nu }(\mathsf{\phi })=\frac{\mathsf{\alpha }}{\sqrt{1-{\mathrm{e}}^2{\sin }^2\mathsf{\phi }}}$$

(7)

where α = 6,378,137 m is the semi-major axis length of the earth. The conversion from the ECEF coordinates to the local ENU coordinates was estimated from the equation Π_n = Tne${\mathrm{T}}_{\mathrm{e}}^{\mathrm{n}}$Π_e, where Π_n = [n, e, d]^T is the position in the local ENU coordinate system, and Tne${\mathrm{T}}_{\mathrm{e}}^{\mathrm{n}}$ is the transformation matrix from the ECEF to the ENU coordinate system
其中，a = 6，378，137 m是地球的半长轴长度。从ECEF坐标到本地ENU坐标的转换是根据方程RST = RST估计的，其中RST = [n，e，d]是本地ENU坐标系中的位置，并且是从ECEF到ENU坐标系的转换矩阵

Tne = ⎡⎣⎢⎢⎢⎢⎢−sin(λ)−sin(φ)cos(λ)cos(φ)cos(λ)0cos(λ)−sin(φ)sin(λ)sin(λ)cos(φ)00cos(φ)sin(φ)00001⎤⎦⎥⎥⎥⎥⎥.$${\mathrm{T}}_{\mathrm{e}}^{\mathrm{n}}=\left[\begin{array}{cccc}-\sin \left(\mathsf{\lambda }\right)& \cos \left(\mathsf{\lambda }\right)& 0& 0\\ -\sin \left(\mathsf{\phi }\right)\cos \left(\mathsf{\lambda }\right)& -\sin \left(\mathsf{\phi }\right)\sin \left(\mathsf{\lambda }\right)& \cos \left(\mathsf{\phi }\right)& 0\\ \cos \left(\mathsf{\phi }\right)\cos \left(\mathsf{\lambda }\right)& \sin \left(\mathsf{\lambda }\right)\cos \left(\mathsf{\phi }\right)& \sin \left(\mathsf{\phi }\right)& 0\\ 0& 0& 0& 1\end{array}\right].$$

(8)

The sensor-frame enabled transformations and translation due to the known distances and angles between the sensors. The LiDAR frame was transformed into the RTK-GNSS frame by several translations (Trans (yaxis ${\mathrm{y}}_{\mathrm{axis}}$, −a) × Trans (zaxis ${\mathrm{z}}_{\mathrm{axis}}$, +b)) and two rotations (Rot (‘roll’,−π2$\frac{\mathsf{\pi }}{2}$) × Rot (‘yaw’, −𝜋2$\frac{\pi }{2}$) in order to georeference each point of the LiDAR
由于传感器之间的已知距离和角度，传感器框架能够实现转换和平移。通过多次平移（Trans（，-a）x Trans（，+b））和两次旋转（Rot（' roll '，-）x Rot（' yaw '，-）将LiDART框架转换为RTK-GNSS框架，以便对LiDART的每个点进行地理参考

[xsensor−frame, ysensor−frame, zsensor−frame, 1]T = [Tsensor−frameLiDAR] × [xLiDAR, yLiDAR, zLiDAR, 1]T$${[{\mathrm{x}}_{\mathrm{sensor}-\mathrm{frame}},{\mathrm{y}}_{\mathrm{sensor}-\mathrm{frame}},{\mathrm{z}}_{\mathrm{sensor}-\mathrm{frame}}, 1]}^{\mathrm{T}}=[{\mathrm{T}}_{\mathrm{LiDAR}}^{\mathrm{sensor}-\mathrm{frame}}]\times [{\mathrm{x}}_{\mathrm{LiDAR}},{\mathrm{y}}_{\mathrm{LiDAR}},{\mathrm{z}}_{\mathrm{LiDAR}}, 1{]}^{\mathrm{T}}$$

(9)

where Tsensor−frameLiDAR${\mathrm{T}}_{\mathrm{LiDAR}}^{\mathrm{sensor}-\mathrm{frame}}$ is the transformation matrix from the coordinate system of the LiDAR sensor to the coordinate system of the RTK-GNSS
其中是LiDART传感器坐标系到RTK-GNSS坐标系的变换矩阵

Tsensor−frameLiDAR = Rot (‘roll’, −π2) × Rot (‘yaw’, −𝜋2)× Trans (yaxis, −a) × Trans (zaxis, +b).$${\mathrm{T}}_{\mathrm{LiDAR}}^{\mathrm{sensor}-\mathrm{frame}}=\mathrm{Rot}(‘\mathrm{roll}’,-\frac{\mathsf{\pi }}{2})\times \mathrm{Rot}(‘\mathrm{yaw}’,-\frac{\pi }{2})\times \mathrm{Trans}({\mathrm{y}}_{\mathrm{axis}},-\mathrm{a})\times \mathrm{Trans}({\mathrm{z}}_{\mathrm{axis}},+\mathrm{b}).$$

(10)

The general matrix form of all rotations (Rot) and translations (Trans) are shown in Equation (11), and by substituting them in Equation (10), a matrix form representation of the equations can be generated. The three-dimensional tilt orientations were implemented in the point cloud considering the LiDAR timestamp (t).
所有旋转（Rot）和平移（Trans）的一般矩阵形式如方程（11）所示，并且通过将它们替换到方程（10）中，可以生成方程的矩阵形式表示。考虑LiDART时间戳，在点云中实现了三维倾斜方向（）。

Rot (‘roll’, φ) = ⎡⎣⎢⎢⎢⎢⎢100 00cos(φt)sin(φt)00−sin(φt)cos(φt)00001⎤⎦⎥⎥⎥⎥⎥, Trans (‘xaxis’, xtranslation) = ⎡⎣⎢⎢⎢⎢10 0 xtranslation010000100001⎤⎦⎥⎥⎥⎥,Rot (‘pitch’, θ) = ⎡⎣⎢⎢⎢⎢⎢cos(θt)0−sin(θt) 00100sin(θt)0cos(θt)00001⎤⎦⎥⎥⎥⎥⎥, Trans (‘yaxis’, ytranslation) = ⎡⎣⎢⎢⎢⎢100 0010𝑦𝑡𝑟𝑎𝑛𝑠𝑙𝑎𝑡𝑖𝑜𝑛00100001⎤⎦⎥⎥⎥⎥,Rot (‘yaw’, ψ) = ⎡⎣⎢⎢⎢⎢⎢cos(𝜓𝑡)sin(𝜓𝑡)0 0−sin(𝜓𝑡)cos(𝜓𝑡)0000100001⎤⎦⎥⎥⎥⎥⎥, Trans (‘zaxis’, ztranslation) = ⎡⎣⎢⎢⎢⎢100 00100001ztranslation0001⎤⎦⎥⎥⎥⎥$$\begin{array}{c}\mathrm{Rot}(‘\mathrm{roll}’,\mathsf{\phi })=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& {\cos (\mathsf{\phi }}_{\mathrm{t}})& {-\sin (\mathsf{\phi }}_{\mathrm{t}})& 0\\ 0& {\sin (\mathsf{\phi }}_{\mathrm{t}})& {\cos (\mathsf{\phi }}_{\mathrm{t}})& 0\\ 0& 0& 0& 1\end{array}\right],\mathrm{Trans}(‘{\mathrm{x}}_{\mathrm{axis}}’,{\mathrm{x}}_{\mathrm{translation}})=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ {\mathrm{x}}_{\mathrm{translation}}& 0& 0& 1\end{array}\right],\\ \mathrm{Rot}(‘\mathrm{pitch}’,\mathsf{\theta })=\left[\begin{array}{cccc}{\cos (\mathsf{\theta }}_{\mathrm{t}})& 0& {\sin (\mathsf{\theta }}_{\mathrm{t}})& 0\\ 0& 1& 0& 0\\ {-\sin (\mathsf{\theta }}_{\mathrm{t}})& 0& {\cos (\mathsf{\theta }}_{\mathrm{t}})& 0\\ 0& 0& 0& 1\end{array}\right],\mathrm{Trans}(‘{\mathrm{y}}_{\mathrm{axis}}’,{\mathrm{y}}_{\mathrm{translation}})=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& y_{translation}& 0& 1\end{array}\right],\\ \mathrm{Rot}(‘\mathrm{yaw}’,\mathsf{\psi })=\left[\begin{array}{cccc}\cos ({\psi }_t)& -\sin ({\psi }_t)& 0& 0\\ \sin ({\psi }_t)& \cos ({\psi }_t)& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right],\mathrm{Trans}(‘{\mathrm{z}}_{\mathrm{axis}}’,{\mathrm{z}}_{\mathrm{translation}})=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& {\mathrm{z}}_{\mathrm{translation}}& 1\end{array}\right]\end{array}$$

(11)

The three-dimensional tilt of the IMU data was implemented to stabilize the orientation of the coordinate system of the sensor-frame (Equation (12))
实施IMU数据的三维倾斜以稳定传感器框架坐标系的方向（方程（12））

Tts−geoLiDAR = Tsensor−frameLiDAR× Rot (‘roll’,−φ) ×Rot (‘pitch’,−θ) × Rot (‘yaw’,ψ).$${\mathrm{T}}_{\mathrm{LiDAR}}^{\mathrm{ts}-\mathrm{geo}}={\mathrm{T}}_{\mathrm{LiDAR}}^{\mathrm{sensor}-\mathrm{frame}}\times \mathrm{Rot}(‘\mathrm{roll}’,-\mathsf{\phi })\times \mathrm{Rot}(‘\mathrm{pitch}’,-\mathsf{\theta })\times \mathrm{Rot}(‘\mathrm{yaw}’,\mathsf{\psi }).$$

(12)

2.5. Point Cloud Rigid Registration and Stitching
2.5.点云刚性配准和缝合

The process started by implementing the SOR filter to each point cloud pair to first compute the average distance of each point to its neighbors and then reject the points that are further than the average distance, while also considering the standard deviation. Furthermore, the random sample consensus (RANSAC) algorithm was used to remove points that belonged to the ground. The distance threshold value between a data point and the plane was defined at 5 mm, while the slope of the field was additionally considered in the plane model through the implementation of the RTK-GNSS height difference of each row. As mentioned above, a plane with 5 mm height was applied in both rows in order to remove the reflectance of the soil. The RANSAC algorithm was applied in each path separately.
该过程首先对每个点云对实施SEN过滤器，首先计算每个点到其邻居的平均距离，然后拒绝比平均距离更远的点，同时还要考虑标准差。此外，还使用随机样本共识（RASAC）算法来删除属于地面的点。数据点与平面之间的距离阈值定义为5 mm，而平面模型中通过实现各行RTK-GNSS高度差额外考虑了场的倾斜度。如上所述，在两行中都应用了一个5毫米高的平面，以去除土壤的反射率。RASAC算法分别应用于每条路径。

The pairwise pathways were pre-aligned due to the aforementioned georeferencing method of point clouds using the RTK-GNSS. An additional rough pre-alignment was carried out in CloudCompare in order to overlap the pair sides in the common coordinate system. The alignment of the pairwise pathways was operated with the iterative closest point algorithm (ICP), which iteratively minimizes the distances of corresponding points. A binary search tree, where each node represents a partition of the k-dimensional space (K-d tree) was used to speed up the ICP process [28]. Moreover, the local surface normal vector and the curvatures were considered as variants in the ICP to select the corresponding point. A weight factor that incorporates the difference of eigenvalues of matched values for the minimization of distances was assigned to x, y, z, and the normal vector. Higher weight was assigned to the x axis, since it was the axis of the driving direction producing the highest deviation between the point cloud pairs. The maximum correspondence distance was set to 0.1 m. Thus, the ICP minimized the distance between a point and the surface of its corresponding point, or the orientation of the normal. Finally, the two sides were merged together into a single point cloud, using a voxel grid filter (1 mm × 1 mm × 1 mm). A similar methodology was followed for the reconstruction of each individual side of the metal box. Initially, a SOR filter was used in order for the point clouds of each individual side to be denoised. Subsequently, the normal vectors and curvatures were estimated and integrated into the ICP algorithm. Thus, the edges of the sides were aligned and later merged. Each side was divided vertically and horizontally into 6 parts using a 0.3 m step in order to acquire the average height and the average width, respectively. The height was measured by the subtraction of maximum and minimum point of each vertical part for each individual side, whereas the width was obtained by the subtraction of maximum and minimum width of each horizontal part.
由于使用RTK-GNSS的上述点云地理参考方法，两条路径预先对准。还在CloudCompare中进行了额外的粗略预对准，以便在共同坐标系中重叠两侧。两两路径的比对采用迭代最近点算法，迭代最小化对应点的距离。使用二叉树(其中每个节点表示k维空间的一个分区(K-d树))来加速比较方案处理[28]。此外，将局部表面法向量和曲率作为ICP值的变量来选择对应点。将用于最小化距离的匹配值的特征值的差值的权重因子分配给x、y、z和法向量。X轴被赋予更高的权重，因为它是在点云对之间产生最大偏差的行驶方向的轴。最大对应距离被设置为0.1m。因此，ICP最小化了点与其对应点的表面之间的距离，或法线的方向。最后，使用体素网格过滤器(1 mm×1 mm×1 mm)将两个边合并成单个点云。对金属盒子的每一面的重建都采用了类似的方法。最初，使用SOR滤波器来对每个单独侧面的点云进行去噪。然后，估计法向量和曲率，并将其整合到ICP算法中。因此，两条边的边缘对齐了，然后又进行了合并。每一面都用0号在垂直和水平方向上分为6个部分。3 m步，分别获得平均高度和平均宽度。高度是通过减去每个单独边的每个垂直部分的最大和最小点来测量的，而宽度是通过减去每个水平部分的最大和最小宽度来获得的。

Furthermore, the height and the width of each individual extended bar (ΕΒ) were measured. The difference of average height between the EB (ΔEB) was estimated, using four vertical planes to obtain the difference between the minimum and maximum point.
此外，测量了每个单独的延伸条（EB）的高度和宽度。使用四个垂直平面来获得最小点和最大点之间的差，估计EB之间的平均高度差（ΔEB）。

2.6. Tree Row Alignment  2.6.树行对齐

The methodology for determining the stem position relied on the previous study by Vazquez-Arellano and co-workers [24]. A first-order polynomial curve considering the least absolute residuals (LAR) was used to connect the reference stem positions in each row (Equation (13))
确定股骨柄位置的方法依赖于Vazquez-Arellano及其同事之前的研究[ 24]。使用考虑最小绝对剩余（RAL）的一阶多项曲线连接各行中的参考股骨柄位置（方程（13））

𝑦 = 𝛼 × 𝑥 + 𝛽$$y =\alpha \times x +\beta$$

(13)

The slope of the line was utilized by rotating the merged point cloud until it was aligned with the x-axis.
通过旋转合并的点云直到它与x轴对齐来利用直线的斜率。

2.7. Tree Stem Estimation
2.7.树干估计

The generated and rotated point clouds were used to acquire 3D point density histograms using a bivariate approach by pairing the x and y values of every point in the point cloud. The bivariate histogram depicts the frequency of points, which fall into each bin of size 3 cm².
生成和旋转的点云通过配对点云中每个点的x和y值，使用双变量方法获得三维点密度直方图。双变量直方图描绘落入大小为3 cm ² 的每个箱中的点的频率。

A bivariate histogram of the entire row was used to detect the stem position with the highest frequency. A cylinder with the coordinates of the highest peak as a base, with a 0.65 m radius and 4 m height was applied in order to bring a point density histogram within this frame. The point with the highest frequency was considered as the starting point, since it was the one with the highest probability of being a stem. The position of the next stem was estimated in both directions of the x-axis by adding the average intra row distance (0.95 m). At this point, a second cylinder, with the aforementioned dimensions, and the coordinates of the next stem as the center of the base was applied. The point density histogram was used within the cylinder frame to test whether the coordinates of the next stem would remain the same or need to be repositioned. The tree point cloud within the later cylinder was segmented and this process was iteratively applied for all trees of each row.
使用整行的双变量直方图来检测出现频率最高的词干位置。采用半径为0.65m、高度为4m的圆柱体，以最高点的坐标为基点，在该帧内得到点密度直方图。频率最高的点被认为是起点，因为它是茎的概率最高的点。通过加上平均行内距离(0.95m)，在x轴的两个方向上估计下一个茎的位置。在这一点上，应用了第二个圆柱体，具有上述尺寸，并将下一个杆的坐标作为底座的中心。在圆柱体框架内使用点密度直方图来测试下一个杆的坐标是保持不变还是需要重新定位。对后一个圆柱体内的树点云进行分割，并对每行的所有树迭代应用此过程。

2.8. Tree Height  2.8.树高

When each tree was segmented based on its stem position, the height was estimated by calculating the difference between the maximum and the minimum point in the z-axis (H_LiDAR). However, the plane that was fitted into the soil points with the RANSAC algorithm falsly classified some plant points, which were close to the soil, as soil points [29,30]. Therefore, the rigid transformations used for registering the tree point clouds were implemented to align and merge the soil point clouds. The merged soil point cloud was filtered utilizing the SOR filter, and merged with the tree point cloud. The merged point cloud capturing the tree and soil information was utilized for the height estimation.
当每棵树根据其树干位置进行分割时，通过计算z轴上最大点和最小点之间的差(H _LiDAR )来估计高度。然而，用RANSAC算法拟合到土壤点的平面错误地将一些靠近土壤的植物点归类为土壤点[29，30]。为此，采用刚体变换对树木点云进行配准，对土壤点云进行对齐和融合。利用SOR滤波器对合并后的土壤点云进行滤波，并与树状点云进行融合。利用融合后的点云获取树木和土壤信息，进行高度估计。

2.9. Stem Diameter   2.9。茎直径

As was mentioned above, the S_manual was measured at a height of 0.30 m above the ground for each tree. Thus, the region between 0.05 m and 0.3 m above the soil was selected in each individual segmented tree to calculate the stem diameter (S_LiDAR). The selected points were plotted in x and y-axis. A line connected the points on the boundaries with each other, creating a polygon. The k-means clustering method was applied in order to define the coordinates of the center inside the polygon. Consequently, the maximum Euclidean distance between the center and the boundary points were estimated and compared with the manually measured reference stem diameter.
如上所述，S _manual 是在每棵树离地0.30米的高度测量的。因此，在每一株切段树中，选择距土壤0.05m到0.3m之间的区域来计算树干直径(S _LiDAR )。选定的点在x轴和y轴上绘制。一条线将边界上的点彼此连接起来，创建了一个多边形。为了确定多边形内中心的坐标，采用了k-均值聚类法。因此，估计了中心到边界点之间的最大欧几里德距离，并与手动测量的参考茎直径进行了比较。

2.10. Canopy Volume  2.10.冠层体积

The H_manual of the tree and distance of the branching point of the tree from the ground were recorded. Canopy width of each tree at height intervals of 0.5 m was measured from both sides. The cross-sectional area of each zone was calculated by multiplying the mean value of the top and bottom width measured with the corresponding height. As was mentioned above, the product of the area and h_manual of the tree was multiplied to determine the V_manual [4].
记录树木的H和树木分支点距地面的距离。以0.5 m高度间隔从两侧测量每棵树的树冠宽度。通过将测量的顶部和底部宽度的平均值乘以相应的高度来计算每个区域的横截面积。如上所述，将树的面积与h的积相乘以确定V [ 4]。

On the other hand, the convex hull was applied in the segmented tree point clouds. The Delaunnay approach was selected for triangulation, whereas a shrink factor equal to 1 was used to create a compact boundaries that envelops the points. This method fits a boundary around the structure specified and calculates the volume for that region. For measuring volume with this method, an algorithm was developed using the boundary function with same shrink factor accordingly in Matlab [31] to estimate the volume.
另一方面，在分割的树点云中应用了凸壳。选择Delaunnay方法进行三角测量，而使用等于1的收缩因子来创建使点重叠的紧凑边界。此方法适合指定结构周围的边界，并计算该区域的体积。为了用这种方法测量体积，在Matlab [ 31]中使用具有相同收缩因子的边界函数开发了一种算法来估计体积。

3. Results and Discussion
3.结果和讨论

3.1. Measuring Uncertainty of the System Applied in the Field
3.1.现场应用系统的不确定度测量

After reconstructing the sides of the metal box, the LiDAR system was able to detect the slope of the soil, and the structural details of the box (Figure 2). The extended metal bars on one side were depicted, whereas their shape was also shown in the main body of S₁. However, gaps of points were visible due to the sun’s reflection from the metal box. The angular resolution of 0.1667° resulted in an average of 20 hits per mm². In total, 565 planes were recorded on S₁, S₃, and S₄ while 643 planes were recorded in S₂. The ICP algorithm was utilized to align and register the sides together and then merged together (Figure 3b). The C2C function was utilized to estimate the average overlapping of the sides, revealing a mean distance error of 12 mm with a standard deviation of 0.5 mm.
重建金属箱的侧面后，LiDART系统能够检测土壤的斜坡和箱的结构细节（图2）。一侧的延伸金属条被描绘出来，而S的主体中也显示了它们的形状。然而，由于太阳从金属盒的反射，点的间隙是可见的。0.1667°的角分辨率导致平均每毫米撞击20次。S、S和S上总共记录了565架飞机，S上记录了643架飞机。利用ICP算法将边对齐和配准在一起，然后合并在一起（图3b）。利用C2C功能来估计两侧的平均重叠，显示平均距离误差为12 mm，标准差为0.5 mm。

Figure 2. (a) Reconstruction of the four sides from the metal box. S₃ and S₄ are the sides parallel to the y axis, while S₁ and S₂ are the sides parallel to the x axis. (b) The resulting merged point cloud after registration and alignment of the four sides of the metal box.
图2. (a)从金属盒中重建四面。S和S是平行于y轴的边，而S和S是平行于x轴的边。(b)配准和对齐金属盒的四个侧面后产生的合并点云。

Figure 3. Representation of the point clouds were reconstructed when the LiDAR system on the tractor drove (a) through the path a, c on the left side and (b) through the path d, b on the right side of apple trees. For illustration purposes, only a part of each row was presented.
图3.当拖拉机上的LiDART系统（a）穿过苹果树左侧的路径a、c和（b）穿过苹果树右侧的路径d、b时，重建了点云的表示。出于说明目的，仅显示了每一行的一部分。

The highest MAE and RMSE of height estimation was found for S₁ with 8.18 mm and 1.63%, respectively, while the highest MBE of 2.75 mm was noted at the same side (Table 2). For S₄ the lowest MAE and RMSE were depicted with 1.06 mm and 0.21%, respectively. On the other hand, the results of the width showed generally lower MAE, while enhanced MBE were found.
S的高度估计MAE和RMES最高，分别为8.18毫米和1.63%，而同侧的MBE最高，为2.75毫米（表2）。对于S，最低MAE和RSSE分别为1.06毫米和0.21%。另一方面，宽度的结果显示MAE普遍较低，而MBE增强。

Table 2. Height and width results for each side (S) of the metal box (n = 6), maximum (Max), minimum (Min), standard deviation (SD), mean absolute error (MAE), bias (MBE), root mean square error (RMSE).

		Mean (m)	Max (m)	Min (m)	SD (mm)	MAE (mm)	MBE (mm)	RMSE (%)
Height	S1	1.83	1.83	1.81	13.98	8.18	2.75	1.63
	S2	1.80	1.81	1.79	13.34	3.00	1.75	0.60
	S3	1.80	1.82	1.79	1.78	3.31	0.68	0.66
	S4	1.79	1.81	1.79	5.37	1.06	−0.81	0.21
Width	S1	0.60	0.60	0.58	2.10	4.51	2.33	0.90
	S2	0.59	0.60	0.58	7.10	1.50	−1.5	0.30
	S3	0.58	0.59	0.58	7.96	2.06	−4.06	0.81
	S4	0.59	0.61	0.59	8.5	2.18	−2.22	0.43

When regarding the extended bars (EB), the S₁ results are becoming more pronounced. The measuring uncertainty increased, when considering the analysis of height and width of EB. EB located closest to the ground, deviating the most from the actual dimensions of the bars (Table 3). The EB₄ and EB₅ were closer to the LiDAR height level and less gaps were recognized. Consequently, the measuring uncertainty was reduced. The mean height of EB₅ was 29 mm with the lowest MAE of 0.2 mm and −0.2 mm for the MBE. The highest MBEs of 2.8 mm and −2.8 mm and RMSEs of 5.5% and 5.3% were depicted in EB₁ and EB₂, respectively. Similar results were found by Palleja et al. [15], who mentioned that the LiDAR height misposition due to the roughness of the ground and the oscillations of the tractor engine’s trajectory during the operation can influence the measurement outcome, producing errors of up to 3.5%.
当关于延伸条（EB）时，S结果变得更加明显。当考虑EB高度和宽度分析时，测量的不确定度增加。EB位于最靠近地面的位置，与棒的实际尺寸偏差最大（表3）。EB和EB更接近LiDART的高度水平，并且发现的差距更小。因此，测量的不确定度降低了。EB的平均高度为29 mm，最低MAE为0.2 mm，MBE为-0.2 mm。EB和EB中的最高MBE分别为2.8毫米和-2.8毫米，RSSE分别为5.5%和5.3%。Palleja等人[ 15]也发现了类似的结果，他们提到，由于地面粗糙度和运行期间拖拉机发动机轨迹的振荡而导致的LiDART高度错位可能会影响测量结果，产生高达3.5%的误差。

Table 3. Estimated height and width of the extended bars (EB) (n = 6), maximum (Max), minimum (Min), standard deviation (SD), mean absolute error (MAE), bias (MBE), root mean square error (RMSE).

		Mean (mm)	Max (mm)	Min (mm)	SD (mm)	MAE (mm)	MBE (mm)	RMSE (%)
Height	EB1	21	30	11	7.9	2.1	−2.1	4.2
	EB2	10	30	10	8.5	2.8	−2.8	5.6
	EB3	27	29	20	4.8	0.8	−0.8	1.5
	EB4	28	30	20	4.8	0.6	−1.5	3.2
	EB5	29	30	28	0.9	0.2	−0.2	3.0
Width	EB1	89	100	81	8.5	2.8	2.8	5.5
	EB2	81	85	79	2.6	0.4	−2.8	5.3
	EB3	104	120	85	4.8	0.8	−0.8	1.5
	EB4	102	103	100	1.9	1.5	0.4	0.8
	EB5	102	104	100	1.9	1.6	0.4	0.7

The mean distance between EB₁ and EB₂ (ΔEB₁) was 312.67 mm illustrating, among the other ΔEB, the lowest SD (±1.74 mm) and MAE of 1.27 mm, but also the highest MBE 3.1 mm. Furthermore, the ΔEB₃ deviate the most from the known dimensions of the metal box, with a 293 mm mean and ±6.68 mm SD. It should be mentioned that the mean values of ΔEB₂ and the ΔEB₃ were rather close with 295.2 mm and 296.8 mm, depicting a 4.77 mm and 3.17 mm MBE, respectively.
EB和EB之间的平均距离（A EB）为312.67 mm，说明除其他外，最低SD（±1.74 mm）和MAE为1.27 mm，但最高MBE为3.1 mm。此外，A EB与金属盒已知尺寸的偏差最大，平均值为293 mm，SD为±6.68 mm。应该指出的是，A EB和A EB的平均值相当接近，分别为295.2 mm和296.8 mm，分别描绘了4.77 mm和3.17 mm MBE。

The scanning process can basically be affected by the scanner mechanism, the atmospheric conditions and environment, the object properties, and the scanning geometry [32]. The highest MBE and MAE have been denoted in the EB measured at a high incident angle. Consequently, the point cloud quality is affected by the scanning geometry, when increased incidence angles and enhanced distances from the scanner occured [33,34]. While the bias may be corrected, the lack of precision remains. It should be also considered that materials of high reflectance can distort the backscattered signal [32]. This was the case for S₃, where 440 missing hits were found. The distortion will be especially pronounced if part of the beam falls outside the target, the likelihood of which additionally depends on the angular resolution [33]. Therefore, the combination of the perturbating factors can possibly increase the MBE and MAE, especially in objects of smaller dimension such as the EB. Whereas this was an expected finding, the quantification points to the measuring uncertainty that can be expected in the analysis of tree stem.
扫描过程基本上可以受到扫描仪机构、大气条件和环境、物体属性以及扫描几何形状的影响[32]。在高入射角下测得的电子束具有最高的MBE和MAE。因此，当入射角增大和与扫描仪的距离增大时，扫描几何形状会影响点云质量[33，34]。尽管这种偏见可能会得到纠正，但精确度的不足依然存在。还应考虑高反射率的材料会使后向散射信号失真[32]。S ₃ 的情况就是这样，有440个点击率缺失。如果部分光束落在目标之外，失真将特别明显，其可能性还取决于角度分辨率[33]。因此，扰动因素的组合可能会增加MBE和MAE，特别是在诸如EB的较小维度的对象中。虽然这是一个意料之中的发现，但量化指向了在树干分析中可以预期的测量不确定度。

3.2. Separation of Trees by Means of the Stem Position
3.2.通过树干位置分离树木

Each side of the row was scanned from both sides by the LiDAR system. The structure of the left side of the trees was described by the a and c pathways, while the right side was described by b and d (Figure 3a,b).
激光雷达系统从两侧扫描了一行的每一面。树的左侧结构由a和c路径描述，而右侧由b和d路径描述(图3a，b)。

The registration and alignment of the point clouds was carried out pairwise, utilizing the ICP algorithm. The ICP algorithm was utilized in two cases: on the pairwise sides with the soil points (Figure 4a) and on the pairwise sides after the removal of soil points using RANSAC (Figure 4b) in order to evaluate which case provides a better alignment. The C2C function was used to quantify the overlapping differences. Alignment including the soil points resulted in a distance error of 163.5 mm with 31.1 mm SD, and 80.5 mm MBE was observed. The removal of soil points by contrast led to a reduced mean distance error of 81.6 mm with a 21.4 mm SD and lower RMSE of 0.71% as well as an MBE of 5.2 mm. As a result, for the case of the removal of soil points, low measuring uncertainty affecting the overlapping process was found [24,35].
点云的配准和对准是利用ICP算法两两进行的。在两种情况下使用了比较方案算法：在与土壤点成对的一侧(图4a)和在使用RANSAC去除土壤点之后的成对侧(图4b)，以评估哪种情况提供了更好的对准。采用C2C函数对重叠差异进行量化。包括土壤点的排列导致距离误差为163.5 mm，标准偏差为31.1mmSD，并观察到80.5mmMBE。相比之下，去除土壤点导致的平均距离误差减少了81.6 mm，SD为21.4 mm，RMSE为0.71%，MBE为5.2 mm。结果，对于去除土壤点的情况，发现影响重叠过程的低测量不确定度[24，35]。

Figure 4. Registration and alignment using the iterative closest point algorithm (ICP) algorithm of a single tree point cloud in two cases (a) tree point cloud with soil points (b) tree point cloud after soil removal with RANSAC. The trees with the best results for each case are presented.
图4.在两种情况下，使用单个树点云的迭代最近点算法（ICC）算法进行配准和对齐（a）具有土壤点的树点云（b）使用RASAC去除土壤后的树点云。给出了每个情况下具有最佳结果的树。

The point clouds of both rows were merged into a single point cloud. The structure of the trees was clearly defined. However, it was observed that weed points above 5 mm remained within the rows (Figure 4b).
两行的点云合并为单个点云。树木的结构被明确定义。然而，观察到5毫米以上的杂草点仍留在行内（图4 b）。

The raw data of the 3D reconstructed point cloud consisted of 18,441,933 points. The application of the SOR filter reduced the point cloud by 12.17% removing 2,241,786 points. The segmentation of the soil points using RANSAC algorithm revealed that a further 58.59% of the data (10,820,983 points) belonged to the soil. The remained 5,382,849 points belonged to the trees representing 29.23% of the total amount of points. Arellano and co-workers [30] already pointed out the dominance of soil hits in their 3D reconstructed point clouds, suggesting that the alignment of plant points is separate from the soil points.
3D重建点云的原始数据由18，441，933个点组成。SEN过滤器的应用使点云减少了12.17%，删除了2，241，786个点。使用RASAC算法对土壤点进行分割后发现，另外58.59%的数据（10，820，983个点）属于土壤。剩余的5，382，849个积分属于树木，占积分总数的29.23%。Arellano和同事[ 30]已经指出了土壤撞击在他们的3D重建点云中的主导地位，这表明植物点的排列与土壤点是分开的。

The bivariate point density histogram enabled the detection of the peak of laser hits for each individual tree (Figure 5). The assumption that the stem points are concentrated in the center of the plant was based on previous studies [8,24]. Consistently, the area closer to the stem position was assumed to appear with enhanced frequency. However, the assumption needed confirmation considering the formation of the tree canopy and, furthermore, the analysis of number and position of peaks can be affected by the quality of the point cloud. The potential noise or the low quality of overlapping of the point cloud could lead to more peaks since the local bin count becomes more uneven. The use of a stable bin size triggered the filtering of the peaks in the bivariate histogram and, therefore, affected the detection of the highest peak.
二元点密度矩形图能够检测每棵树的激光撞击峰值（图5）。茎点集中在植物中心的假设是基于之前的研究[ 8，24]。一致地，假设更接近茎位置的区域出现频率增加。然而，考虑到树冠的形成，该假设需要确认，此外，峰值数量和位置的分析可能会受到点云质量的影响。潜在的噪音或点云重叠的低质量可能会导致更多峰值，因为局部bin计数变得更加不均匀。稳定的bin大小的使用触发了二元矩形图中峰的过滤，因此影响了最高峰的检测。

Figure 5. Result of the registration and alignment of one apple row point cloud after merging and filtering. (a) Intensity image, where the warmer squares indicate the peaks and (b) the bivariate point density histogram of the resulted row point cloud.
图5.合并和过滤后一个苹果行点云的配准和对齐结果。(a)强度图像，其中较暖的正方形指示峰值，以及（b）所得行点云的二元点密度矩形图。

The coordinates of the estimated stem position were compared with the stem position, which was defined with the RTK-GNSS (Figure 6). The estimated stem positions were closely related to the RTK stem positions, revealing 33.7 mm MAE with 36.5 mm MBE. Despite the fact that the alignment occurred with a certain error, the algorithm was able to miss only one stem position out of 224 trees (Figure 6).
将估计的茎位置的坐标与使用RTK-GNSS定义的茎位置进行比较（图6）。估计的股骨柄位置与TEK股骨柄位置密切相关，显示MAE为33.7 mm，MBE为36.5 mm。尽管对齐存在一定误差，但该算法只能错过224棵树中的一个茎位置（图6）。

Figure 6. Estimation of the tree (n = 224) stem position compared to the RTK-GNSS defined stem position.
图6.与RTK-GNSS定义的树干位置相比，估计树木（n = 224）树干位置。

3.3. Estimation of Tree Variables
3.3.树变量的估计

After the determination of stem position, the subtracted soil was merged with the point cloud of trees (Figure 7a) to avoid loosing tree hits. The coordinates of the estimated stem position were utilized as the center for the segmentation cylinders in order to obtain the points that belonged to each individual tree. The points within the boundaries of the cylinder were segmented and considered as the tree points (Figure 7b).
确定树干位置后，将减去的土壤与树木的点云合并（图7a），以避免失去树木撞击。估计的茎位置的坐标被用作分割圆柱体的中心，以获得属于每棵树的点。圆柱体边界内的点被分割并视为树点（图7 b）。

Figure 7. (a) Merged point cloud with estimated stem position and stem position defined by RTK-GNSS. For illustration purposes, only first part of row 1 is represented. (b) Projected cylinders for single tree segmentation based on the estimated stem position.
图7. (a)具有估计干位置和RTK-GNSS定义的干位置的合并点云。出于说明目的，仅表示第1行的第一部分。(b)基于估计的茎位置进行单树分割的投影圆柱体。

The difference between the maximum and the minimum point in the z-axis was considered to estimate the height. As was mentioned above, the radius of the cylinder was based on the manually measured average width of the trees (n = 224). However, some tree points of long branches were laid out of the cylinder, resulting in an underestimation of the canopy volume. Nevertheless, this error did not reveal an underestimation of the height (Table 4). Arellano et al. [24], with a similar method, determined the stem position in maize and acquired the height profile of the plants pointing to a MAE of 7.3 mm from the actual height. Walnuts trees were segmented and their height was estimated considering the Min and Max values in the z axis, resulting in a significant relationship with the manual measurements (R² = 0.95) [36]. Malambo et al. [37] utilized a similar methodology with cylinder boundaries to segment sorghum trees, reporting a correlation of r = 0.86 with RMSE = 11.4 cm.
考虑z轴上最大点和最小点之间的差来估计高度。如上所述，圆柱体的半径基于手动测量的树木平均宽度（n = 224）。然而，一些长枝的树点被布置在圆柱体之外，导致树冠体积被低估。然而，这一错误并没有表明对身高的低估（表4）。Arellano等人[ 24]采用类似的方法确定了玉米的茎位置，并获得了植物的高度轮廓，表明MAE与实际高度的7.3毫米。核桃树被分割，并考虑z轴上的Min和Max值估计其高度，结果与手动测量存在显着关系（R = 0.95）[ 36]。Malambo等人[ 37]利用类似的圆柱边界方法来分段高粱树，报告r = 0.86与RSSE = 11.4厘米的相关性。

Table 4. Descriptive statistics of manual measurements of height (H_manual) (mm), stem diameter (S_manual) (mm), canopy volume (V_manual) (mm³) and LiDAR-based measurements H_LiDAR, S_LiDAR, V_LiDAR (n = 224) regarding maximum (Max), minimum (Min), mean absolute error (MAE) (mm), bias (MBE) (mm), root mean square error (RMSE), coefficient of determination (R²).

	Min	Max	Mean	MAE	MBE	RMSE (%)	R2
Hmanual (mm)	1900	2800	2310	5.55	0.62	5.71	0.87
HLiDAR (mm)	1870	2820	2350
Smanual (mm)	55	132.5	97.1	2.52	−3.75	2.23	0.88
SLiDAR (mm)	63	136	99.5
Vmanual (m3)	0.23	1.12	0.55	5.23	5.93	4.64	0.77
VLiDAR (m3)	0.38	1.05	0.58

In the present study, the plant height profile of each row was measured manually, as the ground-truth in order to validate the results of the point cloud, pointing out a 900 mm variation in height (Table 4). The average values of manual measurements (H_manual) and estimated height, which was derived from LiDAR data (H_LiDAR), varied slightly by 40 mm when comparing the means. Similar values can be expected in commercial apple orchards. The reference and estimated values correlated, revealing an R² = 0.87 with a MAE of 5.55 mm and MBE of 0.62 mm. The RMSE of 5.7% can be explained by the measuring uncertainty of the LiDAR, as well as the uncertainty of defining the highest point of a the stem elongation by means of manual measurements.
在本研究中，手动测量了每一行的植物高度轮廓，作为地面事实，以验证点云的结果，指出高度变化900毫米（表4）。在比较平均值时，手动测量的平均值（H）和根据LiDART数据得出的估计高度（H）的平均值略有变化40 mm。商业苹果园也可以预期类似的价值。参考值和估计值相关，显示R = 0.87，MAE为5.55 mm，MBE为0.62 mm。5.7%的RSSE可以通过LiDART的测量不确定度以及通过手动测量定义股骨柄延伸最高点的不确定度来解释。

The section of stem between 0.10 m and 0.30 m was selected in the segmented tree point clouds to analyse the stem diameter (Figure 8). The grafting region of the cultivar and rootstock appeared always below this section. The S_manual and S_LiDAR depict 2.4 mm difference in their means, whereas an 8 mm difference in minimum values was illustrated. The latter results were described by a considerable difference in the MAE of 2.52 mm and the MBE of −3.75 mm, however, reference and estimated values showed high R² = 0.88 and relatively low 2.23% RMSE recognizing the pertubating effects of enhanced incident angle when measuring on the extended bars of the box (Table 4). The cylinder fitting method based on a stem position for estimating the stem diameter is rather scarce in horticulture. The automatic detection that fits cylinders on stem points, has been investigated broadly in forestry [38,39,40]. However, the majority of these studies employed a 3D terrestrial LiDAR, which is able to provide a 3D representation of stem points without moving the system and as a consequence the limitation of errors during the alignment. Moreover in forests, the tree stem has no coverage from leaves from a few meters above the ground, allowing more hits per stem and segmentation of stem points without imposing coinciding structures.
在分段的树木点云中选择0.10米到0.30米之间的树干截面来分析树干直径(图8)。品种和砧木的嫁接区域总是出现在这一区段以下。S _manual 和S _LiDAR 的平均值相差2.4 mm，而最小值相差8 mm。后者的结果由MAE为2.52 mm和−为3.75 mm的相当大的差异来描述，然而，参考值和估计值显示较高的R ² =0.88和相对较低的2.23%RMSE，当在盒子的延伸杆上测量时，认识到增强的入射角的扰动效应(表4)。园艺中基于茎位置的圆柱体拟合法估算茎直径的方法比较少见。将圆柱体安装在树干尖端的自动检测在林业中得到了广泛的研究[38，39，40]。然而，这些研究中的大多数使用了3D陆地LiDAR，它能够在不移动系统的情况下提供茎点的3D表示，从而限制了对准过程中的误差。此外，在森林中，树干在离地面几米的地方没有树叶的覆盖，允许每个树干更多的命中率和树干点的分割，而不会强加重合的结构。

Figure 8. (a) Stem point selection marked in red between 0.1 m and 0.3 m and (b) selected stem points detected by the LiDAR (S_LiDAR) plotted in 2D with stem boundaries and stem position.
图8. (a)干点选择在0.1 m和0.3 m之间用红色标记，以及（b）由LiDART（S）检测到的选定干点，以2D方式绘制，并带有干边界和干位置。

The convex hull algorithm was used to estimate the canopy volume (Figure 9b). The boundary for each individual volume was illustrated in false colour in order to be easily distinct. The convex hull method does not compute volume of the gaps within the canopy, similar to ground-truth methods but instead considers the overall canopy geometry. The manual measurement of the volume showed a strong correlation with the volume data from the segmented tree point clouds (R² = 0.77) (Table 4). The rough simplification of the canopy structure by the manual method based on the cylinder-fit method mainly overestimated the actual data and the ones computed with the convex-hull algorithm. Consequently, the highest bias was found for the volume analysis (Table 4), while a slightly reduced RMSE was observed compared to the height estimation.
使用凸壳算法来估计树冠体积（图9 b）。为了易于区分，每个单独体积的边界都以假色示出。与地面真实方法类似，凸壳方法不会计算树冠内间隙的体积，而是考虑整体树冠几何形状。体积的手动测量显示与来自分段树点云的体积数据具有强相关性（R = 0.77）（表4）。基于圆柱体配合法的手工方法对树冠结构进行粗略简化，主要高估了实际数据和凸壳算法计算的数据。因此，发现体积分析的偏差最高（表4），而与高度估计相比，观察到RMES略有降低。

Figure 9. (a) Point cloud of four example apple trees and (b) visualization of convex hull method showing the canopy volume (V_LiDAR,T) of singularized tree (T_n).
图9. (a)四棵示例苹果树的点云和（b）凸壳方法的可视化，显示奇异树（T）的树冠体积（V）。

This overestimation of volume can be a result of weed points that fall within the cylinder and have not been removed during the soil removal, consequently the convex-hull boundaries were increased which in turn provided an overestimation. In a former study, the convex-hull approach was utilized for the determination of volume in apple trees and vines, revealing a strong correlation of r = 0.81 and r = 0.82, respectively [2]. Furthermore, it was mentioned that the convex-hull approach does not consider the canopy gaps and can lead to an overestiamation. Cheein et al. [13] developed a cylinder based model for real time characterization of volume in apple trees, segmenting the trees based on their canopy center. The model showed 87% accuracy, however the points in the periphery of the cylinder can cause an exaggerated and unrealistic increase in the calculated volume. The methodolology that was proposed in this study for tree segmentation, based on their stem position, revealed a strong relationship between the V_manual and V_LiDAR (R² = 0.77). However, in some trees, the main drawback was that within the cylinder, tree points from neighbor trees interfered, leading to a possible bias of the volume.
这种对体积的高估可能是由于落在圆柱体内的杂草点在土壤清除过程中没有被去除，因此凸壳边界被增加，这反过来又提供了高估。在以前的研究中，用凸壳方法测定苹果树和藤本植物的体积，其相关性分别为r=0.81和r=0.82[2]。此外，有人提到，凸壳方法没有考虑树冠间隙，可能会导致过度估计。Cheein等人。[13]开发了一种基于圆柱体的模型，用于实时描述苹果树的体积，根据树冠中心对树木进行分割。该模型显示了87%的准确率，但是圆柱体外围的点可能会导致计算体积的夸大和不切实际的增加。这项研究中提出的基于树干位置的树木分割方法表明，V _manual 和V _LiDAR 之间存在很强的关系(R ² =0.77)。然而，在一些树木中，主要的缺点是在圆柱体内，来自相邻树木的树点干扰，导致可能的体积偏差。

4. Conclusions  4.结论

The LiDAR system was able to estimate a simple structure of objects, with the highest RMSE of 1.63 in a situation with missing hits due to reflectance of the object and an increased incident angle of laser scanner and hit density. Application of a bivariate point density histogram allowed the detection of the tree stem position with an MAE of 43.7 mm and MBE of 36.5 mm. It also provided meaningful information about the H_LiDAR, which correlated strongly with the H_manual (R² = 0.87). Furthermore, the approach for the estimation of S_LiDAR indicated a high relationship with the S_manual (R² = 0.88). The V_manual and the V_LiDAR revealed a lower correlation in comparison with the above parameters (R² = 0.77), probably due to errors in V_manual or the presence of weeds within the trees in the point cloud data. The convex-hull method for estimating the canopy volume was confirmed in apple trees.
LiDAR系统能够估计简单的物体结构，在由于物体的反射和激光扫描器入射角度和命中密度增加而导致命中失误的情况下，最高RMSE为1.63。应用双变量点密度直方图可以检测树干位置，MAE为43.7 mm，MBE为36.5 mm。它还提供了关于H _LiDAR 的有意义的信息，这与H _manual (R ² =0.87)密切相关。此外，对S _LiDAR 的估计方法与S _manual 有很高的相关性(R ² =0.88)。与上述参数(R ² =0.77)相比，V _manual 和V _LiDAR 的相关性较低(R ² =0.77)，这可能是由于V _manual 中的错误或点云数据中树内存在杂草。在苹果树上证实了凸壳估测冠层体积的方法。

Further research needs to be done using clustering algorithms for single plant segmentation (e.g., Euclidian, max-flow/min-cut theorem, region growing segmentation) so that parts of the canopy are not cut off from the tree point cloud, and plants at a later growth stage can be analyzed. Also, non-rigid registration and alignment need to be explored to reduce the effect of the miscellaneous errors that reduce the overlapping between point cloud pairs.
需要使用用于单植物分割的聚类算法(例如，欧几里德、最大流/最小切割定理、区域生长分割)进行进一步的研究，以便树冠的部分不会从树的点云中被切断，并且可以分析后期生长阶段的植物。此外，需要探索非刚性配准和对齐，以减少各种误差的影响，这些误差减少了点云对之间的重叠。

In summary, the LiDAR system was able to detect geometric variables in apple trees that can be used in agricultural applications to measure tree growth for tree-individual orchard management, while considering mechanical pruning, irrigation, and spraying.
总而言之，LiDART系统能够检测苹果树中的几何变量，这些变量可用于农业应用，以测量树木生长，以进行树木个体果园管理，同时考虑机械修剪、灌溉和喷洒。