Estimation of Crop Height Distribution for Mature Rice Based on a Moving Surface and 3D Point Cloud Elevation

: Estimation of rice plant height distribution plays a signiﬁcant role in keeping the feed rate of rice combine harvesters stable. This is an effective way to guarantee the working stability of the whole machine, as a consequence, improving threshing and cleaning efﬁciency and reducing loss and damage rates. However, dense growth and leafy and bent branches of mature rice make it difﬁcult to detect the lowest point of aggregated growing plants in three dimensional (3D) point cloud data. Therefore, the objective of this study was to put forward a method to estimate plant height distribution on the basis of a moving surface and 3D point cloud elevation. The statistical outlier removal (SOR) algorithm was used to reduce noise points far away from target point cloud body, and then moving surface ﬁtting elevation was applied to achieve accurate classiﬁcation of ground and crop point cloud data for plant height estimation. Experiments showed that, compared with the actual value, the average square root error (RMSE) of the estimation results was 8.29, the average absolute percentage error (MAPE) was 9.28%, and the average accuracy was 90%. The proposed method could accurately estimate the height of mature rice and is beneﬁcial to calculating the feed rate in advance, which can provide a reference for further investigation in automatic and intelligent harvesting.


Introduction
Rice is an important grain crop that is highly appreciated worldwide owing to its high nutritive value; the use of rice combine harvesters has been increasing rapidly with the increase in the planting area and yield in recent years [1]. Stable feed rate technology plays an important role in maximizing the capacity of a combine harvester and guaranteeing the performance parameters when harvesting, including the impurity rate, broken rate, and loss rate. However, the crop morphology varies spatially; generally the crop height, cultivating density, and yield present a great difference even in the same field, which has become the overriding reason that operators are typically unable to adjust the driving speed on time and therefore are unable to maintain a stable feed rate. In order to assist drivers, an automatic feed rate control system was developed by fixing internal sensors to the machine to constantly measure the crop intake [2]. Since crops have already been cut and are in the delivery stream, a time delay is unavoidable, and it could be impossible to take actions for an instant adjustment. Therefore, the estimation of the feed rate in advance could be a significant help to support drivers, so as to adjust the working parameters of combine harvesters to cope with the variable harvesting conditions [3][4][5]. Plant height of crops in the area to be harvested is one of essential influential factors, and the estimation in advance could provide necessary information for the control of working parameters under optimum conditions [6,7]. 2 of 14 Clumping mature rice in the field has the characteristics of luxuriant branches and leaves, and their panicles are usually bent and intersected. Rice harvesting usually causes plenty of dust and MOG (Materials other than Grains). Many methods for estimation of crop canopy attributes depend on ultrasonic ranging, machine vision, and LiDAR (light detection and ranging). Maertens et al., [8] adopted ultrasonic sensors and investigated crop density measurement systems, and the experiments indicated that the reflection of ultrasonic waves seemed highly dispersed and could be easily affected by wind, dust, and crop distribution. Munch et al., [9] developed a 3D stereo-camera in order to provide point cloud data for crop volume estimation; however, the system was sensitive to the dust conditions and failed to consistently obtain reliable data.
In contract, LiDAR has the advantages of high resolution and seems not sensitive to light intensity; it does not suffer from the above disadvantages. More recent research has been focused on using LiDAR for phenotyping purposes [10]. Saeys et al., [11] estimated the density of small grains relying on the 3-D reconstruction of LiDAR technology, obtaining a good result of an R 2 of 0.8. Blanquart et al., [12] expanded Saeys's work and analyzed the influence of laser mounting position and robust online measurement of height and density from a single LiDAR scan.
Recent literature showed that LiDAR technology has aroused great interest to estimate plant height. It seems that once the crop vertex (the highest point of the crop canopy) and ground profile can be identified, the crop height could be calculated through the height difference between them. Zhang et al., [13] designed a laser scanner system for the height of the stem of giant miscanthus and acquired the elevation value of the miscanthus point cloud data, the ground points, and the highest point of the crop. The least-squares surface fitting algorithm was used to generate the top and ground areas of the coverage area so as to successfully estimate the average crop height in the coverage area with an average margin of 5.08%. Hämmerle et al., [14] established a Crop Surface Model (CSM) using crop point cloud data, and the difference between the average CSM elevation difference and the crop height was less than 4%. Zhu et al., [15] collected the 3D point cloud data of field corn through an airborne LiDAR measuring system and generated the height map of the crop in the target area according to data processing and visualization. Phan et al., [16] used a laser scanner to obtain 3D point cloud data of rice at the seedling stage by determining the reference position on the top of the rice, estimating the height of the rice plant as the difference between the top and bottom positions of the rice plant. Anice et al., [17] used a LiDAR sensor to detect the canopy height and width in vineyards and proposed a Bayesian point cloud algorithm to automate the filtering and classifying of 2D LiDAR data. The above algorithms could be effective for immature crops or single crops, but the application to luxuriant crops that cover the ground still need further investigation in that the ground profile can hardly be identified.
In recent years, experts made great efforts with automatic systems to obtain the crop height in large fields. Tilly et al., [18] established a CSM with a resolution of 10 mm, which relied on terrestrial laser scanning. The ground surface is obtained before the crop emerges, and the information is continuously updated during the growth of different crops. Then, crop height can be calculated, but the information collection process needs a rather long time. Jimenez-Berni et al., [19] mounted the LiDAR on a mobile platform for rapid multi-temporal and non-destructive estimation of canopy height, ground cover, and above-ground biomass. However, this method is intended for ridge-planting types, meaning the ground can be clearly seen. Sritarapipat et al., [20] proposed an automatic detection method for rice plant height calculation based on a digital camera; the average relative error of this method reached 3.34%. However, the markers that provide a height reference are required to be fixed in the field, and the system must be frequently compared with the scale of the marker.
Based on the above, LiDAR technology has been applied to the estimation of crop height, which shows a modern 3D scanning method that could describe the crop biomass information. However, the luxuriant canopy of mature crops with larger density brings the problem of acquiring ground profile information, leading to the traditional data process failing to analyze crop height. Therefore, the objective of this study is to investigate the potential of using a LiDAR measuring system to estimate the plant height of mature rice, and the moving surface engaged in extracting ground information was developed. Finally, the system was tested under field conditions to examine the estimation accuracy. In order to acquire the 3D point cloud data of the mature rice, the LiDAR measuring system established in this study mainly includes three parts: mobile measurement terminal, base station, and computer-controlled terminal ( Figure 1). The VLP-16 (Velodyne, San Jose, CA, USA) laser scanner (parameters shown in Table 1) and Span-IGM-A1 (NovAtel, Beijing, China) integrated navigation module is installed on the mobile measurement terminal. The laser scanner is used to measure the distance between the laser launch point and the reflection point of the ground object. The integrated Navigation module mainly includes the Global Navigation Satellite System (GNSS) receiver and Inertial Measurement Unit (IMU). The GNSS receiver with dynamic differential includes the main and secondary receiving antenna, which is used to determine the spatial position change of the laser scanner during data acquisition. The IMU is used to measure the 3D tilt orientation of the laser scanner. As the measuring system is installed on the guide rail to move and acquire data, there is only slight fluctuation. GNSS and IMU can obtain the position and 3D tilt orientation information of the measuring system in real time, so the fluctuation effect can be eliminated and the measuring accuracy guaranteed [21,22]. are required to be fixed in the field, and the system must be frequently compared with the scale of the marker.

Materials and
Based on the above, LiDAR technology has been applied to the estimation of crop height, which shows a modern 3D scanning method that could describe the crop biomass information. However, the luxuriant canopy of mature crops with larger density brings the problem of acquiring ground profile information, leading to the traditional data process failing to analyze crop height. Therefore, the objective of this study is to investigate the potential of using a LiDAR measuring system to estimate the plant height of mature rice, and the moving surface engaged in extracting ground information was developed. Finally, the system was tested under field conditions to examine the estimation accuracy.

LiDAR Measuring System Construction
In order to acquire the 3D point cloud data of the mature rice, the LiDAR measuring system established in this study mainly includes three parts: mobile measurement terminal, base station, and computer-controlled terminal ( Figure 1). The VLP-16 (Velodyne, USA) laser scanner (parameters shown in Table 1) and Span-IGM-A1 (NovAtel, China) integrated navigation module is installed on the mobile measurement terminal. The laser scanner is used to measure the distance between the laser launch point and the reflection point of the ground object. The integrated Navigation module mainly includes the Global Navigation Satellite System (GNSS) receiver and Inertial Measurement Unit (IMU). The GNSS receiver with dynamic differential includes the main and secondary receiving antenna, which is used to determine the spatial position change of the laser scanner during data acquisition. The IMU is used to measure the 3D tilt orientation of the laser scanner. As the measuring system is installed on the guide rail to move and acquire data, there is only slight fluctuation. GNSS and IMU can obtain the position and 3D tilt orientation information of the measuring system in real time, so the fluctuation effect can be eliminated and the measuring accuracy guaranteed [21,22].    Data acquisition process: mobile measurement terminal, GNSS base station, and PC terminal should be connected to the hardware. Li-acquire software (Beijing Digital Green Earth Technology Co., Ltd., Beijing, China) on the computer is turned on to collect data. The PC terminal communicates with the mobile measurement control terminal. The base station data are collected in advance about ten minutes before the communication is established, and the obtained information is saved on an SD card. Whether the data recording is normal could be checked through the serial port debugging assistant. After the communication is successfully established, start to collect IMU data and GNSS data, and wait for about five minutes to complete the initial calibration. When the azimuth accuracy is high (<0.1), start collecting laser scanner data. IMU measures the altitude information during LiDAR measurement, including the yaw angle, pitch angle, and roll angle [23]. The GNSS receiver provides precise location information of the LiDAR. Then, according to the altitude, position, and distance information, the accurate 3D space coordinates (x, y, z) of each laser footprint are calculated. Finally, the built-in function of Li-acquire software is used to synchronize the data obtained by the laser scanner with the position information obtained by GNSS and the 3D tilt orientation information obtained by the IMU according to the time stamp to obtain 3D point cloud data [24].
In order to analyze the data precisely, an initialization coordinate system is artificially defined before importing the acquired point cloud data. Under this coordinate system, a 10 × 10 grid is divided within the range of 1 × 1 m, and the acquired point cloud data are translated and rotated so that the point cloud coordinate system overlaps the initialization coordinate system. At this point, the point cloud data are divided into 10 × 10 small blocks by the grid. The process is shown in Figure 2.
Data acquisition process: mobile measurement terminal, GNSS base station, and PC terminal should be connected to the hardware. Li-acquire software (Beijing Digital Green Earth Technology Co., Ltd., China) on the computer is turned on to collect data. The PC terminal communicates with the mobile measurement control terminal. The base station data are collected in advance about ten minutes before the communication is established, and the obtained information is saved on an SD card. Whether the data recording is normal could be checked through the serial port debugging assistant. After the communication is successfully established, start to collect IMU data and GNSS data, and wait for about five minutes to complete the initial calibration. When the azimuth accuracy is high (<0.1), start collecting laser scanner data. IMU measures the altitude information during LiDAR measurement, including the yaw angle, pitch angle, and roll angle [23]. The GNSS receiver provides precise location information of the LiDAR. Then, according to the altitude, position, and distance information, the accurate 3D space coordinates (x, y, z) of each laser footprint are calculated. Finally, the built-in function of Li-acquire software is used to synchronize the data obtained by the laser scanner with the position information obtained by GNSS and the 3D tilt orientation information obtained by the IMU according to the time stamp to obtain 3D point cloud data [24].
In order to analyze the data precisely, an initialization coordinate system is artificially defined before importing the acquired point cloud data. Under this coordinate system, a 10 × 10 grid is divided within the range of 1 × 1 m, and the acquired point cloud data are translated and rotated so that the point cloud coordinate system overlaps the initialization coordinate system. At this point, the point cloud data are divided into 10 × 10 small blocks by the grid. The process is shown in Figure 2.

Data Acquisition of Field-Mature Rice Point Cloud
The experiment was carried out in a field near Jingkou District, Zhenjiang, Jiangsu Province (32°10′8″ N, 119°35′30″ E). The rice variety "Nanjing 9108" (Jiangsu Academy of Agricultural Sciences, China) was planted by broadcast sowing. The sowing time was in mid-May, and measurement started in mid-October when the crop was mature. The average panicle length was 112.7 mm, the yield was 9023.4 kg/ha, and the 1000-grain weight

Data Acquisition of Field-Mature Rice Point Cloud
The experiment was carried out in a field near Jingkou District, Zhenjiang, Jiangsu Province (32 • 10 8" N, 119 • 35 30" E). The rice variety "Nanjing 9108" (Jiangsu Academy of Agricultural Sciences, Nanjing, China) was planted by broadcast sowing. The sowing time was in mid-May, and measurement started in mid-October when the crop was mature. The average panicle length was 112.7 mm, the yield was 9023.4 kg/ha, and the 1000-grain weight was 26.2 g. The rice was harvested with 62% stem moisture content and 21% grain moisture content.
A sampling frame of 1 × 1 m was divided into 100 grids in a pattern of 10×10 using a soft rope (as shown in Figure 3). The amount and size of grids are consistent with those of a manually defined initialization coordinate system for convenient comparison between Agronomy 2022, 12, 836 5 of 14 estimated and actual values. The sampling was used to randomly select five groups of field-mature rice samples (Figure 4). At the same time, the average plant height in each small grid was measured manually to compare the accuracy of subsequent detection. Point cloud data of each sample were obtained by the constructed LiDAR measuring system and compared with the measured rice plant height. was 26.2 g. The rice was harvested with 62% stem moisture content and 21% grain moisture content.
A sampling frame of 1 × 1 m was divided into 100 grids in a pattern of 10×10 using a soft rope (as shown in Figure 3). The amount and size of grids are consistent with those of a manually defined initialization coordinate system for convenient comparison between estimated and actual values. The sampling was used to randomly select five groups of field-mature rice samples (Figure 4). At the same time, the average plant height in each small grid was measured manually to compare the accuracy of subsequent detection. Point cloud data of each sample were obtained by the constructed LiDAR measuring system and compared with the measured rice plant height.

Point Cloud Classification Based on Moving Surface Fitting Elevation
The plant height distribution of mature rice in the field was estimated based on the ground and crop point cloud data after accurate classification. During collection of 3D point cloud data of crops by LiDAR measuring systems, the data will inevitably produce extreme points due to the environment change, inherent error of equipment, and the error caused by personnel operation. Actually, extreme points would be unavoidable under bright light conditions or they could appear when the laser hits the edge of the crop leaf. The extreme points could also be understood as noise points, in that they are far from the target point cloud data and need be eliminated to increase measurement accuracy. The existence of noise points makes it difficult to classify 3D point cloud data precisely [25]. Hence, in this study, the statistical outlier removal (SOR) algorithm was used to denoise point cloud data before the ground and crop point cloud classification of mature rice in the field. The moving surface fitting elevation was used to accurately classify the ground and crop point clouds. It laid a foundation for further estimation of rice plant height distribution based on post-classification point cloud elevation.

D Point Cloud Denoising
Noise points away from the target point cloud body in 3D point cloud data are reduced based on the SOR algorithm [26]. The specific process based on the SOR algorithm is as follows: Agronomy 2022, 12, x FOR PEER REVIEW 5 of 14 was 26.2 g. The rice was harvested with 62% stem moisture content and 21% grain moisture content. A sampling frame of 1 × 1 m was divided into 100 grids in a pattern of 10×10 using a soft rope (as shown in Figure 3). The amount and size of grids are consistent with those of a manually defined initialization coordinate system for convenient comparison between estimated and actual values. The sampling was used to randomly select five groups of field-mature rice samples (Figure 4). At the same time, the average plant height in each small grid was measured manually to compare the accuracy of subsequent detection. Point cloud data of each sample were obtained by the constructed LiDAR measuring system and compared with the measured rice plant height.

Point Cloud Classification Based on Moving Surface Fitting Elevation
The plant height distribution of mature rice in the field was estimated based on the ground and crop point cloud data after accurate classification. During collection of 3D point cloud data of crops by LiDAR measuring systems, the data will inevitably produce extreme points due to the environment change, inherent error of equipment, and the error caused by personnel operation. Actually, extreme points would be unavoidable under bright light conditions or they could appear when the laser hits the edge of the crop leaf. The extreme points could also be understood as noise points, in that they are far from the target point cloud data and need be eliminated to increase measurement accuracy. The existence of noise points makes it difficult to classify 3D point cloud data precisely [25]. Hence, in this study, the statistical outlier removal (SOR) algorithm was used to denoise point cloud data before the ground and crop point cloud classification of mature rice in the field. The moving surface fitting elevation was used to accurately classify the ground and crop point clouds. It laid a foundation for further estimation of rice plant height distribution based on post-classification point cloud elevation.

D Point Cloud Denoising
Noise points away from the target point cloud body in 3D point cloud data are reduced based on the SOR algorithm [26]. The specific process based on the SOR algorithm is as follows:

Point Cloud Classification Based on Moving Surface Fitting Elevation
The plant height distribution of mature rice in the field was estimated based on the ground and crop point cloud data after accurate classification. During collection of 3D point cloud data of crops by LiDAR measuring systems, the data will inevitably produce extreme points due to the environment change, inherent error of equipment, and the error caused by personnel operation. Actually, extreme points would be unavoidable under bright light conditions or they could appear when the laser hits the edge of the crop leaf. The extreme points could also be understood as noise points, in that they are far from the target point cloud data and need be eliminated to increase measurement accuracy. The existence of noise points makes it difficult to classify 3D point cloud data precisely [25]. Hence, in this study, the statistical outlier removal (SOR) algorithm was used to denoise point cloud data before the ground and crop point cloud classification of mature rice in the field. The moving surface fitting elevation was used to accurately classify the ground and crop point clouds. It laid a foundation for further estimation of rice plant height distribution based on post-classification point cloud elevation.

D Point Cloud Denoising
Noise points away from the target point cloud body in 3D point cloud data are reduced based on the SOR algorithm [26]. The specific process based on the SOR algorithm is as follows: (1) Determine neighborhood radius r and the standard deviation-multiple threshold m.
(2) Traverse all point cloud data x i within the neighborhood radius r, and calculate the average distance d i between all point cloud data; the calculation of d is as Equation (1): (3) Calculate the mean u and standard deviation σ of each point; the calculation equation is shown in Equations (2) and (3)  Under massive experimental verification, when the neighborhood radius r = 5 and the standard-deviation multiple threshold m = 2, the noise point removal based on the SOR algorithm results in efficient noise removal for the point cloud data of mature rice in the field [27,28], as shown in Figure 5. (1) Determine neighborhood radius r and the standard deviation-multiple threshold m.
(2) Traverse all point cloud data xi within the neighborhood radius r, and calculate the average distance di between all point cloud data; the calculation of d is as Equation (1) Under massive experimental verification, when the neighborhood radius r = 5 and the standard-deviation multiple threshold m = 2, the noise point removal based on the SOR algorithm results in efficient noise removal for the point cloud data of mature rice in the field [27,28], as shown in Figure 5.

Ground and Crop Point Cloud Classification
Ground points can be defined as the cloud data in which the laser beam penetrates the crop and then reaches the ground surface. Crop points include the point cloud data of the rice stem, leaf, and panicle. The elevation values of the ground cloud data and the adjacent crop point cloud data will have an obvious mutation in some regions, especially at the edge of objects. The discontinuity of the partial elevation values of the contour edge can reflect the shape change of the object. Thus, the moving surface fitting was used to classify the 3D point cloud data of field-mature rice in this study. A quadric surface was used to obtain the partial special surface that was fitted. The moving surface fitting method determines a neighborhood range with the inner interpolation point as the center and determines the fitting function according to the number of all or part of the sampling points falling in the neighborhood range.
The elevation value is expressed according to the fitted value, and the definition of a quadric surface is shown in the following equation:

Ground and Crop Point Cloud Classification
Ground points can be defined as the cloud data in which the laser beam penetrates the crop and then reaches the ground surface. Crop points include the point cloud data of the rice stem, leaf, and panicle. The elevation values of the ground cloud data and the adjacent crop point cloud data will have an obvious mutation in some regions, especially at the edge of objects. The discontinuity of the partial elevation values of the contour edge can reflect the shape change of the object. Thus, the moving surface fitting was used to classify the 3D point cloud data of field-mature rice in this study. A quadric surface was used to obtain the partial special surface that was fitted. The moving surface fitting method determines a neighborhood range with the inner interpolation point as the center and determines the fitting function according to the number of all or part of the sampling points falling in the neighborhood range.
The elevation value is expressed according to the fitted value, and the definition of a quadric surface is shown in the following equation: where a 0 , a 1 , a 2 , a 3 , a 4 , a 5 are the multinomial coefficients to be solved, (x i , y i , z i ) represents the 3D coordinates of points on the surface. It can be seen from the equation that a total of six known points could be calculated through the multinomial coefficient a i . Due to the shelter of crops, the obtained ground point cloud is relatively rare. Through the moving surface method, the holes in the ground area could be repaired to obtain an integral ground area, and the integral ground area can also be extracted to distinguish it from the crop point cloud.
When the surface contains only three known points, the fitted surface is a plane. The plane is expressed as follows: As shown in Figure 6, the specific process is as follows: Agronomy 2022, 12, x FOR PEER REVIEW 8 of 1 Figure 6. Flow chart of moving surface classification.

Mesh Generation Point Cloud Elevation
To ensure the plant height estimation resolution, mesh generation was conducted fo rice point cloud data in each region. Set the side length of the rectangular area to b meshed as L1 × L2. The side length of the single square grid is l after segmentation; th number of grids is a × b.
The grid division process is as follows: (1) Calculate the boundaries of point cloud data (L1, L2) = (Xmax − Xmin, Ymax − Ymin); then the coordinate XOY is established according to the data boundaries.

Estimation of Rice Plant Height Distribution
According to the mesh generation of point cloud data, the average elevation value z of ground point cloud data in each grid was calculated. The equation is as follows:

1.
All point cloud data have 3D coordinates (X i , Y i , Z i ). In order to speed up the process of finding adjacent points and reduce the time spent in subsequent classification, all point cloud data are sorted based on the K-D tree algorithm.

2.
Among all point cloud data, the three points with the smallest height value are taken as the initial ground points. The x and y values of these three points are input into Equation (5), and the surface equation of the initial ground is obtained by fitting; the initial ground height is Z o .

3.
In order to reduce the amount of data processing, a threshold T 1 is set, which is determined according to the fluctuation of the ground to remove the points that do not belong to the ground points, determined according to the following rules:

4.
For the remaining points that need to be further classified (that is, the target area), find the adjacent points of the initial ground point according to the K-D tree and substitute the x and y values of the adjacent points into the initial ground equation obtained in step 2 to obtain the fitting height Z i '. If the difference between the fitted value and the actual value is within the set threshold range, it is considered that the point belongs to the ground point cloud data. Otherwise, it does not belong to the ground point data. 5.
The newly determined ground points are fitted with the initial ground points. When the number of ground points is greater than six, quadratic surface fitting is performed according to Equation (4). Ground points and discriminate new data are added continuously if the adjacent new point cloud data are ground points. Once the number of fitted ground points reaches the set number, the number of ground point cloud data would be retained. As the new ground point is being determined, the old one would be removed, and the surface equation would be updated as well. 6.
Step 5 would be repeated until all point cloud data are determined; then the ground and non-ground point cloud data would be saved separately.

Mesh Generation Point Cloud Elevation
To ensure the plant height estimation resolution, mesh generation was conducted for rice point cloud data in each region. Set the side length of the rectangular area to be meshed as L 1 × L 2 . The side length of the single square grid is l after segmentation; the number of grids is a × b.
The grid division process is as follows: (1) Calculate the boundaries of point cloud data (L 1 , L 2 ) = (X max − X min , Y max − Y min ); then, the coordinate XOY is established according to the data boundaries. (2) Set the grid side length l according to the requirements of plant height estimation resolution, and the number of grids a and b are determined based on Equations (8) and (9).

Estimation of Rice Plant Height Distribution
According to the mesh generation of point cloud data, the average elevation value z ij of ground point cloud data in each grid was calculated. The equation is as follows: where i ∈ [1,2,3 . . . . . . a], j ∈ [1,2,3 . . . . . . b], k is the number of point clouds in the grid in the row i, column j. After calculating the average elevation value z ij of ground point cloud data and the maximum elevation of crop point z ij , the elevation values of all the crop point cloud data in the grid in row i and column j are sorted from the largest to the smallest, and the top ten-point cloud data values in the ranking of elevation values are taken to calculate the average value: where i ∈ where i ∈ The difference between estimated and measured values in each grid was analyzed. The root mean square error (RMSE) and mean absolute percentage error (MAPE) were used to evaluate the estimation accuracy of rice plant height [29].

Results
As shown in Figure 3, the sampling frame was divided into 10 ×10 grids with cotton thread, and the average height in each grid was measured manually. A group of measured data is shown in Table 2.  1  775  700  780  790  795  793  765  790  780  795  2  750  750  760  800  820  800  780  785  790  790  3  790  810  810  810  800  805  800  810  805  810  4  810  815  810  830  845  845  845  830  830  820  5  820  860  860  875  890  890  890  830  820  800  6  830  835  840  820  820  885  895  840  830  815  7  840  845  845  845  850  850  850  850  850  815  8  830  850  835  825  800  810  810  860  760  800  9  800  800  820  810  800  830  825  815  815  790  10  785  770  775  770  800  760  780  765  780  790 Correlation analysis was conducted between the estimated results and the actual results, as shown in Table 3. As can be seen from the table, the measured plant height in the field is similar to the plant height calculated by the method in this study. The average RMSE and MAPE of the estimated and measured values in the field were 8.29% and 9.28% respectively, with an average accuracy of 90%. The average correlation coefficient R 2 of the estimated plant height and the actual results was 0.9, reaching a significant level. Thus, the method proposed in this study could be used to accurately estimate the distribution of height in a field of mature rice. Classification results of ground points and crop points of two groups by using moving surface fitting are shown in Figure 7. Table 1 shows collected images as an example, and the number of ground points in each grid is shown in Figure 8a. There are few ground points in the middle area of the image, which is caused by the occlusion of crops. The average elevation of ground points is shown in Figure 8b. It can be seen from the figure that the elevation of the ground point is about 20 mm. Since there is a certain height difference between the coordinate system used and the ground, the ground point height is not 0 mm. Classification results of ground points and crop points of two groups by using moving surface fitting are shown in Figure 7.  Table 1 shows collected images as an example, and the number of ground points in each grid is shown in Figure 8a. There are few ground points in the middle area of the image, which is caused by the occlusion of crops. The average elevation of ground points is shown in Figure 8b. It can be seen from the figure that the elevation of the ground point is about 20 mm. Since there is a certain height difference between the coordinate system used and the ground, the ground point height is not 0 mm. The number of crop points in each grid is shown in Figure 8c; the distribution of crop points is more even, with an average of more than 200 points per grid. Grids with fewer point clouds are marked in the pictures. The maximum elevation of crop points in each grid is shown in Figure 8d. Due to the small number of crop points in the surrounding The number of crop points in each grid is shown in Figure 8c; the distribution of crop points is more even, with an average of more than 200 points per grid. Grids with fewer point clouds are marked in the pictures. The maximum elevation of crop points in each grid is shown in Figure 8d. Due to the small number of crop points in the surrounding area, the elevation of crop points is low. The surface height of crops can be calculated by the elevation of ground points and crop points. The rice plant height in each grid was estimated for five groups of field-mature rice samples; estimation results of plant height distribution are shown in Figure 9. The number of crop points in each grid is shown in Figure 8c; the distribution of crop points is more even, with an average of more than 200 points per grid. Grids with fewer point clouds are marked in the pictures. The maximum elevation of crop points in each grid is shown in Figure 8d. Due to the small number of crop points in the surrounding area, the elevation of crop points is low. The surface height of crops can be calculated by the elevation of ground points and crop points. The rice plant height in each grid was estimated for five groups of field-mature rice samples; estimation results of plant height distribution are shown in Figure 9.

Discussion
Experimental results showed that the proposed method based on the moving surface and 3D point cloud elevation is applicable to estimate rice height. James et al., [30] used a camera to obtain 3D point cloud data. Through Pix4D's (Pix4D, Co., Ltd., Switzerland) inbuilt measurement tool, the height was measured at four randomly selected points within the top layer of canopy, and the R 2 between the measured height and the actual height was 0.84. The R 2 of the method proposed in this paper was higher, which may result from the accuracy of the reconstructed point cloud by camera surround shooting being lower than that obtained directly by the LiDAR. A lack of noise removal made the reconstructed 3D cloud data retain abnormal points that affected the accuracy. Chen et al., [31] used a combination of LiDAR and IMU to obtain a transformation matrix between a LiDAR coordinate system and ground coordinate system through calibration and then used the DBSCAN adaptive point cloud clustering method to identify the one with the largest average height in the grain area point cloud clustering category as the final grain height; its clustering method was easily affected by stubble and crushed stalks and could not accurately segment the crop area and the harvested area. Besides, the method relies on the transformation matrix between the LiDAR coordinate system and the ground coordinate system being accurate and then to guarantee the point whose elevation value is 0 is correct on the ground surface. The terrain in the farmland was bumpy, and the movement of the harvester may have slight angle changes of pitch and roll. The transformation matrix obtained by its calibration was extremely sensitive to this change, which easily resulted in errors in the plant height measurement. In contrast, the proposed method in this study could be less affected by rough terrain in farmland.
In addition, compared with the previous literature, the method applied to the crops planted sparsely or planted with ridges could not be effective for rice plants that are cultivated densely and have a luxuriant canopy when mature. Obtaining ground and canopy data simultaneously has become the fundamental reason that hinders the further applica-

Discussion
Experimental results showed that the proposed method based on the moving surface and 3D point cloud elevation is applicable to estimate rice height. James et al., [30] used a camera to obtain 3D point cloud data. Through Pix4D's (Pix4D, Co., Ltd., Prilly, Switzerland) inbuilt measurement tool, the height was measured at four randomly selected points within the top layer of canopy, and the R 2 between the measured height and the actual height was 0.84. The R 2 of the method proposed in this paper was higher, which may result from the accuracy of the reconstructed point cloud by camera surround shooting being lower than that obtained directly by the LiDAR. A lack of noise removal made the reconstructed 3D cloud data retain abnormal points that affected the accuracy. Chen et al., [31] used a combination of LiDAR and IMU to obtain a transformation matrix between a LiDAR coordinate system and ground coordinate system through calibration and then used the DBSCAN adaptive point cloud clustering method to identify the one with the largest average height in the grain area point cloud clustering category as the final grain height; its clustering method was easily affected by stubble and crushed stalks and could not accurately segment the crop area and the harvested area. Besides, the method relies on the transformation matrix between the LiDAR coordinate system and the ground coordinate system being accurate and then to guarantee the point whose elevation value is 0 is correct on the ground surface. The terrain in the farmland was bumpy, and the movement of the harvester may have slight angle changes of pitch and roll. The transformation matrix obtained by its calibration was extremely sensitive to this change, which easily resulted in errors in the plant height measurement. In contrast, the proposed method in this study could be less affected by rough terrain in farmland.
In addition, compared with the previous literature, the method applied to the crops planted sparsely or planted with ridges could not be effective for rice plants that are cultivated densely and have a luxuriant canopy when mature. Obtaining ground and canopy data simultaneously has become the fundamental reason that hinders the further application of this detection method to estimate the height of rice crops. Tilly's [18] finding based on CSM and DEM also needs the collection of ground information when the crops are in the young sprout period; in addition, it may be difficult to determine whether the ground profile would experience great change considering the influence of external natural factors and the continuous operation of other plant machines. Using auxiliary markers to provide a reference to calculate crop height is a direct way to solve the problem. Sritarapipat's [20] method effectively simplified the process of plant height detection, but its detection range is narrow, requiring the manual placement of a large number of markers in advance, which is not beneficial in small fields in Asia that need combine harvesters to transfer frequently.
In this study, the average correlation coefficient R 2 reached 0.9. It showed a high accurate, but estimated plant height always seemed to be smaller than the actual value. There may be a small number of point clouds at the root of crops that are misclassified as ground points, resulting in the fitted ground height being higher than the actual ground height, but the estimated accuracy could be an effective reference for further feed-rate calculation. Other experiments are still needed in large fields. For some rice varieties, lodging is a common phenomenon, and the feasibility to apply this method needs further investigation. In addition, it is expected that crop height could be estimated when the harvester is driven at a high speed, so future research is recommended to improve the system that could be fixed on the harvester for real-time estimation.

Conclusions
Our investigation on the height estimation of mature rice plants was based on the moving surface and 3D point cloud elevation by means of a LiDAR phenotyping process. In this study, the ground profile that covered by the luxuriant rice canopy was effectively classified, and the elevation of the ground and canopy was obtained. SOR was used for noise point removal, and useful cloud point data for several samples of mature rice plants in the field were extracted to enhance the estimation accuracy. A comparison of the estimated crop height and actual values indicated the average RMSE and MAPE of the estimated and measured values in the field were 8.29% and 9.28%, respectively. The average accuracy reached 90%, and the average correlation coefficient R 2 of the estimated plant height and the actual results was 0.9, reaching a significant level. The field validation experiments indicated that the moving surface and 3D point cloud elevation could be used in the LiDAR measurement system to estimate mature rice plant height when the ground profile is difficult to identify, which would be a feasible strategy to guide feed rate estimation in advance to improve the harvesting efficiency.

Data Availability Statement:
The data presented in this study are available on request to the corresponding author.