A field wheat stem and tiller number extraction method based on a voxel interpolation mean shift algorithm

Abstract

The application discloses a kind of field wheat stem tiller number extraction method, comprising the following steps:First, the sample of wheat stem tiller number and point cloud data are collected, and the data is preprocessed, in data preprocessing, introduce Kalman filtering algorithm, to remove redundant noise, improve initial point cloud data quality;Second, after data is voxelized, the mathematical relationship between porosity and voxel point cloud density is used to interpolate missing voxels in the canopy, to reduce the influence of occlusion on the algorithm, to obtain relatively complete canopy point cloud information;Finally, the mean shift algorithm is used to cluster the interpolated canopy, and the cluster number is the stem tiller number.The wheat stem tiller number extracted by the method is compared with the field measured stem tiller number, and the feasibility of the algorithm is verified.The application solves the noise and occlusion problem in the application of ground-based laser radar to a certain extent, and provides technical support for the estimation of other growth parameters by ground-based laser radar in the future.

Description

TECHNICAL FIELD

[0001] The present application relates to a field wheat stem tiller number extraction method, in particular to a field wheat stem tiller number extraction method based on a mean shift algorithm of voxel interpolation, and belongs to the field of non-destructive monitoring of crop life information in precision agriculture. BACKGROUND

[0002] The wheat stem tiller number determines the number of effective spikes per mu, and therefore has a close relationship with yield. A reasonable wheat stem tiller number is the basis of high and stable yield. At the same time, the wheat stem tiller number is also an important population structure information of wheat, which reflects the adaptability of wheat to the environment. Therefore, real-time and accurate acquisition of the wheat stem tiller number is helpful for reasonable setting of population density, precise fertilizer and water management, and exploration of crop high-yield and stable-yield technology.

[0003] The traditional manual wheat stem tiller number and spike number counting method is time-consuming and laborious, and is greatly affected by human error, and is no longer suitable for high-throughput wheat phenotype parameter research. The development of remote sensing technology provides a new efficient and non-destructive monitoring means. At present, the research on wheat stem tiller number and spike number has gradually developed from one-dimensional spectrum to two-dimensional image, and from linear model based on vegetation index to deep learning of wheat stem tiller number and spike number characteristics using rich color and texture information in RGB image. However, optical remote sensing methods are generally greatly affected by light and environment, and lack more closely related spatial and volumetric information of plant yield and other related traits.

[0004] As an active remote sensing sensor, laser radar not only overcomes the disadvantage of optical remote sensing affected by light conditions, but also has the advantage of penetrating the canopy to obtain accurate three-dimensional coordinate information of the canopy surface and interior. At present, there is still much room for development in the use of ground-based laser radar to estimate wheat stem tiller number, and the method of estimating wheat stem tiller number based on three-dimensional point cloud has important significance for promoting the application of ground-based laser radar in crop yield phenotype parameter acquisition. SUMMARY

[0005] The purpose of the present application is to provide a field wheat stem tiller number extraction method. In view of the problems of excessive noise and serious occlusion of wheat canopy point cloud, the Kalman filter algorithm is used to remove redundant noise to improve the quality of original point cloud data, the voxel information in the unsampled voxel field is valued to reduce the influence of occlusion, and then the mean shift algorithm is used to cluster the interpolated canopy in combination with the characteristics of point cloud density. Finally, the improvement of the field wheat stem tiller number estimation algorithm is realized, the algorithm accuracy is improved, and the algorithm robustness is enhanced.

[0006] To achieve the above purpose, the technical scheme adopted by the present application is as follows: a field wheat stem tiller number extraction method based on a mean shift algorithm of voxel interpolation, comprising the following steps:

[0007] Step 1, collecting wheat stem tiller number samples and point cloud data;

[0008] Step 2, pre-process the data, and introduce Kalman filtering algorithm to the processed data;

[0009] Step 3, voxelize the point cloud data;

[0010] Step 4, calculate the wheat canopy porosity;

[0011] Step 5, calculate the wheat canopy point cloud density;

[0012] Step 6, calculate the mathematical relationship between the wheat canopy porosity and the point cloud density;

[0013] Step 7, interpolate the missing voxels in the canopy;

[0014] Step 8, cluster the interpolated canopy using the mean shift algorithm.

[0015] Further, in step 1, a three-dimensional laser scanner is used to obtain wheat point cloud, and wheat tiller number samples are obtained synchronously.

[0016] Further, in step 2, the point cloud preprocessing includes point cloud registration and denoising processing, wherein,

[0017] The point cloud registration includes: in the point cloud processing software RiSCAN PRO, using ICP algorithm to register different coordinate systems to the same coordinate system;

[0018] The denoising processing includes: (1) removing drift points: using setting elevation threshold to remove noise points floating in the air; (2) removing mixed points: using setting Deviation value to remove deviation points, wherein the Deviation value is an index for measuring the change between the transmitted pulse waveform and the echo pulse waveform.

[0019] Further, in step 3, the voxelization of point cloud data includes:

[0020] a, establishing initial voxels, calculating the maximum coordinate value point and the minimum value point of all points in the data set;

[0021] b, dividing the initial voxels, dividing the initial voxels into single small voxels along the XYZ axis direction according to a certain step length, and completing the voxelization of the point cloud data surface;

[0022] c, removing invalid voxels, after completing the voxelization of the point cloud surface, using Bresenham algorithm to voxelize the internal space of the point cloud, and removing the blank voxels within the space bounding box outside the plant.

[0023] Further, in step 4, calculating the wheat canopy porosity includes:

[0024] a. The porosity calculation formula is:

[0025] Where Ntotal is the total number of all received / recorded pulses, NoR is the number of returns per pulse, and 1 / NoR is the weight assigned to each loop in this given pulse;

[0026] b. For the three returns recorded by the ground-based laser radar, the point cloud data marks the number of returns for each pulse for each point. Let the number of points marked as 1 be a, the number of points marked as 2 be b, and the number of points marked as 3 be c. Then the number of complete pulses is c, with a weight of 1(3 / 3); the number of two-echo pulses is (b-c), with a weight of 2 / 3; and the number of one-echo pulses is (a-b), with a weight of 1 / 3.

[0027] Substituting the above data into the porosity calculation formula gives:

[0028] c. According to the extinction model, infinitesimal leaves are randomly distributed in space, and it is assumed that the leaf azimuth angle is randomly distributed. Through the Poisson model, P gap is obtained in the direction θ, which is associated with LAI, and the formula is: gap P (θ)=e -G(θ)·LAI / cos(θ) ,

[0029] Where G(θ) is the leaf projection function, which is equal to the projection of a unit area of the plant in the direction perpendicular to the plane, and is uniform in all directions.

[0030] d. Under the condition of canopy aggregation, the Poisson model is modified, and the formula is:

[0031]

[0032] Where Ω is the aggregation index, and k is the extinction coefficient.

[0033] Further, in step 5, the wheat canopy point cloud density calculation formula is:

[0034]

[0035] Where D represents the point cloud density, P total represents the total number of points, x max , y max , and z max represent the maximum values of the point cloud data in the x, y, and z directions, respectively, and x min , y min , and z min represent the minimum values of the point cloud data in the x, y, and z directions, respectively.

[0036] Further, in step 6, the wheat canopy is divided into height layers of different heights H (H = 1, 2, 3, 4 cm). The relationship between the overall porosity of the canopy at different heights and the point cloud density corresponding to the layer height is described by the same form of mathematical expression, which is:

[0037]

[0038] where a1, b1, c1, a2, b2, and c2 are generated by Gaussian curve fitting according to the input values. When H = 3 cm, the fitting effect is the best, so 3 cm is set as the layer height in subsequent interpolation.

[0039] Further, in step 7, the interpolation of the missing voxels in the canopy includes:

[0040] a. Fitting expression: Calculate the density D and porosity Pgap of the point cloud of each 3 cm layer, and obtain the expression f(x) according to Gaussian curve fitting.

[0041] b. Calculate the porosity of the neighborhood of the unsampled voxel i. The voxel located inside the canopy but not containing data points is an unsampled voxel. Retrieve eight voxels in the neighborhood of any unsampled voxel and calculate the porosity of the local canopy formed by these eight voxels.

[0042] c. Calculate the point cloud density of the unsampled voxel i. Substitute the porosity calculated in (b) into the expression f(x) obtained in (a) to obtain the point cloud density value Di of the unsampled voxel.

[0043] d. Repeat steps (a), (b), and (c) until each unsampled voxel is assigned a value, complete the interpolation of the unsampled voxels, and obtain the complete canopy information. The final voxel is represented by , where G represents the gravity point of the voxel, and D represents the point cloud density value of the voxel.

[0044] Further, in step 8, the interpolated canopy is clustered using the mean shift algorithm, including:

[0045] a. Randomly select the gravity point of any voxel in the unmarked data points as the center point.

[0046] b. Find all voxels within the bandwidth distance from the center point and classify these voxels into cluster 1.

[0047] c. Compare the point cloud density from the center point to the gravity point of each voxel in the set, and move the center point to the center point with the maximum density. The moving vector is the shift vector shift.

[0048] ​​d. moving the center point along the shift direction of the shift vector ||shift||;

[0049] f. repeating steps b-d until convergence is reached, and recording the center point position at this time; all points traversed during the iteration process are classified as a cluster;

[0050] g. when convergence is reached, if the distance between the center point of the current cluster and the center point of other existing clusters is less than the threshold value, i.e. bandwidth / 2, the clusters are merged, otherwise, the current cluster is taken as a new cluster, and the number of categories is increased;

[0051] h. repeating the above steps a-g until all voxels are labeled and accessed;

[0052] i. according to each category, taking the category with the maximum point access frequency as the category to which the current voxel set belongs;

[0053] j. counting the number of categories, i.e. the final number of tillers.

[0054] Further, the feasibility of the extraction method is verified, including: verifying the feasibility and accuracy of the test algorithm by using independent year tiller sample data, and comprehensively evaluating the results by using the determination coefficient R2 and the root mean square error RMSE;

[0055] R 2 and the formula of RMSE is:

[0056] Wherein, Q i is the estimated value, Pi is the true value, the average value of the true value, and n represents the total number of samples.

[0057] Compared with the prior art, the beneficial effects of the present application are:

[0058] The present application aims at the noise problem, and introduces a Kalman filter algorithm in the point cloud preprocessing step to further remove redundant noise points. For the problem of missing point cloud data caused by the shielding effect, a mathematical relationship between the voxel point cloud density and the canopy porosity is established to try to interpolate the missing voxels. The field wheat tiller number extracted by the method of the present application is compared with the field measured tiller number, and the verification result shows that the determination coefficient R2 is 0.69, and the root mean square error RMSE is 79 / m2, which proves the feasibility of the method of the present application. BRIEF DESCRIPTION OF DRAWINGS

[0059] Figure 1 is the field wheat tiller number extraction method flow chart of the voxel interpolation based mean shift algorithm of the present application.

[0060] Figure 2 is a schematic diagram of introducing Kalman filtering for wheat point cloud preprocessing; wherein,Figure 2 a is the point cloud data before introducing Kalman filter, Figure 2 b is the point cloud data after introducing Kalman filter processing.

[0061] Figure 3 is the point cloud data voxelization process; wherein, Figure 3 a is the original point cloud, Figure 3 b is the initial voxel, Figure 3 c is the wheat point cloud after voxelization.

[0062] Figure 4 is the relationship between the overall canopy porosity of the test plot and the point cloud density.

[0063] Figure 5 is the relationship between the layered canopy porosity and the point cloud density; Figure 5 a is the layer height H equal to 1cm, Figure 5 b is H equal to 2cm, Figure 5 c is H equal to 3cm, Figure 5 d is H equal to 4cm.

[0064] Figure 6 is the mean shift algorithm clustering process. DETAILED DESCRIPTION

[0065] The present application will be described in detail below in conjunction with the drawings and specific embodiments.

[0066] As shown in Figure 1 , the present application provides a mean shift algorithm based on voxel interpolation method for extracting field wheat stem and tiller number, comprising the following steps:

[0067] S101: obtaining original point cloud data;

[0068] S102: point cloud data preprocessing;

[0069] S103: point cloud data voxelization;

[0070] S104: calculating the wheat canopy porosity;

[0071] S105: calculating the wheat canopy point cloud density;

[0072] S106: interpolating the missing voxels in the canopy;

[0073] S107: using the mean shift algorithm to cluster the interpolated canopy.

[0074] The steps will be described in detail as follows:

[0075] S101: Obtain field wheat point cloud by RIEGL VZ-1000 three-dimensional laser scanner, and synchronously obtain field wheat tiller number measured data. The field wheat test is selected as follows:

[0076] 2 wheat varieties: Shengxuan No. 6 (V1), Yangmai 16 (V2).

[0077] 2 density levels: D1 row spacing is 25 cm (2.4 x 10 6 seedlings / ha); D2 row spacing is 40 cm (1.5 x 10 6 seedlings / ha).

[0078] 3 nitrogen levels: pure nitrogen 0 kg / ha (N0), 150 kg / ha (N1), 300 kg / ha (N2); nitrogen, phosphorus and potassium are urea, superphosphate and potassium chloride, respectively, of which 50% of nitrogen is base fertilizer, and 50% is applied at the jointing stage, and the others are all applied as base fertilizer.

[0079] Wheat planting method: strip seeding. The test adopts random block design with 3 repetitions, a total of 36 plots. The area of each plot is 30 m 2 (6 m x 5 m), and the total area is about 1080 m 2 .

[0080] Wheat point cloud acquisition method: 8-station laser radar test scheme is adopted, and the scanning mode is 60 mode, i.e. the angular resolution is 0.06°.

[0081] Method for obtaining field wheat tiller number measured data: select 2 rows x 1 m of wheat test area in each plot, and test the artificial tiller number, record the tiller number of one meter double row wheat in each plot, take the average value, i.e. the tiller number of one meter single row wheat m, then the whole row wheat measured tiller number is 6*m. The measured tiller number of the whole plot is 6*m*the number of rows in the plot (row spacing is 25 cm, then there are 26 rows; row spacing is 40 cm, then there are 17 rows), to obtain the wheat tiller number of the whole plot. The specific conditions of field test are as shown in Table 1.

[0082] Table 1

[0083]

[0084] S102: Data registration and denoising processing are performed on the wheat point cloud obtained in step S101.

[0085] Point cloud registration: Each scanning station will have its own coordinate system-scanner coordinate system. In the point cloud processing software RiSCAN PRO, the ICP algorithm is used to register different coordinate systems to the same coordinate system (Project coordinate system). The ICP algorithm is an optimal matching method based on the least square method, repeatedly determines the nearest point pair, and calculates the coordinate conversion parameters until the given convergence accuracy is met, and the iteration ends. For two point cloud data sets A and B with some overlapping parts, A is the reference point cloud data set, and point cloud B is the coordinate system to be registered. It needs to be converted to the coordinate system of point cloud A through rotation and translation.

[0086] Denoising processing: (1) Remove floating points: Set the elevation threshold to remove noise points floating in the air. After elevation screening, abnormal point cloud data can be clearly seen and deleted. In this embodiment, the elevation threshold is set to 2m. (2) Remove mixed points: Set the Deviation value to remove deviation points. The Deviation value is an index for measuring the change between the transmitted pulse waveform and the echo pulse waveform. When the Deviation value is 0, it means that the waveform is not distorted, and the received pulse shape is completely consistent with the transmitted pulse shape; when the Deviation value is very large, it indicates that the waveform is severely distorted, and these point cloud data need to be deleted. In this embodiment, the Deviation is set to 250.

[0087] The point cloud data after point cloud registration and denoising processing is introduced into the Kalman filtering algorithm.

[0088] The calculation process of Kalman filtering is as follows:

[0089] a. Calculate the state estimate value:

[0090] b. Calculate the one-step prediction of the state:

[0091] c. Calculate the information:

[0092] d. Calculate the Kalman filter gain: K(t+1) = P(t+1|t)H T [HP(t+1|t)H T +R] -1

[0093] e. Calculate the next step prediction mean square error: P(t+1|t) = ΦP(t|t)Φ T +ΓQΓ T

[0094] f, calculate the next step prediction estimate mean square error: P(t+1|t+1)=[I n -K(t+1)H]P(t+1|t)

[0095] For Kalman filter, only need to give the initial value And P0, according to the observation value Zk at k time, the state estimation at k time can be obtained According to the variance estimation between the estimated value and the observation value, the estimated value is selected.

[0096] S103: voxelizing the point cloud data after S102 step.

[0097] a, establish an initial voxel, calculate the maximum coordinate value point Pmax(Xmax, Ymax, Zmax) and the minimum value point Pmin(Xmin, Ymin, Zmin) of all points in the data set. The initial voxel (space bounding box) containing all data points is calculated, and the length, width and height are a=Xmax-Xmin, b=Ymax-Ymin and c=Zmax-Zmin.

[0098] b, divide the initial voxel, and divide the initial voxel into single small voxel along the XYZ axis direction according to a certain step length, and complete the voxelization of the point cloud data surface. The step length is called voxel resolution, and the size is generally determined according to the original data point cloud spacing, and different sizes of voxel resolution can be set in different directions. The voxel resolution of XYZ three directions of the present application is consistent, and the specific size is selected by trial according to the specific situation, and here, in order to adapt to the subsequent voxel interpolation step, the height is adjusted, and 1cm voxel resolution is adopted.

[0099] c, remove invalid voxel, after completing the voxelization of the point cloud surface, the Bresenham algorithm is used to voxelize the internal space of the point cloud, and the blank voxel inside the space bounding box of the plant is removed, so as to reduce unnecessary search and calculation. Bresenham algorithm is the most widely used straight line scanning conversion method in computer graphics field, which is extended to three-dimensional space, and P min is the starting point and each data point is the end point, and the voxel space of the data point is determined.

[0100] S104: calculate the porosity of the point cloud data processed in S102 step.

[0101] a, the porosity calculation formula is:

[0102] Where Ntotal is the total number of all received / recorded pulses, and NoR is the number of returns per pulse. 1 / NoR is the weight assigned to each loop in this given pulse.

[0103] b、According to the extinction model, the infinitesimal leaf blade is randomly distributed in space, and the leaf azimuth angle is randomly distributed. By the Poisson model, P gap associated with LAI. The formula is: P gap (θ) = e -G(θ)·LAI / cos(θ) .

[0104] Where G(θ) is the leaf projection function, which is equal to the projection of unit area of the plant in the direction perpendicular to the plane, and is uniform in all directions.

[0105] c、Under the condition of canopy aggregation, the Poisson model can be modified. The formula is:

[0106] Where Ω is the aggregation index, and k is the extinction coefficient.

[0107] S105: Calculate the point cloud density of the point cloud data processed in step S102.

[0108] The calculation formula is:

[0109] Where D represents the point cloud density, P total represents the number of all points, x max , y max , and z max represent the maximum values of the point cloud data in the x, y, and z directions, respectively, and x min , y min , and z min represent the minimum values of the point cloud data in the x, y, and z directions, respectively.

[0110] S106: Interpolate the missing voxels in the canopy in step S103.

[0111] The wheat canopy is divided into height layers of different heights H (H = 1, 2, 3, 4 cm). Under the division of different height layers, the relationship between the overall porosity of the canopy at different heights and the point cloud density corresponding to the layer height can be described by a mathematical expression of the same form.

[0112] The formula is:

[0113] Where a1, b1, c1, a2, b2, and c2 are generated by Gaussian curve fitting according to the input values. When H = 3 cm, the fitting effect is the best (R 2 highest). Therefore, 3 cm is set as the layer height in subsequent interpolation.

[0114] The specific steps for interpolating the missing voxels in the canopy are as follows:

[0115] a. Fitting expression. Calculate the density D and the porosity P gap of each 3cm layer point cloud, and fit the expression f(x) according to Gaussian curve.

[0116] b. Calculate the porosity of the neighborhood of the unsampled voxel i. The voxels that are located in the canopy but do not contain data points are called unsampled voxels. Retrieve the eight voxels within the arbitrary neighborhood of the unsampled voxel, and calculate the porosity of the local canopy composed of these eight voxels

[0117] c. Calculate the point cloud density of the unsampled voxel i. Substitute the porosity calculated in (b) into the expression f(x) obtained in (a), and the point cloud density value Di of this unsampled voxel can be obtained, that is, the neighborhood information is assigned to it.

[0118] d. Repeat steps (a), (b), and (c) until each unsampled voxel is assigned, complete the interpolation of the unsampled voxels, and obtain the more complete canopy information. The final voxel can be represented as , , where the center of the voxel is represented, and D represents the point cloud density value of the voxel.

[0119] S107: Cluster the canopy after S106 using the mean shift algorithm.

[0120] a. Randomly select the center point of any voxel in the unmarked data points as the center point;

[0121] b. Find all voxels within the bandwidth distance from the center point, and classify these voxels into a cluster. The bandwidth is the only parameter that needs to be defined in this algorithm, and the specific value can be selected according to the results obtained from different values, and the bandwidth value corresponding to the value closest to the reference value is selected. The bandwidth value used in this method is between 0.5 and 1;

[0122] c. Compare the point cloud density from the center point and the center point of each voxel in the set, and move the center point to the center point with the maximum density. The moving vector is the shift vector shift;

[0123] d. Move the center point along the shift vector shift direction ||shift||;

[0124] f. Repeat steps b-d until the iteration converges, and record the position of the center point at this time. The points traversed in the iteration process are classified into clusters;

[0125] g. When the convergence is reached, if the distance between the center point of the current cluster and the center point of other existing clusters is less than the threshold (bandwidth / 2), the clusters are merged. Otherwise, the current is taken as a new cluster, and the number of categories is increased;

[0126] ​h、repeat the above steps until all the voxels are labeled access;

[0127] i、according to each class, take the class with the maximum point access frequency as the current voxel set belonging class;

[0128] j、statistics class number, namely the final tiller number.

[0129] After the above steps are executed, feasibility verification can be performed, including:

[0130] Adopting the coefficient of determination R 2 , root mean square error RMSE and relative root mean square error RRMSE to comprehensively evaluate the results;

[0131] R 2 and RMSE formula is:

[0132] Wherein Q i is the estimated value, P i is the true value, Indicate the average value of true value, n indicates the total number of samples.

[0133] The method constructed can effectively calculate the field wheat tiller number precision, and the overall precision reaches R2=0.69, RMSE=79 tills / m 2 .

[0134] The above shows and describes the basic principle, main features and advantages of the present application. Those skilled in the art should understand that the above examples do not limit the protection scope of the present application in any form, and any technical solutions obtained by equivalent replacement or the like fall within the protection scope of the present application.

[0135] The part not involved in the present application is the same as or can be realized by using the prior art.

Claims

1. A method for extracting the number of stems and tillers of wheat in the field based on a voxel interpolation-based mean shift algorithm, characterized in that The method comprises the following steps: Step 1, collecting wheat stem tiller number samples and point cloud data; Step 2, preprocessing the data, and introducing Kalman filtering algorithm to the processed data; Step 3, voxelizing the point cloud data; Step 4, calculating the wheat canopy porosity; Step 5, calculating the wheat canopy point cloud density; Step 6, calculating the mathematical relationship between the wheat canopy porosity and the point cloud density; in this step, the wheat canopy is divided into height layers of different heights H, H = 1, 2, 3, 4 cm, and the relationship between the overall porosity of the canopy at different heights and the point cloud density corresponding to the layer height is described by a mathematical expression of the same form, the formula being: ， Wherein, a1, b1, c1, a2, b2, c2 are generated according to the input values by Gaussian curve fitting, and when H = 3 cm, the fitting effect is the best, and 3 cm is set as the layer height in subsequent interpolation; Step 7, interpolating the missing voxels in the canopy; in this step, interpolating the missing voxels in the canopy includes: a, fitting expression; respectively calculating the density D and porosity Pgap of the point cloud of every 3 cm layer, and obtaining the expression f(x) according to Gaussian curve fitting; b. Calculate the porosity of the neighborhood of the unsampled voxel i. Voxel inside the canopy but not containing data points is unsampled voxel. Retrieve eight voxels in the neighborhood of any unsampled voxel, and calculate the porosity of the local canopy formed by these eight voxels ; c. Calculate the point cloud density of the unsampled voxel i; Substitute the expression f(x) obtained in a into the porosity calculated in b, and the point cloud density value of this unsampled voxel is obtained c. Calculate the point cloud density of the unsampled voxel i; Substitute the expression f(x) obtained in a into the porosity calculated in b, and the point cloud density value of this unsampled voxel is obtained c. Calculate the point cloud density of the unsampled voxel i; Substitute the expression f(x) obtained in a into the porosity calculated in b d, repeating steps a, b, c until each unsampled voxel is assigned a value, completing the interpolation of the unsampled voxels, obtaining a more complete crown layer information; the final voxel is represented by , represents the voxel center point, and D represents the voxel point cloud density value; Step 8, clustering the interpolated canopy by using the mean shift algorithm.

2. The method according to claim 1, wherein, In the step 1, a three-dimensional laser scanner is used to obtain the wheat point cloud, and the wheat stem tiller number samples are obtained synchronously.

3. The method of claim 1, wherein the method is characterized by, In the step 2, the point cloud preprocessing is point cloud registration and denoising processing, wherein, The point cloud registration includes: in the point cloud processing software RiSCAN PRO, the ICP algorithm is used to register different coordinate systems to the same coordinate system; The denoising processing includes: (1) removing drift points: using a height threshold to remove noise points floating in the air; (2) removing mixed points: using a Deviation value to remove deviation points, wherein the Deviation value is an index for measuring the change between the transmitted pulse waveform and the echo pulse waveform.

4. The method of claim 1, wherein the method is characterized by, In the step 3, the voxelization of the point cloud data includes: a, establishing an initial voxel, calculating the maximum coordinate value point and the minimum value point of all points in the data set; b, dividing the initial voxel, dividing the initial voxel into single small voxels along the XYZ axis direction according to a certain step, and completing the voxelization of the point cloud data surface; c, removing invalid voxels, after completing the voxelization of the point cloud surface, using the Bresenham algorithm to voxelize the internal space of the point cloud, and removing the blank voxels within the bounding box outside the plant.

5. The method of claim 1, wherein the method is characterized by, In the step 4, calculating the wheat canopy porosity includes: a The porosity calculation formula is: ​ wherein, is the total number of all received / recorded pulses, NoR is the number of returns per pulse, and 1 / NoR is the weight assigned to each loop in this given pulse; b, for the three returns recorded by the ground-based laser radar, the return number of pulses for each point is marked, assuming that the number of points marked as 1 is a, the number of points marked as 2 is b, and the number of points marked as 3 is c, then the number of complete pulses is c, the weight is 1; the number of two-way echo pulses is (b-c), the weight is 2 / 3; the number of one-way echo pulses is (a-b), the weight is 1 / 3; Substituting the above data into the calculation formula of porosity gives: ; c. According to the extinction model, the infinitesimal leaves are randomly distributed in space, and the leaf azimuth angle is assumed to be randomly distributed. The leaf area index (LAI) is defined as the total leaf area per unit ground surface area. The LAI is related to the extinction coefficient (K) by the formula: LAI = 1 / K , Wherein, G(θ) is a leaf projection function, equal to the projection of a plant per unit area in the direction perpendicular to the plane, uniform in all directions; d. Under the condition of canopy aggregation, the Poisson model is modified, and the formula is: ， wherein, Ω is the aggregation index, and k is the extinction coefficient.

6. The method of claim 1, wherein the method is characterized by, In step 5, the formula for calculating the point cloud density of the wheat canopy is: ， wherein D denotes the point cloud density, denotes the total number of points, denotes the maximum value in the x, y, z direction of the point cloud data, respectively, denotes the minimum value in the x, y, z direction of the point cloud data, respectively.

7. The method of claim 1, wherein the method is characterized by, In step 8, the mean shift algorithm is used to cluster the interpolated canopy, including: a. Randomly select the center point of any voxel in the unmarked data points as the center point; b. Find all voxels within the bandwidth distance from the center point, and classify these voxels into cluster 1; c. Compare the point cloud density from the center point with the center point of each voxel in the set, move the center point to the center point with the maximum density, and the moving vector is the shift vector shift; d. Move the center point along the shift vector shift direction ||shift||; f. Repeat steps b-d until iteration converges, and record the center point position at this time; the points traversed in the iteration process are classified as clusters; g. When convergence is reached, the distance between the center point of the current cluster and the center point of other existing clusters is less than the threshold, i.e. bandwidth / 2, the clusters are merged, otherwise, the current cluster is taken as a new cluster, and the class number is increased; h. Repeat the above steps a-g until all voxels are marked and accessed; i. According to each class, take the class with the maximum point access frequency as the class to which the current voxel set belongs; j. Count the number of classes, which is the final tiller number.

8. The method of claim 1, wherein the method is a method of extracting wheat stem and tiller number in field based on a voxel interpolation mean shift algorithm, characterized in that, The feasibility verification of the extraction method comprises: verifying the feasibility and accuracy of the test algorithm by using the sample data of stem tiller number in independent years, and comprehensively evaluating the results by using the determination coefficient R 2 and the root mean square error RMSE. R 2 and the formula for RMSE is: , wherein is an estimated value, is a true value, denotes the average of the true values, n denotes the total number of samples.

Images