A kind of automatic measurement method of coal storage quantity in coal bunker based on laser radar

Abstract

This invention discloses an automated method for measuring the amount of coal stored in a coal bunker based on lidar, belonging to the field of lidar point cloud technology. The specific technical solution is as follows: A full-coverage scan of the coal bunker is performed using lidar to acquire the original point cloud data; then, the lidar point cloud data undergoes preprocessing such as thinning and noise reduction; next, the rotation matrix and translation vector between each pair of lidar points are calculated using a spherical target to synthesize the point cloud; local voids in the point cloud are detected and repaired by interpolating 3D data points using neighborhood information; then, the point cloud is segmented according to regions, and the projection surface is fitted using the least squares method to obtain two periods of data: the segmented region point cloud and the corresponding projection surface; finally, the volume of coal piles in each region of the coal bunker is calculated using surface integral, and the weight of coal is calculated based on the density of different coal types. The weight of the entire coal bunker is then summed, representing the total amount of coal stored in the bunker during that period, meeting the needs for safe, efficient, and accurate coal inventory.

Description

Technical Field

[0001] This invention belongs to the field of radar point cloud technology, specifically relating to an automated method for measuring the amount of coal stored in a coal bunker based on lidar. Background Technology

[0002] Coal, as a vital energy resource, is costly and abundant, thus requiring rapid and accurate inventory checks. Coal yard inventory checks involve statistically analyzing the amount of coal in stock and comparing it with the recorded inventory data. This is essential for verifying the amount of coal to be fed into the furnace and understanding the total coal stock in the yard; it is a crucial component of coal yard supervision and management. Due to its large workload and low accuracy, coal yard inventory checks have always been a challenge for coal-fired power plants, severely impacting production planning. Therefore, quickly and accurately obtaining coal yard inventory levels has become a pressing issue.

[0003] Existing methods for calculating coal inventory mainly include two types. One is manual measurement, which involves first shaping the piled coal using a bucket wheel stacker-reclaimer (the carrier used to stack and reclaim coal), typically into a relatively regular trapezoidal or rectangular shape, and then measuring it manually with a tape measure. This method has a large measurement error and is labor-intensive, seriously affecting work efficiency. The other method is drone measurement. Although drone measurement solves the problems of high labor intensity and poor working environment, it is still difficult to accurately count the coal inventory inside the coal yard due to its difficult operation and the high dust concentration inside the coal yard. Summary of the Invention

[0004] To meet the needs of coal storage bases for real-time inventory of coal storage in coal bunkers, this invention provides an automated measurement method for coal storage in coal bunkers based on lidar. By using lidar to perform a full-coverage scan of the coal yard, the change volume of the coal yard is calculated based on fine three-dimensional point cloud data, thereby realizing the calculation of the storage amount of all coal in the coal bunker.

[0005] To achieve the above objectives, the technical solution adopted by this invention is: an automated method for measuring the amount of coal stored in a coal bunker based on lidar, the specific steps of which are as follows:

[0006] S1. Acquire lidar point cloud data, arrange the laser scanner reasonably in the internal space of the coal yard, perform full-coverage scanning of the coal yard, collect point cloud data at regular intervals, and generate three-dimensional coordinate data.

[0007] S2. Point cloud preprocessing: preprocessing the original lidar point cloud data by downsampling, noise reduction, etc.

[0008] Furthermore, step S2 includes the following steps:

[0009] S21. Point cloud downsampling: The original LiDAR point cloud data is thinned using a voxelization downsampling method. This method creates a uniform downsampled point cloud from the input point cloud by using a regular voxel grid. Specifically: ① The point cloud is loaded into a voxel grid; ② The midpoints of each occupied voxel are averaged to obtain an accurate point.

[0010] S22. Point cloud denoising: Statistical filtering is used to denoise the downsampled point cloud data. This method calculates the average distance from each point to all its neighboring points. It is assumed that the result is a Gaussian distribution, the shape of which is determined by the mean and standard deviation. Points whose average distance is outside the standard range (defined by the global distance mean and variance) are defined as outliers and removed from the dataset.

[0011] S3. Point cloud registration: Using a spherical target, the rigid transformation is obtained through singular value decomposition, and the rotation matrix is ​​calculated. Translation matrix This allows the source point cloud to rotate and translate to overlap with the target point cloud, thus achieving point cloud registration.

[0012] Furthermore, step S3 includes the following steps:

[0013] S31. Establish a spatial rectangular coordinate system for coordinate transformation between radars. Utilize the principle of constructing a rectangular coordinate system using three mutually perpendicular edges sharing a common vertex, and establish a spatial rectangular coordinate system for transformation using five sets of points in space.

[0014] S32. Automatic identification of spherical targets: First, the spherical points of the spherical target are extracted from the global point cloud collected by each lidar, i.e., feature extraction of the spherical target; then, the points of the extracted spherical points are identified one-to-one according to the radar relationship, i.e., automatic correspondence of target points.

[0015] For spherical target feature extraction, an automatic spherical point extraction algorithm is used. All points with average curvature and Gaussian curvature greater than 0 are selected, and a random point is chosen. The neighborhood is calculated, and a local spherical surface is fitted. The distance from each point in the neighborhood to the center of the sphere is considered. If the given conditions are met, the point is considered a spherical point; calculate all points whose coordinate distances to the center of the sphere meet the radius condition, and output the spherical points.

[0016] The automatic mapping of points on the spherical target involves designing algorithms based on the target's structure to identify the extracted spherical points according to radar relationships. The identification results are then sequentially stored in a database for use in subsequent rotation matrix R and translation vector mappings between the source and target point clouds. Parameter calculation.

[0017] S33. The spherical fitting algorithm calculates the center of the sphere. Considering the influence of noise, statistical filtering preprocessing is performed on the point cloud data before fitting the sphere, the process being the same as step S22 above. Then, the least squares method is used to fit the three-dimensional spherical data, aiming to minimize the error of the fitted radius in the mean square sense. The objective function is as follows:

[0018]

[0019] In the formula, It is the error of the fitted radius. Given spherical coordinates, and These are the fitted coordinates of the sphere's center and radius.

[0020] The sum of squares of the errors is:

[0021]

[0022] In the formula, yes The function, therefore, let To each The center of the sphere can be determined by the partial derivative being equal to 0. and radius .

[0023] S34. Solve for the rotation matrix R and translation vector. Find the rotation and translation matrix between two point clouds to transform the source point cloud into the same coordinate system as the target point cloud.

[0024] Among them, solving for the rotation matrix R and the translation vector Using singular value decomposition and least-squares rigid registration, the optimization solution of the rotation matrix R and translation vector t is equivalent to minimizing The objective function is as follows:

[0025]

[0026] In the formula, R is the rotation matrix. It is a translation vector. express A set of rotation groups in 3D space, It is a transformation space. and These are two sets of points in the transformation space. ( >0) represents the weight of each pair of points, which is determined using eigenvalue decomposition. Calculate the rotation matrix, which is as follows:

[0027]

[0028] in, and They are respectively and an orthogonal unit matrix It is called a left singular matrix. It is called a right singular matrix. It is a diagonal matrix, and the elements on the diagonal are singular values.

[0029] Translation vector ,in, , It is the weighted centroid of two point sets.

[0030] S35. Using point cloud synthesis sequences, multiple radar point clouds are stitched together. First, a breadth-first search approach is employed to synthesize the point clouds from the multiple radars to be synthesized. The output is a sequence. In point cloud stitching, the composite sequence and the rotation matrix R and translation vector between radars calculated in step S34 above are used. This transforms the source point cloud to the same coordinate system as the target point cloud.

[0031] Regarding the breadth-first search algorithm: [This refers to a specific radar system / radar]. As the initial point of the sequence, the radars are accessed sequentially. Each adjacent vertex is visited until all adjacent vertices of the visited vertex have been visited. If there is still radar at this point... If it has not been accessed, select this unaccessed radar. As a new starting point, repeat the above process so that all radars All were accessed.

[0032] S4. Point Cloud Hole Repair: Due to the complex shape of the object being measured, limitations of the scanning equipment, or external occlusion, holes may appear. Failure to address these holes could lead to algorithm failure. This invention first uses Poisson reconstruction to triangulate the point cloud, then estimates the normal vectors of the point cloud based on a local surface fitting method, and finally repairs the holes through point cloud upsampling.

[0033] Furthermore, step S4 includes the following steps:

[0034] S41. Point cloud 3D reconstruction: First, the data is downsampled and preprocessed, and the processing is the same as step S21 above. The synthesized data is processed into a uniform point cloud, and on this basis, the Poisson reconstruction algorithm is used to reconstruct the point cloud surface.

[0035] The Poisson surface reconstruction method can transform discrete point cloud data into a continuous surface model. The main idea of ​​this method is to solve for the shape of the surface using the Poisson equation, thereby realizing the reconstruction of point cloud data.

[0036] S42. Point cloud normal vector estimation, as an important geometric surface feature of 3D point clouds, is key input information for subsequent point cloud upsampling. This embodiment uses the local surface fitting method to estimate the point cloud normal vector, specifically: for each point in the point cloud... With radius Or you can directly calculate the nearest neighbor. To obtain the neighborhood of a point, solving for the normal vector is equivalent to minimizing... ,in, For the centroid of the neighborhood, ( ) represents the decentralized point. It is the matrix transpose symbol. For normal vector, ( )and The dot product represents the product of the points on the map. Projection in the direction. Let ,but:

[0037]

[0038] in, Points within the neighborhood Given the covariance matrix, solve for the eigenvectors, select the eigenvector corresponding to the smallest eigenvalue, and normalize it. This vector is then the point cloud normal vector.

[0039] S43. Point Cloud Upsampling: This algorithm progressively upsamples the point cloud, simultaneously moving towards edge points to generate a denser point set from the input set. It primarily uses three parameters: the edge sensitivity parameter controls the insertion position of new points. A higher edge sensitivity value prioritizes insertion along sharp features, with a range of [0,1]; the sharpness angle parameter controls the smoothness of sharp features, with a higher value resulting in smoother edges, and a range of [0,90]; and the neighborhood radius parameter represents the number of neighboring points, calculated using a fixed sphere. Used to limit the number of return points in the neighborhood of each ball.

[0040] S5. Convert the radar coordinate system to ground coordinates, that is, convert the composite point cloud coordinates after hole repair to ground coordinates based on the horizontal ground of the coal yard. Coordinate plane, ground upwards A spatial rectangular coordinate system along the axis facilitates subsequent point cloud clipping.

[0041] Specifically: Place three spherical targets at the three vertices of an acute triangle, and calculate the centers of the three targets. , , The process of calculating the center of the sphere is the same as described in steps S32 and S33 above; then, the distance between the three centers of the sphere is calculated using the distance calculation formula. , , ), at any point ( (Taking point 1 as an example) is the origin of the coordinate system, and its coordinates are represented as (0,0,0). Point for Points on the axis, with coordinates represented as (0, ... The coordinates of point C are represented as (0, 0); , );

[0042] in, for by The height corresponding to the base is calculated using Heron's formula: , in the formula for The half perimeter, , , , These are the side lengths of the triangle. Then... , , The coordinates of the sphere's center are transformed with the calculated corresponding ground coordinates, i.e., the rotation matrix is ​​calculated. Translation vector The calculation process is the same as step S34 above.

[0043] S6. Point Cloud Cropping: Considering the different types of coal in the coal bunker, the coordinates of the boundary points of the coal pile scanned by radar are taken and recorded according to the region. Then, the point cloud of the coal pile in the coal bunker is cropped according to the specified region for subsequent coal weight calculation. If there are obstacles such as pillars in the middle of the cropped coal pile, a region growing segmentation algorithm is used to segment the pillar point cloud to avoid the influence of obstacles on subsequent algorithms.

[0044] S7. Point Cloud Projection: The clipped point cloud is projected onto the horizontal ground of the coal bunker. The coal pile volume is calculated using the point cloud data from both the clipping and projection phases. The least squares method is used to fit the projection plane, i.e., using the data from the horizontal ground... points ( The fitted plane (≥3 and not collinear) has the following equation: .

[0045] Assumption ≠0, then, .make , , ,

[0046] According to the least squares principle, to minimize the sum of squared errors, the objective function is as follows:

[0047]

[0048] In the formula, , , Given the coordinates of a point on the ground. , , It is an unknown, , , Treating it as the independent variable and differentiating it, we can obtain the coefficients of the plane equation. , , , The value of .

[0049] S8. Storage Calculation: When inventorying a large, irregularly shaped coal pile, the principles of calculus can be used to divide the large coal pile into several smaller coal piles at predetermined steps. The height of each smaller coal pile is then measured to calculate its volume. The volumes of all smaller coal piles are then summed to determine the total volume of the entire coal pile. The mass of the coal pile in each area is then calculated using a mass formula, and finally, the summation yields the amount of coal stored in the coal bunker.

[0050] Furthermore, step S8 includes the following steps:

[0051] S81. Unit block discretization: The discretization step size is determined based on the point cloud density, and the irregular coal pile 3D model is divided into... Each unit block;

[0052] S82. Calculation of coal storage in the coal bunker; calculation of the top area of ​​each unit block. and base area The average of the top area and the bottom area is taken as the grid area. The calculation formula is as follows:

[0053]

[0054] in, The serial number of each unit block. .

[0055] The average height of the elevation values ​​of the point clouds on the top and bottom surfaces of each unit block in the 3D model is taken as the maximum elevation value of the unit block. and minimum value The relative height H of each unit block is obtained by subtracting the top elevation value from the bottom elevation value. The calculation formula is as follows:

[0056]

[0057] The volume of each discretized block is obtained by multiplying the grid area by the relative height of the block. The volume of the coal pile is then calculated and summed to obtain the change volume. This refers to the volume of coal pile in a certain area of ​​the coal bunker during that period, calculated using the following formula:

[0058]

[0059] Based on the coal pile cut out in step S6 above, and considering the coal pile density in different areas, the weight of the coal pile in the coal bunker is calculated using the mass formula as follows:

[0060]

[0061] in, These are the serial numbers of the various areas within the coal bunker. .

[0062] Compared with the prior art, the beneficial effects of the present invention are:

[0063] I. This invention introduces lidar technology into the dynamic monitoring of coal storage in coal yards, enabling refined information management of coal yards, real-time monitoring of coal yard dynamic operations, providing advanced intelligent means to ensure the accuracy of coal inventory results, and improving fuel operation and maintenance efficiency and intelligent management level.

[0064] Second, this invention uses a method to obtain the center of a spherical target for point cloud stitching. By using the target center coordinates to calculate rotation and translation parameters, the point cloud data is converted to the target coordinate system. This method can leverage the advantages of easy deployment and high accuracy, and can effectively improve the efficiency and quality of point cloud stitching.

[0065] Third, this invention employs a series of methods, such as point cloud thinning, three-dimensional reconstruction, point cloud normal vector estimation, and point cloud upsampling, to fill in the holes in the point cloud after the registration of the entire coal yard. This avoids the influence of external occlusion and other factors on the point cloud data, making the calculation results of the coal storage in the coal yard more reliable.

[0066] Fourth, this invention discretizes the three-dimensional model of irregular coal piles into unit blocks. By determining the optimal step size, the volume change of each unit block can be accurately extracted, thereby realizing the calculation of coal storage in the coal yard and providing a basis for real-time monitoring of coal storage in the coal yard. Attached Figure Description

[0067] Figure 1 This is a flowchart illustrating the automated measurement method for coal storage in a coal bunker based on radar equipment, as per the present invention.

[0068] Figure 2 This is a real-life image of the coal bunker in a specific implementation method.

[0069] Figure 3The image shows the registered and cropped point cloud data of the coal pile in a specific implementation.

[0070] Figure 4 This is a point cloud image of the coal pile after the hole surface has been repaired in a specific implementation method. Detailed Implementation

[0071] To make the technical problems, solutions, and beneficial effects of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0072] like Figure 1 As shown, an automated method for measuring the amount of coal stored in a coal bunker based on lidar includes the following steps:

[0073] S1. LiDAR point cloud data acquisition: First, the dimensions of the coal storage area are measured. Then, based on performance indicators such as the effective scanning range of the LiDAR, the LiDAR is rationally deployed, and data is collected periodically. In this embodiment, point cloud data is acquired by periodically collecting point cloud data using radar deployed in a coal bunker of a coal mine. The radar is installed at a height of approximately 11m, and its coverage range extends up to 65m within the bunker. The specific LiDAR parameters used in the implementation are shown in Table 1 below:

[0074] Table 1 Technical parameters of lidar equipment

[0075]

[0076] S2. Point cloud preprocessing: preprocessing the original lidar point cloud data by downsampling, noise reduction, etc.

[0077] Furthermore, step S2 includes the following steps:

[0078] S21. Point Cloud Downsampling: The acquired raw LiDAR point cloud data is thinned using a voxel-based downsampling method. This method creates a uniform downsampled point cloud from the input point cloud using a regular voxel grid. Specifically: ① The point cloud is loaded into a voxel grid; ② The midpoints of each occupied voxel are averaged to obtain a precise point. In this embodiment, the voxel size is specified as 1 cm for voxel-based downsampling.

[0079] S22. Point Cloud Denoising: Statistical filtering is used to denoise the downsampled point cloud data. This method calculates the average distance from each point to all its neighbors. It is assumed that the result follows a Gaussian distribution, the shape of which is determined by the mean and standard deviation. Points whose average distance falls outside the standard range (defined by the global average and variance) are defined as outliers and removed from the dataset. In this embodiment, the number of neighboring points is specified as 30, and the average distance is calculated based on this. The standard deviation multiplier for the average distance of the point cloud is selected as 0.05.

[0080] S3. Point cloud registration: Using a spherical target, the rigid transformation is obtained through singular value decomposition, and the rotation matrix is ​​calculated. Translation matrix This allows the source point cloud to rotate and translate to overlap with the target point cloud, thus achieving point cloud registration.

[0081] Furthermore, step S3 includes the following steps:

[0082] S31. Establish a spatial rectangular coordinate system for inter-radar coordinate transformation. Utilizing the principle of constructing a rectangular coordinate system with three mutually perpendicular edges sharing a common vertex, a spatial rectangular coordinate system for transformation is established using four sets of points in space (each three points defining a mutually perpendicular plane). This embodiment uses five sets of target points to construct the coordinate system; adding points on the diagonal is to improve the robustness of the target and increase the synthesis accuracy.

[0083] S32. Automatic identification of spherical targets. First, the spherical points of the spherical targets are extracted from the global point cloud collected by each lidar, i.e., feature extraction of spherical targets; then, the positions of the extracted spherical points are identified one-to-one according to the radar relationship, i.e., automatic target position correspondence.

[0084] For feature extraction of the spherical target, an automatic point extraction algorithm is used. The average curvature and Gaussian curvature of the points are calculated, and all points with average curvature K > 0 and Gaussian curvature H > 0 are selected. A point is randomly selected, its neighborhood is calculated, and a local spherical surface is fitted. The distance from each point in the neighborhood to the center of the sphere is calculated. If the conditions are met: ( For radius, If the radius error threshold is used, then the point is considered a spherical point; calculate all points whose coordinate distances to the center of the sphere meet the radius condition, and output the spherical points.

[0085] The automatic mapping of spherical target points involves identifying the extracted spherical points one-to-one according to radar relationships, which is then used for the subsequent rotation matrix between the source and target point clouds. Translation vector Parameter calculation. The relevant algorithm steps designed in this embodiment are:

[0086] 1) The dbscan clustering algorithm is used to cluster the data into labeled classes. In this embodiment, the maximum distance between points in the same cluster is 100mm, and the minimum number of points in an effective cluster is 50.

[0087] 2) Classify the clustered data. Use distance to determine the high point, the origin (the point directly below the high point), and the points on the diagonal; then use the right-hand rule of the cross product to determine the remaining two points. The cross product formula is:

[0088]

[0089] when When, vector to vector It rotates counterclockwise; when When, vector to vector The rotation is clockwise, ultimately yielding one-to-one target data.

[0090] S33. The spherical fitting algorithm calculates the center of the sphere. Considering the influence of noise, statistical filtering preprocessing is performed on the point cloud data before fitting the sphere, the process being the same as step S22 above. Then, the least squares method is used to fit the three-dimensional spherical data, aiming to minimize the error of the fitted radius in the mean square sense. The objective function is as follows:

[0091]

[0092] In the formula, It is the error of the fitted radius. Given spherical coordinates, and These are the fitted center coordinates and radius of the sphere, respectively.

[0093] The sum of squares of the errors is:

[0094]

[0095] In the formula, yes The function, therefore let To each The center of the sphere can be determined by the partial derivative being equal to 0. and radius .

[0096] S34. Solving for the rotation matrix Translation vector Find the rotation and translation matrix between two point clouds to transform the source point cloud into the same coordinate system as the target point cloud.

[0097] The embodiment includes scanning data from two adjacent lidars A and B, and common spherical targets C1, C2, C3, C4, and C5. The 3D coordinates of C1, C2, C3, C4, and C5 in lidar A are extracted using point fitting. Similarly, the 3D coordinates of the five spherical targets in lidar B are extracted. A point cloud registration formula is used. (in (That is, the corresponding points in the target point cloud and the source point cloud), transforming the point cloud coordinates of radar B to the coordinate system of radar A.

[0098] Among them, solving the rotation matrix Translation vector We employ least-squares rigid registration using singular value decomposition to optimize the rotation matrix R and translation vector. Equivalent to minimizing The objective function is as follows:

[0099]

[0100] in, It's a rotation matrix. It is a translation vector. express A set of rotation groups in 3D space, It is a transformation space. and These are two sets of points in the transformation space. ( >0) represents the weight of each pair of points, which is determined using eigenvalue decomposition. Calculate the rotation matrix, which is as follows:

[0101]

[0102] in, and They are respectively and an orthogonal unit matrix It is called a left singular matrix. It is called a right singular matrix. It is a diagonal matrix, and the elements on the diagonal are singular values;

[0103] Translation vector ,in, , It is the weighted centroid of two point sets.

[0104] S35. Using point cloud synthesis sequences, multiple radar point clouds are stitched together. First, a breadth-first search approach is adopted to synthesize the point clouds from the multiple radars to be synthesized. The output is a sequence. In point cloud stitching, the rotation matrix between radars is calculated according to the synthesized sequence and step S34 above. Translation vector This transforms the source point cloud to the same coordinate system as the target point cloud.

[0105] Regarding the breadth-first search algorithm: [This refers to a specific radar system / radar]. As the initial point of the sequence, the radars are accessed sequentially. Each adjacent vertex of the visited vertex has been visited until all adjacent vertices of the visited vertex have been visited. If there are still radars at this point... If it has not been accessed, select this unaccessed radar. As a new starting point, repeat the above process so that all radars All were accessed.

[0106] S4. Point cloud hole repair: Due to the complex shape of coal piles in the coal bunker, as well as the limitations of the scanning equipment itself or external obstruction, holes may appear in the point cloud. If hole repair is not performed, the algorithm may fail. Figure 4 The diagram shows the effect of point cloud hole repair in this embodiment. The specific implementation process is as follows: first, the point cloud is triangularly meshed using the Poisson reconstruction algorithm; then, the normal vector of the point cloud is estimated based on the local surface fitting method; and finally, the point cloud hole repair is achieved through upsampling.

[0107] Furthermore, step S4 includes the following steps:

[0108] S41. Point cloud 3D reconstruction: First, the data is downsampled and preprocessed, and the processing is the same as step S21 above. Voxel downsampling processes the synthesized data into a point cloud with relatively uniform density, and on this basis, the Poisson reconstruction algorithm is used to reconstruct the point cloud surface.

[0109] The Poisson surface reconstruction method can transform discrete point cloud data into a continuous surface model. The main idea is to solve for the shape of the surface using the Poisson equation, thereby reconstructing the point cloud data. Specific implementation steps include defining an octree, setting the function space, creating a vector field, solving the Poisson equation, and extracting isosurfaces. In this embodiment, the point cloud data sampling ratio is 0.2, the number of iterations in the algorithm is set to 9, the octree depth is 8, and the minimum number of points required in each octree node is set to 40.

[0110] S42. Point cloud normal vector estimation, as an important geometric surface feature of 3D point clouds, is key input information for subsequent point cloud upsampling. This embodiment uses the local surface fitting method to estimate the point cloud normal vector, specifically: for each point in the point cloud... With radius Or you can directly calculate the nearest neighbor. Each point obtains its neighborhood (nearest neighbors in this embodiment). (Setting the value to 50), solving for the normal vector is equivalent to minimizing ,in For the centroid of the neighborhood, ( ) represents the decentralized point. It is the matrix transpose symbol. For normal vector, ( )and The dot product represents the product of the points on the map. Projection in the direction. Let ,but:

[0111]

[0112] in, Points within the neighborhood Given the covariance matrix, solve for the eigenvectors, select the eigenvector corresponding to the smallest eigenvalue, and normalize it. This vector is then the point cloud normal vector.

[0113] S43. Point cloud upsampling: The point cloud is progressively upsampled, while simultaneously approaching edge singularities (detected through normal changes), generating a denser point set from the input point set. In this embodiment, the specific settings for the three parameters of this algorithm are as follows: the edge sensitivity parameter is set to 0, indicating uniform insertion of the point cloud; the sharpness angle parameter controls the smoothness of sharp features, with a value of 25, which can better preserve the shape of the coal pile; the neighborhood radius parameter is the number of neighboring points, calculated using the radius of a fixed sphere. This is used to limit the number of neighborhood return points for each sphere (to avoid too many points in high-density areas). If the parameter setting is too small, it may not be possible to insert new points. In this embodiment, the neighborhood point search radius is set to 15mm, which can better fill the holes in the point cloud.

[0114] S5. Convert the radar coordinate system to ground coordinates, that is, convert the composite point cloud coordinates after hole repair to ground coordinates based on the horizontal ground of the coal yard. Coordinate plane, ground upwards A spatial rectangular coordinate system along the axis facilitates subsequent point cloud clipping.

[0115] Specifically: Place three spherical targets at the three vertices of an acute triangle, and calculate the centers of the three targets. , , The process of calculating the center of the sphere is the same as described in steps S32 and S33 above; then, the distance between the three centers of the sphere is calculated using the distance calculation formula. , , ). At any point ( (Taking point 1 as an example) is the origin of the coordinate system, and its coordinates are represented as (0,0,0). Point for Points on the axis, with coordinates represented as (0, ... The coordinates of point C are represented as (0, 0); , ).

[0116] in, for by The height corresponding to the base is calculated using Heron's formula: , in the formula for The half perimeter, , , , These are the side lengths of the triangle. Then... , , The coordinates of the sphere's center are transformed with the calculated corresponding ground coordinates, i.e., the rotation matrix is ​​calculated. Translation vector The calculation process is the same as step S34 above.

[0117] S6. Point Cloud Cropping: Considering the different types of coal in the coal bunker, the coordinates of the boundary points of the coal pile scanned by radar are taken and recorded according to the region. Then, the point cloud of the coal pile in the coal bunker is cropped according to the specified region for subsequent coal weight calculation. If there are obstacles such as pillars in the middle of the cropped coal pile, a region growing segmentation algorithm is used to segment the pillar point cloud to avoid the influence of obstacles on subsequent algorithms.

[0118] Furthermore, step S6 includes the following steps:

[0119] S61. Crop the point cloud of a specified region. In this embodiment, the point cloud is... The coordinate input parameters determine a cube bounding box (i.e., the coordinates of the boundary points of the coal pile). The axis defines the height of the selected area (in the coal pile point cloud). (Determine the maximum and minimum values ​​of the coordinates), and cut the coal pile area in the coal bunker according to the bounding box.

[0120] S62. Obstacle segmentation: Here, the obstacles are the pillars in the cut coal pile area. In this embodiment, a region growth segmentation algorithm is used to segment the pillar point cloud. Specifically: (1) Calculate the normal and curvature of the point cloud, and sort the curvature values ​​from largest to smallest; (2) Select the point with the smallest curvature and put it into the seed point sequence, and search its neighborhood; (3) Determine the growth conditions, calculate the angle between the normal of each adjacent point and the normal of the current seed point. If the angle is less than the threshold, add the current point to the current region, and then judge whether the point less than the normal threshold is less than the curvature threshold. If it is less, add it to the seed sequence; (4) Growth terminates. Repeat the above steps until there are no points that meet the conditions, and then the growth ends. In this embodiment, the number of neighborhood search points is set to 30, the smoothing threshold is 30°, and the curvature threshold is 0.05.

[0121] S7. Point cloud projection: The clipped point cloud is projected onto the horizontal ground of the coal bunker. The coal pile volume is calculated using the point cloud data from both the clipping and projection phases. This embodiment uses the least squares method when calculating the projection plane.

[0122] Using the horizontal ground points ( The plane equation obtained by fitting (≥3 and not collinear) is: Assuming ≠0, then, .make , , According to the least squares principle, to minimize the sum of squared errors, the objective function is as follows:

[0123]

[0124] In the formula, , , Given the coordinates of a point on the ground. , , It is an unknown, , , Treating it as the independent variable and differentiating it, we obtain the system of equations:

[0125]

[0126] The coefficients of the plane equation can be obtained by calculating the determinant using Cramer's rule. , , , The value of .

[0127] S8. Storage Calculation: When inventorying a large, irregularly shaped coal pile, the surface integral principle can be used to divide the large coal pile into several smaller coal piles at a defined step size. The height of each smaller coal pile is then measured to calculate its volume. Finally, the volumes of all smaller coal piles are summed to complete the measurement of the entire coal pile volume. In this embodiment, the coal pile in the coal bunker is divided into regions. The volume of each region is calculated using the surface integral principle, and the mass of the coal pile in each region is calculated using a mass formula. Finally, the volumes are summed to obtain the amount of coal stored in the coal bunker.

[0128] Furthermore, step S8 includes the following steps:

[0129] S81. Unit block discretization: The discretization step size is determined based on the point cloud density, and the irregular coal pile 3D model is divided into... Each unit block.

[0130] S82. Calculation of coal storage in the coal bunker; calculation of the top area of ​​each unit block. and base area The average of the top area and the bottom area is taken as the grid area. The calculation formula is as follows:

[0131]

[0132] in, The serial number of each unit block. .

[0133] The average height of the elevation values ​​of the top and bottom point clouds corresponding to each unit block in the 3D model is taken as the maximum elevation value of the unit block. and minimum value The relative height H of each unit block is obtained by subtracting the top elevation value from the bottom elevation value. The calculation formula is as follows:

[0134]

[0135] The volume of each discretized block is obtained by multiplying the grid area by the relative height of the block. The volume of the coal pile is then calculated and summed to obtain the change volume. This refers to the volume of coal pile in a certain area of ​​the coal bunker during that period, calculated using the following formula:

[0136]

[0137] Based on the coal pile cut out in step S6 above, and considering the coal pile density in different areas, the weight of the coal pile in the coal bunker is calculated using the mass formula as follows:

[0138]

[0139] in, These are the serial numbers of the various areas within the coal bunker. .

[0140] Using this measurement method to Figure 2 The coal bunker and coal pile shown are subjected to automated inventory checks. Figure 3 The image shows point cloud data of the cut coal pile surface, as well as the projected point cloud file of the coal pile projected onto the ground. Figure 4 This embodiment illustrates the preprocessed and partially repaired 3D point cloud data of the coal pile (including both point cloud and projection data). Based on this, using the surface integral principle and considering the point cloud density of the coal pile, the 3D model of the coal pile is discretized into blocks according to a certain step size, and then the volume is calculated. The corresponding actual total weight of the coal pile is 2100. (The density of medium-strength coal is calculated as 1.3) The calculation results of the coal bunker's changing volume and coal storage amount obtained by selecting different time lengths are shown in the table below:

[0141] Table 2 Calculation Results of Coal Inventory

[0142]

[0143] The calculation results in the table above show that when the step size is small, the unit matching ratio is low, making it impossible to accurately calculate the mining volume. When the step size is large, although the matching degree of each unit block between the two models is high, the relative height error of the selected blocks is large, and the estimated mining volume is larger than the actual mining volume. When the step size of the processed coal bunker point cloud data is set to 0.10m to 0.13m, the error is within 5%. From the perspective of calculation error and running time, this method is a reliable automated coal inventory method.

[0144] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included within the scope of the present invention.

Claims

1. An automated method for measuring the amount of coal stored in a coal bunker based on lidar, characterized in that, Specific steps: S1. Acquire lidar point cloud data, perform a full-coverage scan of the coal yard using a laser scanner, collect point cloud data periodically, and generate three-dimensional coordinate data; S2. Point cloud preprocessing: The original LiDAR point cloud data undergoes downsampling and noise reduction preprocessing, specifically including the following steps: S21. Point cloud downsampling: The original LiDAR point cloud data is thinned using a voxelization downsampling method. This method creates a uniform downsampled point cloud from the input point cloud by using a regular voxel grid. Specifically: ① The point cloud is loaded into a voxel grid; ② The midpoints of each occupied voxel are averaged to obtain an accurate point. S22. Point cloud denoising: Statistical filtering is used to denoise the downsampled point cloud data. This method calculates the average distance from each point to all its neighboring points. It is assumed that the result is a Gaussian distribution, the shape of which is determined by the mean and standard deviation. Points whose average distance is outside the standard range are defined as outliers and removed from the dataset. S3. Point cloud registration: Using a spherical target, the rigid transformation is obtained through singular value decomposition, and the rotation matrix is ​​calculated. Translation matrix This causes the source point cloud to rotate and translate until it coincides with the target point cloud, thus achieving point cloud registration. S4. Point cloud hole repair: Poisson reconstruction is used to triangular mesh the point cloud, and the normal vector of the point cloud is estimated based on the local surface fitting method. Hole repair is achieved by upsampling the point cloud. Specifically, the following steps are included: S41. Point cloud 3D reconstruction: First, the data is downsampled and preprocessed, and the processing is the same as step S21 above. The synthesized data is processed into a uniform point cloud, and on this basis, the Poisson reconstruction algorithm is used to reconstruct the point cloud surface. The Poisson surface reconstruction method transforms discrete point cloud data into a continuous surface model, and solves for the shape of the surface using the Poisson equation, thereby reconstructing the point cloud data. S42. Point cloud normal vector estimation: The local surface fitting method is used to estimate the point cloud normal vector. Specifically, for each point in the point cloud... With radius Or you can directly calculate the nearest neighbor. To obtain the neighborhood of a point, solving for the normal vector is equivalent to minimizing... ,in, For the centroid of the neighborhood, ( ) represents the decentralized point. It is the matrix transpose symbol. For normal vector, ( )and The dot product represents the product of the points on the map. Projection in the direction, let ,but: ； in, Points within the neighborhood Given the covariance matrix, solve for the eigenvectors, select the eigenvector corresponding to the smallest eigenvalue, and normalize it. This vector is then the point cloud normal vector. S43. Point Cloud Upsampling: This algorithm progressively upsamples the point cloud, simultaneously moving towards edge points, to generate a denser point set from the input set. It primarily uses three parameters: an edge sensitivity parameter, which controls the insertion position of new points (higher sensitivity values ​​indicate higher priority for insertion along sharp features; the value range is [0,1]); an sharpness angle parameter, which controls the smoothness of sharp features (higher values ​​result in smoother edges; the value range is [0,90]); and a neighborhood radius parameter, which represents the number of neighboring points and is calculated using a fixed sphere. Used to limit the number of return points in the neighborhood of each ball; S5. Transform the ground coordinates. Convert the radar coordinate system into ground coordinates. That is, transform the coordinates of the synthesized point cloud after hole repair, and finally transform it into a spatial rectangular coordinate system with the horizontal ground of the coal yard as the yz coordinate plane and the upward direction of the ground as the x-axis. S6. Point cloud cropping: Take the coordinates of the boundary points of the coal pile scanned by the radar, record them according to the region, and crop the point cloud of the coal pile in the coal bunker according to the specified region. S7. Point Cloud Projection: The clipped point cloud is projected onto the horizontal ground of the coal bunker. The coal pile volume is calculated using the point cloud data from both the clipping and projection phases. The projection plane is fitted using the least squares method, i.e., using the horizontal ground... On the surface A plane is fitted at points, and the equation of the plane is: ; Assumption ≠0, then, ,make , , According to the least squares principle, to minimize the sum of squared errors, the objective function is as follows: ； In the formula, , , Given the coordinates of a point on the ground. , , It is an unknown, , , Treating it as the independent variable and differentiating it, we can obtain the coefficients of the plane equation. , , , The value; S8. Storage Calculation: When inventorying a large, irregular coal pile, the principle of calculus is used to divide the large coal pile into several smaller coal piles at a determined step size. The height of each smaller coal pile is then measured to calculate the volume of each smaller coal pile. The volumes of the smaller coal piles are then summed to complete the measurement of the total coal pile volume. The mass of the coal pile in each area is then calculated using the mass formula. Finally, the volumes are summed to obtain the amount of coal stored in the coal bunker.

2. The automated measurement method for coal storage in a coal bunker based on lidar according to claim 1, characterized in that, Step S3 includes the following steps: S31. Establish a spatial rectangular coordinate system for coordinate transformation between radars. Utilize the principle of constructing a rectangular coordinate system using three mutually perpendicular edges sharing a common vertex, and establish a spatial rectangular coordinate system for transformation using five sets of points in space. S32. Automatic identification of spherical targets: First, the spherical points of the spherical targets are extracted from the global point cloud collected by each lidar, i.e., feature extraction of spherical targets. Then, the points of the extracted spherical points are identified one by one according to the radar relationship, i.e., automatic correspondence of target points. An automatic point extraction algorithm is used to extract features from the spherical target, filtering out all points with average curvature and Gaussian curvature greater than 0; a point is randomly selected, its neighborhood is calculated, and a local spherical surface is fitted; the distance from each point in the neighborhood to the center of the sphere is calculated. If the given conditions are met, the point is considered a spherical point; calculate all points whose coordinate distances to the center of the sphere meet the radius condition, and output the spherical points; Regarding the automatic mapping of spherical target points, based on the structure of the spherical target, the extracted spherical points are identified one-to-one according to radar relationships, and the identification results are sequentially stored in the database for subsequent rotation matrix mapping between the source point cloud and the target point cloud. Translation vector Parameter calculation; S33. The spherical fitting algorithm calculates the center of the sphere. Before fitting the sphere, the point cloud data is preprocessed by statistical filtering, the process being the same as step S22 above. Then, the least squares method is used to fit the three-dimensional spherical data, aiming to minimize the error of the fitted radius in the mean square sense. The objective function is as follows: ； In the formula, It is the error of the fitted radius. Given spherical coordinates, These are the fitted coordinates of the sphere's center. It is the radius of the fitted area; The sum of squares of the errors is: ； In the formula, yes The function, therefore, let To each The center of the sphere can be determined by the partial derivative being equal to 0. and radius ; S34. Solving for the rotation matrix Translation vector Find the rotation and translation matrix between two point clouds to transform the source point cloud into the same coordinate system as the target point cloud; Among them, solving the rotation matrix Translation vector We employ least-squares rigid registration using singular value decomposition to optimize the rotation matrix R and translation vector. Equivalent to minimizing The objective function is as follows: ； in, It's a rotation matrix. It is a translation vector. express A set of rotation groups in 3D space, It is a transformation space. and These are two sets of points in the transformation space. ( >0) represents the weight of each pair of points, which is determined using eigenvalue decomposition. Calculate the rotation matrix, which is as follows: ； in, and They are respectively and an orthogonal unit matrix It is called a left singular matrix. It is called a right singular matrix. It is a diagonal matrix, and the elements on the diagonal are singular values; Translation vector ,in, , It is the weighted centroid of two point sets; S35. Using point cloud synthesis sequences, multiple radar point clouds are stitched together. First, a breadth-first search approach is adopted to synthesize the point clouds from the multiple radars to be synthesized. The output is a sequence. In point cloud stitching, the rotation matrix between radars is calculated according to the synthesized sequence and step S34 above. Translation vector Transform the source point cloud into the same coordinate system as the target point cloud; Regarding the breadth-first search algorithm: [This refers to a specific radar...] As the initial point of the sequence, the radars are accessed sequentially. Each adjacent vertex of the visited vertex has been visited until all adjacent vertices of the visited vertex have been visited. If there are still radars at this point... If it has not been accessed, select this unaccessed radar. As a new starting point, repeat the above process so that all radars All were accessed.

3. The automated measurement method for coal storage in a coal bunker based on lidar according to claim 2, characterized in that, Specifically: Place three spherical targets at the three vertices of an acute triangle, and calculate the centers of the three targets. , , The process of calculating the center of the sphere is the same as described in steps S32 and S33 above; then, the distance between the three centers of the sphere is calculated using the distance calculation formula. , , ),by The point is the origin of the coordinate system, and its coordinates are represented as (0,0,0). Point for Points on the axis, with coordinates represented as (0, ... The coordinates of point C are represented as (0, 0); , ); in, for by The height corresponding to the base is calculated using Heron's formula: , in the formula for The half perimeter, , , , Let these be the side lengths of the triangle, and then... , , The coordinates of the sphere's center are transformed with the calculated corresponding ground coordinates, i.e., the rotation matrix is ​​calculated. Translation vector The calculation process is the same as step S34 above.

4. The automated measurement method for coal storage in a coal bunker based on lidar according to claim 1, characterized in that, Step S8 includes the following steps: S81. Unit block discretization: The discretization step size is determined based on the point cloud density, and the irregular coal pile 3D model is divided into... Each unit block; S82. Calculation of coal storage in the coal bunker; calculation of the top area of ​​each unit block. and base area The average of the top area and the bottom area is taken as the grid area. The calculation formula is as follows: ； in, The serial number of each unit block. ; The average height of the point cloud elevation values ​​on the top surface of each unit block in the 3D model is taken as the maximum elevation value of the unit block. The average height of the elevation values ​​of the point cloud on the bottom surface of each unit block in the 3D model is taken as the minimum elevation value of the unit block. The relative height H of each unit block is obtained by subtracting the top elevation value from the bottom elevation value. The calculation formula is as follows: ； The volume of each discretized block is obtained by multiplying the grid area by the relative height of the block. The volume of the coal pile is then calculated and summed to obtain the change volume. This refers to the volume of coal pile in a certain area of ​​the coal bunker during that period, calculated using the following formula: ； Based on the coal pile cut out in step S6 above, and considering the coal pile density in different areas, the weight of the coal pile in the coal bunker is calculated using the mass formula as follows: ； in, These are the serial numbers of the various areas within the coal bunker. .

Images