CN110363861B - Laser radar point cloud-based field crop three-dimensional reconstruction method - Google Patents
Laser radar point cloud-based field crop three-dimensional reconstruction method Download PDFInfo
- Publication number
- CN110363861B CN110363861B CN201910632756.0A CN201910632756A CN110363861B CN 110363861 B CN110363861 B CN 110363861B CN 201910632756 A CN201910632756 A CN 201910632756A CN 110363861 B CN110363861 B CN 110363861B
- Authority
- CN
- China
- Prior art keywords
- point cloud
- point
- crop
- model
- clusters
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
- G01B11/00—Measuring arrangements characterised by the use of optical techniques
- G01B11/24—Measuring arrangements characterised by the use of optical techniques for measuring contours or curvatures
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S17/00—Systems using the reflection or reradiation of electromagnetic waves other than radio waves, e.g. lidar systems
- G01S17/88—Lidar systems specially adapted for specific applications
- G01S17/89—Lidar systems specially adapted for specific applications for mapping or imaging
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/60—Analysis of geometric attributes
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2200/00—Indexing scheme for image data processing or generation, in general
- G06T2200/08—Indexing scheme for image data processing or generation, in general involving all processing steps from image acquisition to 3D model generation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10028—Range image; Depth image; 3D point clouds
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
- G06T2207/10044—Radar image
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30181—Earth observation
- G06T2207/30188—Vegetation; Agriculture
Abstract
The invention discloses a field crop three-dimensional reconstruction method based on laser radar point cloud, which comprises the steps of obtaining point cloud data of field crops through a laser radar, and carrying out data preprocessing on the point cloud data; extracting a small amount of point clouds to serve as key control points, and establishing a digital surface model and a digital ground model by utilizing the key control points and an irregular triangular network; constructing an initial spike layer space model according to the digital surface model, and acquiring point cloud data in the actual spike layer space model by an Otsu threshold method; partitioning the point cloud in the actual spike-layer space model into a plurality of point cloud clusters by adopting a density-based clustering algorithm; correcting the number of the point cloud clusters, wherein the number of the point cloud clusters is the number of the crop ears, and calculating the height of the crop according to the digital ground model; and constructing a three-dimensional model of the crops in the field according to the obtained crop quantity, the crop height and the crop ear coordinates. The method realizes the construction of a three-dimensional model of field crops and provides enough data support for the digital management of the field.
Description
Technical Field
The invention belongs to the field of agricultural machinery intellectualization, and particularly relates to a field crop three-dimensional reconstruction method based on laser radar point cloud.
Background
The three-dimensional construction of the crop model has important research significance for crop plant type analysis and plant canopy physiological and ecological index calculation, and provides theoretical basis for crop digital management, accurate fertilization, pesticide application, yield estimation and the like. The current research on the reconstruction of the three-dimensional model of the crop mainly focuses on the geometrical modeling and visualization research of the overground organ of a single crop, including the morphological construction of stems, leaves and ears of the crop or the three-dimensional morphological modeling of a root system. However, the three-dimensional model of the single plant cannot provide enough data basis for field management.
Disclosure of Invention
The field crop three-dimensional reconstruction method based on the laser radar point cloud comprises the steps of scanning field crops by using an unmanned aerial vehicle carrying the laser radar, calculating the number of the crops, the height of the crops and the coordinates of crop ears by using point cloud data, realizing the construction of a three-dimensional model of the crops, and providing data support for digital management of the field.
In order to achieve the technical purpose, the technical scheme adopted by the invention is as follows:
a field crop three-dimensional reconstruction method based on laser radar point cloud comprises the following steps:
(1) Scanning field crops through a laser radar to obtain point cloud data, and performing data preprocessing on the point cloud data output by the laser radar;
(2) Extracting a small amount of point clouds from the preprocessed point cloud data to serve as key control points, and establishing a digital surface model s and a digital ground model by using the key control points and an irregular triangular network;
(3) Presetting a height value h to be expanded above the digital surface model s u And the height value h to be expanded below d S + h u And s-h d The upper boundary and the lower boundary of the initial spike layer space model are formed together, and point cloud data which do not belong to the initial spike layer space model are filtered from the preprocessed point cloud data; calculating the optimal height threshold h by adopting an Otsu threshold method t Removing the part in the space model of the initial spike layer which is lower than the height threshold h t Obtaining the point cloud in the actual spike layer space modelData;
(4) Dividing the point cloud data in the actual spike layer space model into a plurality of point cloud clusters by adopting a density-based clustering algorithm DBSCAN method, and obtaining the number of the point cloud clusters and each point cloud data in the point cloud clusters;
(5) Correcting the number of the point cloud clusters, taking the corrected number of the point cloud clusters as the number of crop ears, and taking the average value of each point cloud coordinate in each point cloud cluster as the coordinate of each crop ear, wherein the number of the crop ears is the number of crops;
(6) Bringing the X-axis coordinate value and the Y-axis coordinate value of the crop ears into the digital ground model to obtain corresponding ground height values, and subtracting the corresponding ground height values from the Z-axis coordinate values of the crop ears to obtain the actual heights of the crops corresponding to the crop ears;
(7) And drawing a three-dimensional view according to the obtained number of the crops, the actual height of each crop and the coordinates of the crop ears at the top of each crop to obtain a three-dimensional model of the crops in the field.
As a further improved technical scheme of the invention, the step (1) comprises the following steps:
unmanned aerial vehicle carries on laser radar, establishes laser radar coordinate system O respectively L -X L Y L Z L And unmanned aerial vehicle coordinate system O V -X V Y V Z V ；
The unmanned aerial vehicle flies through the field and the laser radar on the unmanned aerial vehicle scans field crops to obtain point cloud data of the field crops, and the point cloud data output by the laser radar is converted from a polar coordinate to a laser radar Cartesian coordinate system O L -X L Y L Z L In the laser radar coordinate system O L -X L Y L Z L Point cloud data in the unmanned aerial vehicle is converted into an unmanned aerial vehicle coordinate system O V -X V Y V Z V In, the unmanned plane coordinate system O V -X V Y V Z V Converting the point cloud data in the space into a world coordinate system O-XYZ;
setting an intensity threshold, filtering point cloud data with the reflection intensity smaller than the intensity threshold in the point cloud data, and removing noise points in the point cloud data by adopting a KD-Tree point cloud denoising method.
As a further improved technical scheme of the invention, the step (2) comprises the following steps:
dividing an XOY horizontal plane of a coordinate system O-XYZ where the preprocessed point cloud data are located into a plurality of cells with consistent sizes, establishing a plurality of cuboids by taking each cell as a bottom surface, respectively searching point clouds with maximum elevation values in each cuboid, taking the found point clouds as first key control points, mapping all the first key control points onto the XOY plane, and establishing a digital surface model s based on the first key control points on the XOY plane by adopting a Delaunay irregular triangulation network generation method;
dividing an XOY horizontal plane of a coordinate system O-XYZ where preprocessed point cloud data are located into a plurality of cells with the same size, establishing a plurality of cuboids by taking each cell as a bottom surface, respectively searching a point cloud with the minimum height value in each histogram, taking the found point cloud as a second key control point, mapping all the second key control points onto the XOY plane, and establishing a digital ground model based on the second key control points on the XOY plane by adopting a Delaunay irregular triangulation network generation method.
As a further improved technical solution of the present invention, the establishing of the digital surface model s based on the first key control point on the XOY plane and by using the Delaunay irregular triangulation network generation method specifically includes:
(a) Putting all the first key control points mapped on the XOY plane into a point set, and dividing the point set into two sub-point sets with equal points;
(b) Constructing a Delaunay triangulation network in the sub-point set;
(c) Respectively connecting the bottom line and the top line of a Delaunay triangulation network formed by the two sub-point sets;
(d) Merging the Delaunay triangulation networks of the two sub-point sets from the bottom line to the top line, and optimizing the merged triangulation networks by utilizing a Lawson local optimization algorithm to obtain the triangulation network according with the Delaunay rule;
(e) And (d) mapping the triangulation network finally generated in the step (d) into an O-XYZ three-dimensional space coordinate system to obtain a digital surface model s.
As a further improved technical solution of the present invention, the establishment of the digital ground model based on the second key control point on the XOY plane and by using the Delaunay triangulation method specifically comprises:
(a) Putting all second key control points mapped on the XOY plane into a point set, and dividing the point set into two sub-point sets with equal points;
(b) Constructing a Delaunay triangulation network in the sub-point set;
(c) Respectively connecting the bottom line and the top line of a Delaunay triangulation network formed by the two sub-point sets;
(d) Merging the Delaunay triangulation networks of the two sub-point sets from the bottom line to the top line, and optimizing the merged triangulation networks by utilizing a Lawson local optimization algorithm to obtain the triangulation network according with the Delaunay rule;
(e) And (d) mapping the triangulation network finally generated in the step (d) to an O-XYZ three-dimensional space coordinate system to obtain the digital ground model.
As a further improved technical scheme of the invention, the step (5) specifically comprises the following steps:
(a) Projecting each point cloud data in the point cloud clusters on a horizontal plane, acquiring convex closure bags projected by each point cloud cluster in the horizontal plane, if the convex closure bags of a certain two point cloud clusters are overlapped, setting the convex closure bags of the two point cloud clusters as A and B respectively, and calculating the overlapped areas of the convex closure bags and the B;
(b) Setting a coincidence area threshold T overlap If the overlapping area of the convex closure of some two point cloud clusters is equal to max { S } A ,S B The ratio of } is greater than or equal to T overlap Respectively calculating the average value of the Z-axis coordinates of each point cloud in the two point cloud clusters, and removing the point cloud clusters with low average value of the Z-axis coordinates, wherein S A ，S B Areas of A, B, respectively; if the overlapping area of the convex closure bags of two point cloud clusters is equal to max { S } A ,S B The ratio of is less than T overlap Then continuing to reserve the two point cloud clusters;
(c) The removed point cloud clusters are not used as crop ears, the rest point cloud clusters are used as crop ears, and the corrected number of the point cloud clusters is the number of the crop ears;
(d) And respectively taking the average value of the X-axis coordinates, the average value of the Y-axis coordinates and the average value of the Z-axis coordinates of each point cloud in the point cloud cluster as the coordinate values of the X-axis, the Y-axis and the Z of the corresponding crop ear.
The invention has the beneficial effects that:
the invention directly identifies the ears of the crops through the three-dimensional point cloud data, can obtain the coordinates of the ears of the crops, the actual heights and the number of the crops, has greatly lower labor cost compared with a manual measurement method, can improve the accuracy of the estimation of the crop parameters compared with a method for measuring the information of single crop and counting the crop information of a field, and provides enough data support for the digital management of the field.
Drawings
Fig. 1 is a flowchart of the operation of the present embodiment.
Fig. 2 is a schematic diagram of a relationship between polar coordinates and cartesian coordinates of the laser radar point cloud data according to this embodiment.
Fig. 3 is a schematic diagram of the irregular triangulation network of the present embodiment.
Fig. 4 is a point cloud distribution histogram of the present embodiment.
Fig. 5 is a diagram of the DBSCAN point cloud clustering process after KD-Tree division in this embodiment.
Detailed Description
The following further describes embodiments of the present invention with reference to fig. 1 to 5:
the field crop of the embodiment takes the wheat field as an object, and the field wheat three-dimensional model reconstruction method is shown in fig. 1, and specifically comprises the following steps:
acquiring data and preprocessing the data:
and controlling the unmanned aerial vehicle carrying the laser radar to fly through the wheat field, scanning the wheat by the laser radar to obtain point cloud data, and then preprocessing the data.
The data preprocessing process of the embodiment is as follows:
(1) Point cloud coordinate conversion:
(a) Establishing a laser radar coordinate system O L -X L Y L Z L And unmanned aerial vehicle coordinate system O V -X V Y V Z V . Wherein the radar coordinate system is indicated by the subscript L, origin of coordinates O L Arranged in the radar interior receiver at the centre, X L In front of the radar, Y L And X L Vertical and oriented to the left, Z L Perpendicular to X L O L Y L The plane is upward. The coordinate system of the unmanned aerial vehicle is indicated by subscript V, and the origin of coordinates O V Coincident with the unmanned aerial vehicle centroid, when the unmanned aerial vehicle is in a stationary state on the horizontal plane, X V Pointing forwards parallel to the ground, Y V Pointing to the left, Z V Perpendicular to X V O V Y V The plane is upward. The laser radar output data is in polar coordinate form, and as shown in fig. 2, the polar coordinates are converted into cartesian coordinates by equation (1):
wherein R is from the detection Point Data Point to the origin of coordinates O L The angles ω and α are shown in FIG. 2, ω is the detection point and the origin of coordinates O L Connecting line to it at X L O L Y L The included angle of plane projection, alpha is the detection point and the origin of coordinates O L Is connected to X L O L Y L Planar projection and Y L The angle of the axes.
(b) Laser radar coordinate system O L -X L Y L Z L With unmanned aerial vehicle coordinate system O V -X V Y V Z V The laser radar coordinate system is not overlapped, a certain rotation angle exists between each axis of the radar coordinate system and each axis of the unmanned aerial vehicle coordinate system, the origin of coordinates is not overlapped, and the laser radar coordinate system needs to be converted into the unmanned aerial vehicle coordinate system. According to the installation position of the laser radar on the unmanned aerial vehicle, converting any point P under a laser radar coordinate system into the unmanned aerial vehicle coordinate system by using the formula (2):
in the formula, Q is a rotation matrix, T is a translation matrix, and R and T are obtained by a calibration method.
(c) According to the flight path and attitude information of the unmanned aerial vehicle, the coordinate system O of the unmanned aerial vehicle is determined V -X V Y V Z V The point cloud data in the space is converted into a world coordinate system which is expressed by O-XYZ.
(2) Removing low-intensity points and noise points:
in order to filter point cloud with too low intensity, an intensity threshold value is set, and point cloud data with reflection intensity lower than the intensity threshold value are filtered.
The crop point cloud is in a scattered and disordered state, and in order to quickly remove noise, a spatial topological relation among point cloud data needs to be established. In the embodiment, a KD-Tree method is adopted to divide point cloud, then a near point of a certain point is searched, and a point which is too far away from the near point is marked as a noise point to be removed. The algorithm is described as follows:
(a) Reading in three-dimensional point cloud data;
(b) Expanding the point cloud data into KD-Tree;
(c) After the point cloud space division is finished, searching any point p in the point cloud space division i K nearest neighbor of
Wherein->
Is p i A set of K nearest neighbors;
(d) Calculating p i In close proximity of K
Mean value of the distances of the points in>
(e) Will be provided with
And a set threshold value D σ Compared, if->
Then consider point p i And deleting the noise points, otherwise, keeping the noise points.
(f) And repeating the steps until all the points in the point cloud are processed.
Through the steps, the noise points in the point cloud space can be removed quickly.
(II) establishing a digital surface model and a digital ground model:
in the embodiment, a small amount of point clouds are extracted from the preprocessed point cloud data to serve as key control points, and a digital surface model s and a digital ground model are established by utilizing the key control points and through an irregular triangulation network.
(1) Establishing a digital surface model s:
dividing an XOY horizontal plane of a coordinate system O-XYZ where preprocessed wheat field point cloud data are located into a plurality of cells with the same size, establishing a plurality of cuboids by taking each cell as a bottom surface, respectively searching point clouds with the maximum elevation values in each cuboid, taking the searched point clouds as first key control points, mapping all the first key control points to the XOY plane, and establishing a digital surface model s based on the first key control points on the XOY plane by adopting a divide-and-conquer Delaunay irregular triangulation network generation method; the method comprises the following specific steps:
(a) Putting all the first key control points mapped on the XOY plane into a point set, and dividing the point set into two sub-point sets with approximately equal points;
(b) Recursively constructing a Delaunay triangulation network in the sub-point set;
(c) Respectively connecting the bottom line and the top line of a Delaunay triangulation network formed by the two sub-point sets;
(d) Merging the Delaunay triangulation networks of the two sub-point sets from the bottom line to the top line, and optimizing the merged triangulation networks by utilizing a Lawson local optimization algorithm to obtain the triangulation network according with the Delaunay rule;
(e) And (d) mapping the triangulation network finally generated in the step (d) into an O-XYZ three-dimensional space coordinate system (i.e. remapping the first key control point on the triangulation network into the O-XYZ three-dimensional space coordinate system) to obtain a digital surface model s.
Through the steps, the irregular triangular net s representing the plant height can be generated as shown in figure 3.
(2) Establishing a digital ground model:
since the ground point cloud is sparse, the scanned wheat field point cloud data is divided into cells of smaller size in the horizontal plane, namely: dividing an XOY horizontal plane of a coordinate system O-XYZ where preprocessed wheat field point cloud data are located into a plurality of cells with consistent sizes and smaller sizes, establishing a plurality of cuboids by taking each cell as a bottom surface, respectively searching point clouds with minimum height values in each cuboid, taking the searched point clouds as second key control points, mapping all the second key control points to the XOY plane, establishing a digital ground model based on the second key control points on the XOY plane by adopting a divide-and-conquer Delaunay irregular triangulation network generation method, wherein the establishing process of the digital ground model is similar to the establishing process of the digital surface model, and the specific steps are as follows:
(a) Putting all second key control points mapped on the XOY plane into a point set, and dividing the point set into two sub-point sets with approximately equal points;
(b) Recursively constructing a Delaunay triangulation network in the set of sub-points;
(c) Respectively connecting a bottom line and a top line of a Delaunay triangulation network formed by the two sub-point sets;
(d) Merging the Delaunay triangulation networks of the two sub-point sets from the bottom line to the top line, and optimizing the merged triangulation network by using a Lawson local optimization algorithm to obtain a triangulation network according with the Delaunay rule;
(e) And (d) mapping the triangulation network finally generated in the step (d) to an O-XYZ three-dimensional space coordinate system to obtain the digital ground model.
(III) establishing a spike layer space model:
after obtaining the digital surface model s, according toPerforming statistical analysis on the field measured value to determine the height h required to be expanded above and below the digital surface model s u And h d Then s + h u And s-h d And (3) forming the upper and lower boundaries of the initial spike layer space model together, and removing all point clouds which do not belong to the initial spike layer space model from the preprocessed point cloud data without clustering.
Calculating an optimal height threshold h on the preliminarily established spike layer space model by adopting an Otsu threshold method t From h by t Separating the point clouds of the stem leaves and the ears of the crops, removing the point clouds of the stem leaves of the crops, and further reducing the point clouds of the non-ears.
And drawing a histogram according to the height values of all the point clouds, wherein the abscissa is the height value, and the ordinate is the number of the point clouds with corresponding heights. Most of the laser mainly irradiates the crop ears, part irradiates the stem leaves, and the other part irradiates the ground, so the point cloud is distributed on the height histogram as shown in figure 4.
Optimal height threshold h t The calculation method comprises the steps of setting a smaller initial value for a height threshold h, defining two classes on two sides of the threshold on a point cloud distribution histogram, wherein the left side is a stem leaf point cloud with less point clouds, the right side is a crop ear point cloud (namely the wheat ear point cloud) with more point clouds, and the proportions of the total number of the stem leaf point clouds and the crop ear point cloud in the histogram are respectively marked as omega 1 ，ω 2 The mean values are respectively recorded as μ 1 And mu 2 By calculating the between-class variance σ b 2 To find the optimal height threshold h t Between-class variance σ b 2 As shown in the formula:
the value of h is from small to large, and all values are traversed by a certain step length to ensure that sigma is increased b 2 Taking h of the maximum value as the optimal height threshold value h t . Below the height threshold h in the initial spike layer spatial model t The point clouds are regarded as the point clouds reflected by stems, leaves and the ground, all the point clouds are removed, and the rest point is the actual spike layer space modelPoint cloud data within.
(IV) confirming the ear of wheat (i.e. the ear of the crop):
(1) Point cloud clustering:
the wheat ears are represented as point cloud clusters with high density in the point clouds, the point clouds in the actual ear layer space model are divided into the point cloud clusters by adopting a density-based clustering algorithm DBSCAN method, and the number of the point cloud clusters and each point cloud data in the point cloud clusters are obtained. In order to improve the operation efficiency and avoid calculating the distance between a certain point and all the points when judging whether the certain point is a core point, the point cloud in the spike layer space model is divided by using KD-Tree, then clustering is performed by using DBSCAN, and when k adjacent points of the certain point need to be calculated, searching is performed by using KD-Tree, so that the operation efficiency of the system is improved. And outputting the data of the point cloud cluster by the DBSCAN algorithm, and primarily obtaining the number of the crop ears. In the DBSCAN calculation process, a neighborhood density threshold MinPts is directly set, and then a radius parameter epsilon is calculated in a self-adaptive mode. The point cloud clustering process is shown in fig. 5, and the specific algorithm is as follows:
(a) Giving a specific k value, dividing point clouds in an actual spike layer space model by using the KD-Tree, establishing a KD-Tree index, calculating the distance between all the point clouds and the kth adjacent point cloud thereof, and calculating the k value according to d k (p) d represents arranging all point clouds in descending order k (p), search for d k (p) points of inflection p in a curve formed by descending order t D is corresponded to k (p t ) As the radius parameter ∈. d is a radical of k (p) inflection points p in a curve formed by descending order t The point on the left side is a larger distance from its own k-th neighbor, let p t All point clouds on the left side are recorded as noise, and point clouds on the right side are recorded as core points. The radius parameter e of all points is set as d k (p t ) The neighborhood density threshold MinPts is set to k.
(b) P is to be t All core points on the right side are marked as unaccessed points.
(c) And randomly extracting one non-access point p, marking the non-access point p as an access point, creating a point cloud cluster C, and adding the point cloud cluster C into the point cloud cluster C.
(d) A set N is established, incorporating points within p neighborhoods (i.e., regions that are less than e from point p) into the set N.
(e) And randomly drawing a point p 'in N, and if the p' is an unaccessed point, marking the point as an access point. Comparing the distance between p ' and the k-th adjacent point with the radius parameter belonging to the family, if the distance value is larger than the radius parameter belonging to the family, the p ' is a boundary point, otherwise, the p ' is a core point; and if p 'is a core point, adding p' neighborhood points into the set N, and if p 'is not classified into the point cloud cluster C, classifying the p' neighborhood points into the point cloud cluster C.
(f) And (e) if no point has been accessed in N, continuing the step (e) until no point has been accessed in N, and then executing the step (g).
(g) And storing the information of all the points in the point cloud cluster C and outputting the point cloud cluster C.
(h) If p is t Returning to the step (c) when there is no access point on the right side, and executing the steps (c) to (g) until p t And outputting all point cloud clusters C without accessing points on the right side.
Finally, outputting the point cloud clusters and the number of the point cloud clusters and the coordinates of each point cloud in the point cloud clusters by the DBSCAN, temporarily marking all the point cloud clusters as wheat ear point cloud clusters, namely, one point cloud cluster represents one wheat ear, namely, the number of the preliminary wheat ears is obtained through clustering.
(2) Correcting the number of the wheat head point cloud clusters:
firstly, correcting the number of point cloud clusters, taking the corrected number of point cloud clusters as the actual number of wheat ears, and taking the average value of point cloud coordinates in each point cloud cluster as each wheat ear coordinate, wherein the number of wheat ears is the number of wheat plants; the correction steps are as follows:
(a) Projecting each point cloud data in the point cloud clusters on a horizontal plane, acquiring convex closure bags projected by each point cloud cluster in the horizontal plane, if the convex closure bags of a certain two point cloud clusters are overlapped, setting the convex closure bags of the two point cloud clusters as A and B respectively, and calculating the overlapped areas of the convex closure bags and the B;
(b) Setting a coincidence area threshold T overlap If the overlapping area of the convex closure of some two point cloud clusters is equal to max { S } A ,S B The ratio of (A) to (B) is greater than or equal to T overlap Respectively calculating the average value of Z-axis coordinates of all point clouds in the two point cloud clusters, and regarding the point cloud cluster with the lower average value of the Z-axis coordinates as the stem leaf of the plant of wheatTo be eliminated, wherein S A ，S B Areas of A, B, respectively; if the overlapping area of the convex closure of some two point cloud clusters is equal to max { S } A ,S B The ratio of is less than T overlap If the lower point cloud cluster is the wheat ear of a short plant, continuously keeping the two point cloud clusters;
(c) And (3) not taking the point cloud clusters removed in the step as the wheat ears (namely the crop ears), taking the rest of the point cloud clusters as the wheat ears, taking the number of the point cloud clusters left after removal as the actual number of the wheat ears, and taking the average value of X-axis coordinates, the average value of Y-axis coordinates and the average value of Z-axis coordinates of each point cloud in the point cloud clusters as the coordinate values of X-axis, Y-axis and Z of the corresponding wheat ears respectively.
The wheat comprises stem leaves and wheat ears, the wheat ears are positioned at the topmost end of the wheat, the stem leaves are positioned below the wheat ears, the wheat ear point cloud clusters marked by mistake can be removed through the steps, the number correction of the wheat ear point cloud clusters, namely the correction of the number of the wheat ears, is realized, the number of the wheat ears is obtained, and the number of the wheat ears is the number of the wheat plants.
(3) Calculating the actual height of wheat corresponding to the wheat ear:
and substituting the X-axis coordinate value and the Y-axis coordinate value of the wheat head into the digital ground model to obtain a corresponding ground height value, and subtracting the corresponding ground height value from the Z-axis coordinate value of the wheat head to obtain the actual height of the wheat corresponding to the wheat head.
(V) establishing a three-dimensional model of field wheat:
and drawing a three-dimensional view according to the obtained number of the wheat plants, the actual height of the wheat and the coordinates of the wheat ears at the top of the wheat to obtain a three-dimensional model of the wheat field, namely a three-dimensional model of crops in the field, and displaying the position, the height and the total number of each wheat plant.
In the embodiment, the ears of the crops are directly identified through the three-dimensional point cloud data, the height, the number and the coordinate information of the ears of the crops can be obtained, the labor cost is greatly reduced compared with a manual measurement method, the accuracy of crop parameter estimation can be improved compared with a method for measuring the information of single crop and counting the crop information of a field, and an enough data basis is provided for field management.
The scope of the present invention includes, but is not limited to, the above embodiments, and the present invention is defined by the appended claims, and any alterations, modifications, and improvements that may occur to those skilled in the art are all within the scope of the present invention.
Claims (5)
1. A field crop three-dimensional reconstruction method based on laser radar point cloud is characterized by comprising the following steps:
(1) Scanning field crops through a laser radar to obtain point cloud data, and performing data preprocessing on the point cloud data output by the laser radar;
(2) Extracting a small amount of point clouds serving as key control points from the preprocessed point cloud data, and establishing a digital surface model s and a digital ground model by utilizing the key control points and an irregular triangular network;
(3) Presetting a height value h to be expanded above the digital surface model s u And the height value h to be expanded below d S + h u And s-h d The upper boundary and the lower boundary of the initial spike layer space model are formed together, and point cloud data which do not belong to the initial spike layer space model are filtered from the preprocessed point cloud data; calculating the optimal height threshold h by adopting an Otsu threshold method t Removing the part in the space model of the initial spike layer which is lower than the height threshold h t Obtaining point cloud data in the actual spike layer space model;
(4) Dividing the point cloud data in the actual spike layer space model into a plurality of point cloud clusters by adopting a density-based clustering algorithm DBSCAN method, and obtaining the number of the point cloud clusters and each point cloud data in the point cloud clusters;
(5) Correcting the number of the point cloud clusters, taking the corrected number of the point cloud clusters as the number of the crop ears, and taking the average value of each point cloud coordinate in each point cloud cluster as the coordinate of each crop ear, wherein the number of the crop ears is the number of the crops;
(6) Bringing the X-axis coordinate value and the Y-axis coordinate value of the crop ear into the digital ground model to obtain a corresponding ground height value, and subtracting the corresponding ground height value from the Z-axis coordinate value of the crop ear to obtain the actual height of the crop corresponding to the crop ear;
(7) Drawing a three-dimensional view according to the obtained number of the crops, the actual height of each crop and the coordinates of the crop ears at the top of each crop to obtain a three-dimensional model of the crops in the field;
the step (2) comprises the following steps:
dividing an XOY horizontal plane of a coordinate system O-XYZ where preprocessed point cloud data are located into a plurality of cells with consistent sizes, establishing a plurality of cuboids by taking each cell as a bottom surface, respectively searching point clouds with the largest elevation value in each cuboid, taking the found point clouds as first key control points, mapping all the first key control points onto the XOY plane, and establishing a digital surface model s by adopting a Delaunay irregular triangulation network generation method on the basis of the first key control points on the XOY plane;
dividing an XOY horizontal plane of a coordinate system O-XYZ where the preprocessed point cloud data are located into a plurality of cells with the same size, establishing a plurality of cuboids by taking each cell as a bottom surface, respectively searching a point cloud with the minimum height value in each cuboid, taking the found point cloud as a second key control point, mapping all the second key control points onto the XOY plane, and establishing a digital ground model based on the second key control points on the XOY plane and by adopting a Delaunay irregular triangulation network generation method.
2. The method for three-dimensional reconstruction of field crops based on lidar point cloud as claimed in claim 1, wherein the step (1) comprises:
unmanned aerial vehicle carries on laser radar, establishes laser radar coordinate system O respectively L -X L Y L Z L And unmanned aerial vehicle coordinate system O V -X V Y V Z V ；
The unmanned aerial vehicle flies through the field and the laser radar on the unmanned aerial vehicle scans field crops to obtain point cloud data of the field crops, and the point cloud data output by the laser radar is converted from a polar coordinate to a laser radar Cartesian coordinate system O L -X L Y L Z L In the laser radar coordinate system O L -X L Y L Z L Inner pointCloud data is converted into unmanned aerial vehicle coordinate system O V -X V Y V Z V In, the unmanned plane coordinate system O V -X V Y V Z V Converting the point cloud data in the space into a world coordinate system O-XYZ;
setting an intensity threshold, filtering point cloud data with the reflection intensity smaller than the intensity threshold in the point cloud data, and removing noise points in the point cloud data by adopting a KD-Tree point cloud denoising method.
3. The method for three-dimensional reconstruction of field crops based on lidar point cloud according to claim 1, wherein the establishing of the digital surface model s based on the first key control point on the XOY plane and by using the Delaunay irregular triangulation network generation method specifically comprises:
(a) Putting all the first key control points mapped on the XOY plane into a point set, and dividing the point set into two sub-point sets with equal points;
(b) Constructing a Delaunay triangulation network in the sub-point set;
(c) Respectively connecting the bottom line and the top line of a Delaunay triangulation network formed by the two sub-point sets;
(d) Merging the Delaunay triangulation networks of the two sub-point sets from the bottom line to the top line, and optimizing the merged triangulation networks by utilizing a Lawson local optimization algorithm to obtain the triangulation network according with the Delaunay rule;
(e) And (d) mapping the triangulation network finally generated in the step (d) to an O-XYZ three-dimensional space coordinate system to obtain a digital surface model s.
4. The method for three-dimensional reconstruction of field crops based on lidar point cloud of claim 1, wherein the establishing of the digital ground model based on the second key control point on the XOY plane and by using the Delaunay irregular triangulation network generation method specifically comprises:
(a) Putting all the second key control points mapped on the XOY plane into a point set, and dividing the point set into two sub-point sets with equal points;
(b) Constructing a Delaunay triangulation network in the sub-point set;
(c) Respectively connecting the bottom line and the top line of a Delaunay triangulation network formed by the two sub-point sets;
(d) Merging the Delaunay triangulation networks of the two sub-point sets from the bottom line to the top line, and optimizing the merged triangulation networks by utilizing a Lawson local optimization algorithm to obtain the triangulation network according with the Delaunay rule;
(e) And (d) mapping the triangulation network finally generated in the step (d) into an O-XYZ three-dimensional space coordinate system to obtain the digital ground model.
5. The method for three-dimensional reconstruction of field crops based on lidar point cloud as claimed in claim 1, wherein the step (5) comprises:
(a) Projecting each point cloud data in the point cloud clusters on a horizontal plane, acquiring convex closure bags projected by each point cloud cluster in the horizontal plane, if the convex closure bags of a certain two point cloud clusters are overlapped, setting the convex closure bags of the two point cloud clusters as A and B respectively, and calculating the overlapped areas of the convex closure bags and the B;
(b) Setting a coincidence area threshold T overlap If the overlapping area of the convex closure bags of some two point cloud clusters is equal to max { S } A ,S B The ratio of } is greater than or equal to T overlap Respectively calculating the average value of Z-axis coordinates of each point cloud in the two point cloud clusters, and removing the point cloud clusters with low average value of the Z-axis coordinates, wherein S A ，S B A, B, respectively; if the overlapping area of the convex closure of some two point cloud clusters is equal to max { S } A ,S B The ratio of is less than T overlap Then continuing to reserve the two point cloud clusters;
(c) The removed point cloud clusters are not used as crop ears, the rest point cloud clusters are used as crop ears, and the corrected number of the point cloud clusters is the number of the crop ears;
(d) And respectively taking the average value of the X-axis coordinates, the average value of the Y-axis coordinates and the average value of the Z-axis coordinates of each point cloud in the point cloud cluster as the coordinate values of the X-axis, the Y-axis and the Z of the corresponding crop ear.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910632756.0A CN110363861B (en) | 2019-07-14 | 2019-07-14 | Laser radar point cloud-based field crop three-dimensional reconstruction method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910632756.0A CN110363861B (en) | 2019-07-14 | 2019-07-14 | Laser radar point cloud-based field crop three-dimensional reconstruction method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110363861A CN110363861A (en) | 2019-10-22 |
| CN110363861B true CN110363861B (en) | 2023-04-07 |
Family
ID=68219401
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910632756.0A Active CN110363861B (en) | 2019-07-14 | 2019-07-14 | Laser radar point cloud-based field crop three-dimensional reconstruction method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110363861B (en) |
Families Citing this family (14)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113496461A (en) * | 2020-03-18 | 2021-10-12 | 广州极飞科技股份有限公司 | Point cloud data processing method and device, computer equipment and storage medium |
| US11668799B2 (en) * | 2020-03-20 | 2023-06-06 | Aptiv Technologies Limited | Histogram based L-shape detection of target objects |
| CN111487646A (en) * | 2020-03-31 | 2020-08-04 | 安徽农业大学 | Online detection method for corn plant morphology |
| CN111723866A (en) * | 2020-06-19 | 2020-09-29 | 新石器慧通（北京）科技有限公司 | Point cloud clustering method and device, unmanned vehicle and readable storage medium |
| CN111915721B (en) * | 2020-07-28 | 2024-01-26 | 广州中科智巡科技有限公司 | Method and system for extracting horizontal section of power transmission line corridor building based on laser point cloud |
| WO2022095035A1 (en) * | 2020-11-09 | 2022-05-12 | 深圳市大疆创新科技有限公司 | Data processing method and movable platform |
| CN112669462B (en) * | 2021-01-18 | 2023-02-28 | 新拓三维技术(深圳)有限公司 | Model processing method and system suitable for scanning point cloud |
| CN113343917A (en) * | 2021-06-30 | 2021-09-03 | 上海申瑞继保电气有限公司 | Histogram-based substation equipment identification method |
| CN113780144A (en) * | 2021-09-06 | 2021-12-10 | 广西大学 | Crop plant number and stem width automatic extraction method based on 3D point cloud |
| CN116540255A (en) * | 2022-01-26 | 2023-08-04 | 上海飞机制造有限公司 | System and method for measuring and obtaining plane shape by using multiple laser radars |
| CN114782650A (en) * | 2022-05-09 | 2022-07-22 | 京能（锡林郭勒）发电有限公司 | Three-dimensional reconstruction method and system for closed coal yard |
| CN115063677B (en) * | 2022-06-10 | 2023-10-10 | 安徽农业大学 | Wheat Tian Daofu degree identification method and device based on point cloud information |
| CN115100272B (en) * | 2022-06-17 | 2023-09-22 | 浙江大学 | Prefabricated part point cloud data set manufacturing method for deep learning segmentation network |
| CN115901621A (en) * | 2022-10-26 | 2023-04-04 | 中铁二十局集团第六工程有限公司 | Digital identification method and system for concrete defects on outer surface of high-rise building |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105574929A (en) * | 2015-12-15 | 2016-05-11 | 电子科技大学 | Single vegetation three-dimensional modeling method based on ground LiDAR point cloud data |
| CN106845360A (en) * | 2016-12-27 | 2017-06-13 | 郑州大学 | High-resolution crop surface model construction method based on unmanned aerial vehicle remote sensing |
| CN108198230A (en) * | 2018-02-05 | 2018-06-22 | 西北农林科技大学 | A kind of crop and fruit three-dimensional point cloud extraction system based on image at random |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107918753B (en) * | 2016-10-10 | 2019-02-22 | 腾讯科技（深圳）有限公司 | Processing Method of Point-clouds and device |
-
2019
- 2019-07-14 CN CN201910632756.0A patent/CN110363861B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105574929A (en) * | 2015-12-15 | 2016-05-11 | 电子科技大学 | Single vegetation three-dimensional modeling method based on ground LiDAR point cloud data |
| CN106845360A (en) * | 2016-12-27 | 2017-06-13 | 郑州大学 | High-resolution crop surface model construction method based on unmanned aerial vehicle remote sensing |
| CN108198230A (en) * | 2018-02-05 | 2018-06-22 | 西北农林科技大学 | A kind of crop and fruit three-dimensional point cloud extraction system based on image at random |
Non-Patent Citations (2)
| Title |
|---|
| Development of High-Throughput Field Phenotyping System Using Imagery from Unmanned Aerial Vehicle;Ryo Sugiura 等;《ASABE》;20151231;第1页-第8页 * |
| 基于车载三维激光雷达的玉米叶面积指数测量;张漫 等;《农业机械学报》;20190630;第50卷(第6期);第12页-第21页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110363861A (en) | 2019-10-22 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110363861B (en) | Laser radar point cloud-based field crop three-dimensional reconstruction method | |
| CN106815847B (en) | Trees dividing method and single tree extracting method based on laser radar point cloud | |
| CN112489212B (en) | Intelligent building three-dimensional mapping method based on multi-source remote sensing data | |
| CN106529469B (en) | Unmanned aerial vehicle-mounted LiDAR point cloud filtering method based on self-adaptive gradient | |
| CN107369161A (en) | A kind of workpiece point cloud segmentation method at random based on the European cluster of improvement | |
| CN105513127A (en) | Rod-shaped object regular three-dimensional modeling method and rod-shaped object regular three-dimensional modeling system based on density peak clustering | |
| CN106599915B (en) | A kind of vehicle-mounted laser point cloud classifications method | |
| CN111340723B (en) | Terrain-adaptive airborne LiDAR point cloud regularization thin plate spline interpolation filtering method | |
| CN112070769A (en) | Layered point cloud segmentation method based on DBSCAN | |
| CN106651863A (en) | Point cloud data based automatic tree cutting method | |
| CN116310192A (en) | Urban building three-dimensional model monomer reconstruction method based on point cloud | |
| CN114332366A (en) | Digital city single house point cloud facade 3D feature extraction method | |
| CN109166145A (en) | A kind of fruit tree leaf growth parameter(s) extracting method and system based on cluster segmentation | |
| CN115205690B (en) | Method and device for extracting street tree in monomer mode based on MLS point cloud data | |
| CN111598780B (en) | Terrain adaptive interpolation filtering method suitable for airborne LiDAR point cloud | |
| CN110208815A (en) | A kind of large area maturation crop harvest information fast acquiring method based on airborne laser radar | |
| CN113538264B (en) | Denoising method and device for point cloud data and storage medium | |
| CN111950589B (en) | Point cloud region growing optimization segmentation method combined with K-means clustering | |
| CN109766824A (en) | Main passive remote sensing data fusion classification method based on Fuzzy Evidence Theory | |
| CN115187803A (en) | Positioning method for picking process of tender shoots of famous tea | |
| CN116523898A (en) | Tobacco phenotype character extraction method based on three-dimensional point cloud | |
| CN115410036A (en) | Automatic classification method for key element laser point clouds of high-voltage overhead transmission line | |
| CN116051841A (en) | Roadside ground object multistage clustering segmentation algorithm based on vehicle-mounted LiDAR point cloud | |
| CN115661398A (en) | Building extraction method, device and equipment for live-action three-dimensional model | |
| CN105321168B (en) | A kind of method of the automatic compilation mountain area raised path through fields in three-dimensional laser point cloud |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |