A three-dimensional surface reconstruction method for coal mine tunnel arch surface

Abstract

The invention discloses a three-dimensional surface reconstruction method of a coal mine tunnel arch surface, the steps are: flip the laser radar 90 degrees and install it on the guide rail of the mechanical arm base, drive the laser radar collector by the hydraulic motor moving device to collect the point cloud data of the coal mine tunnel arch surface according to the preset multiple sampling points; use multiple filters in the PCL library to remove outliers and simplify the point cloud data; use the PCA algorithm to estimate the point cloud normal; use the SAC-IA rough splicing and ICP fine splicing algorithm to splice and fuse multiple point cloud images; use the Greedy PT algorithm to realize the three-dimensional surface reconstruction of the coal mine tunnel arch surface, and perform hole identification and online repair on the reconstructed triangular mesh, so as to complete the reconstruction of the arch surface within the specified range, and the algorithm implementation of the whole process is automatically completed by the computer; wherein the Greedy PT algorithm has a simple principle, strong real-time performance, and fast calculation speed, and can meet the practical application of engineering scenes. The method has high point cloud splicing accuracy and good coal mine tunnel arch surface reconstruction effect.

Description

Technical Field

[0001] The invention relates to the technical field of three-dimensional curved surface reconstruction, in particular to a three-dimensional curved surface reconstruction method for a coal mine tunnel arch surface. Background Art

[0002] In recent years, the in-depth development of the integration of information technology and industrialization has made the three-dimensional reconstruction technology based on the engineering environment widely used in many engineering fields such as human-computer interaction, unmanned driving, and cultural relics protection. The three-dimensional reconstruction technology based on the PCL (Point Cloud Library) library involves the integration of multiple disciplines such as computer graphics, geometric calculation, and sensors. Since most of the traditional coal mine tunnel spraying work still adopts the method of manual spraying, there are problems such as low work efficiency and frequent safety accidents. To study an automatic anchor spraying robot; in order to meet the actual requirements of the spraying of the coal mine tunnel arch, it is necessary to first reconstruct the tunnel arch in three dimensions; among them, an explosion-proof three-dimensional laser radar is used to collect data on the tunnel arch, but due to the interference of the tunnel environment and the limitations of the laser radar itself, the point cloud data of the coal mine tunnel arch cannot be obtained all at once, and the tunnel arch will dilute some point cloud information, so it is impossible to efficiently complete the subsequent three-dimensional reconstruction of the coal mine tunnel arch. Therefore, the complex tunnel environment has put forward higher requirements on the installation position of the laser radar and the processing method of the point cloud data; the existing three-dimensional surface reconstruction algorithm is usually easily interfered by complex environmental factors, which makes the three-dimensional surface reconstruction effect of the coal mine tunnel arch poor and the algorithm robustness low.

[0003] In view of the shortcomings of the existing technology, there is a need for a method that can effectively collect coal mine tunnel arch surface data in real time and accurately, quickly and robustly complete the three-dimensional surface reconstruction of the coal mine tunnel arch surface. Summary of the invention

[0004] In order to solve the above problems, the present invention proposes a three-dimensional surface reconstruction method for the arch surface of a coal mine tunnel. The method can realize real-time acquisition of point cloud data of the tunnel arch surface at multiple sites and perform point cloud data preprocessing and point cloud stitching, and can also quickly and accurately complete the three-dimensional surface reconstruction of the arch surface of the coal mine tunnel.

[0005] In order to achieve the above technical objectives, the technical solution adopted by the present invention is as follows:

[0006] A method for reconstructing a three-dimensional surface of a coal mine tunnel arch surface uses a laser radar to scan the coal mine tunnel arch surface in real time to obtain point cloud data. The collected point cloud data is reconstructed into a three-dimensional surface according to the following steps:

[0007] Preprocessing: Flip the laser radar 90° and install it upside down on the guide rail of the robot arm base; collect tunnel arch data and set up multiple data collection stations; use the hydraulic motor moving device to drive the laser radar to move at a constant speed on the guide rail at a certain time interval to obtain complete tunnel arch point cloud data, and set the initial point and moving speed of the laser radar on the guide rail respectively; set the time interval for collecting point cloud data at each measuring station during the movement of the laser radar;

[0008] S1: Initialization, the laser radar on the guide rail of the robot base moves to the initial point, and the computer and the laser radar establish communication;

[0009] S2: The guide rail of the robot arm base drives the laser radar collector to move to the data sampling points on the guide rail at a certain speed. The laser collector collects real-time data on the arch surface of the coal mine tunnel at a preset time interval and transmits the collected point cloud data to the computer;

[0010] S3: The computer analyzes the point cloud data to determine whether there are noise, outliers and invalid points, and uses a filter to pre-process the point cloud data;

[0011] S4: Based on the results of the data preprocessing in the previous step, the computer uses the principal component analysis method to perform plane fitting on the N points in the K neighborhood of each point to estimate the normal vector and calculate the curvature at the same time;

[0012] S5: Multiple point cloud images collected from different sampling points on the guide rail are stitched according to the sequence of key point extraction, SAC-IA rough stitching based on FPFH feature descriptor, and ICP fine stitching, so that the intersecting areas between them overlap perfectly;

[0013] S6: The point cloud shape obtained from the arch surface of the coal mine tunnel is non-closed point cloud data. According to its characteristics, the greedy triangulation algorithm is used to realize the 3D surface reconstruction;

[0014] S7: After completing the 3D surface reconstruction of the tunnel arch, the polygonal holes generated by diluting part of the point cloud data on the tunnel arch are repaired online using a hole repair algorithm based on a triangular mesh;

[0015] S8: The guide rail at the base of the robot arm drives the laser radar to move, complete the three-dimensional surface reconstruction of the arch surface of the coal mine tunnel, and return to the initial point. After the grouting is completed, it enters the next arch surface reconstruction.

[0016] With the above design, the computer collects point cloud data of the coal mine tunnel arch surface through the laser radar collector, and performs point cloud data preprocessing based on the point cloud data obtained from each measuring station to remove noise, discrete points and invalid points caused by interference from factors such as the environment and the laser radar itself, thereby reducing the amount of point cloud data; after obtaining multiple streamlined point cloud images, the multiple point clouds are spliced ​​using point cloud data stitching technology to obtain complete point cloud data within a certain range of the coal mine tunnel arch surface; the computer performs three-dimensional reconstruction of the tunnel arch surface based on the complete point cloud data, and performs online repair of the polygonal voids generated after the reconstruction; the use of this method for three-dimensional reconstruction of complex coal mine tunnel arch surfaces can realize an automated process, and can repair triangular meshes that do not meet the requirements after reconstruction, prevent uneven thickness and missed spraying caused by direct grouting after reconstruction, reduce the re-spraying rate, and thus improve the speed and efficiency of the anchor spraying robot's spraying work.

[0017] As a preferred solution of the present invention, in steps S1 and S2, the initial moving point of the guide rail of the robot base is (x0, y0), the moving speed is S, the time interval for the laser radar collector to collect point cloud data is T, and the collector establishes communication between the computer and the laser radar through the IP address;

[0018] Among them, the guide rail of the base of the robotic arm drives the laser radar collector to collect point cloud data at a certain speed and time interval through a hydraulic motor moving device, and obtains the distance d between different data collection stations and the tunnel arch scanning width w. Multiple point cloud images will be generated within the entire guide rail range; the robotic arm is installed on a mobile platform trolley.

[0019] As a preferred solution of the present invention, in step S3, the specific steps of point cloud data preprocessing are:

[0020] S3-1: Use statistical filters to perform statistical analysis on the K neighborhood points of each point in the point cloud, and calculate the average distance from it to all neighboring points; then calculate the average value μ and standard deviation σ of the average distance of each point to determine the distance threshold thresh_d; according to the point cloud density distribution, remove the point cloud whose average neighborhood distance of a point is lower or higher than its threshold range;

[0021] Assume that the original point cloud dataset P = {p i (x i ,y i , z i )|1≤i≤n}, then the distance threshold is:

[0022] thresh_d=μ+m·σ

[0023] The point cloud dataset after outlier filtering is:

[0024] P′=(um*σ,u+m*σ)

[0025] Where m is the standard deviation multiple (usually 1-3);

[0026] S3-2: Import P′ into the pass filter, set multiple dimensional directions and different point cloud thresholds, filter out points outside the parameter range, and obtain point cloud data within the specified range;

[0027] S3-3: Import the point cloud data obtained in S3-2 into a voxel filter or a conditional filter, set a variety of voxel grid thresholds and other conditions of different specifications to simplify the point cloud data, obtain the point cloud preprocessing results, and greatly reduce the amount of computer operations.

[0028] By adopting the above scheme, the computer filters out the noise, outliers and invalid points in the point cloud data set, simplifies the original point cloud data, and calculates the point cloud data within a specific range, effectively reducing the number of point clouds and improving computing efficiency.

[0029] As a preferred solution of the present invention, in step S4, the specific steps of estimating the point cloud normal using the PCA algorithm are:

[0030] S4-1: For any point p in the point cloud i (1≤i≤n) query its k-domain point p ij (1≤j≤k), calculate p i The centroid of its k-neighborhood points

[0031]

[0032] S4-2: Construct the covariance matrix of local features of point cloud:

[0033]

[0034] S4-3: Calculate the eigenvalues ​​λ0, λ1, λ2 of the covariance matrix C, λ0≤λ1≤λ2, and the eigenvectors v0, v1, v2, where λ and v correspond one to one;

[0035] Among them, v0, v1, v2 are orthogonal, and the eigenvectors v1 and v2 determine the point p i An optimal tangent plane at, v0 is orthogonal to the tangent plane, so point p i Normal at n i It can be approximately represented by the eigenvector v0; λ0, λ1, λ2 are the degree of change in the direction of their respective eigenvectors;

[0036] S4-4: Calculate point p i The curvature τ pi :

[0037]

[0038] Using the above scheme, the computer estimates the normal and curvature within the k-neighborhood of the preprocessed point cloud and obtains the geometric features of the point cloud, so that feature point detection, point cloud stitching and surface reconstruction can be better performed later.

[0039] As an improved solution of the present invention, in step S5, the computer completes the point cloud image stitching according to the key point extraction, SAC-IA rough stitching based on FPFH feature descriptor and ICP fine stitching sequence for multiple point cloud images collected from different sampling points on the guide rail. The specific steps are:

[0040] S5-1: Extract key points of point cloud images using internal morphological descriptor algorithm:

[0041] Let point cloud P = {p i (x i ,y i , z i )|1≤i≤n}, for each point p i Establish a local coordinate system and set the search radius p r , determine p i is the center of the sphere, p r are all points within a sphere of radius ;

[0042] Calculate the weight ω ij :

[0043] |p i -p j |<p r

[0044] Calculate each point p i The covariance matrix of is:

[0045]

[0046] The eigenvalues ​​of the covariance matrix Arrange in descending order; set thresholds η1 and η2 to satisfy and The point is the key point, and iterate until all key points are found;

[0047] S5-2: Calculate the fast point feature histogram features of the point cloud to be spliced ​​and the target point cloud, and obtain each calculation point M p The relative relationship between all its neighboring points is used to establish a simple point feature histogram; the FPFH feature is calculated based on the SPFH feature, denoted as F(M P ):

[0048]

[0049] Among them, d i is the Euclidean distance of corresponding point pairs;

[0050] S5-3: Use the sampling consistency initial registration algorithm to stitch the model point cloud and the target point cloud. The algorithm steps are as follows: (1) Select m feature points to be registered in the model point cloud; (2) Find points in the target point cloud that are similar to the FPFH features of the model point cloud, and select the points with the closest distance as the corresponding relationship points; (3) Calculate the rotation and translation matrices of the corresponding point pairs; Use the Huber function to represent the distance error and function after the corresponding point pairs are rotated and translated, denoted as

[0051]

[0052] Where r is the distance threshold, ||m i || represents the Euclidean distance between the i-th group of corresponding points after transformation, minH(m i ) corresponds to the optimal transformation matrix after rough splicing;

[0053] S5-4: The roughly stitched source point cloud P' and the target point cloud G are used as input point clouds for ICP fine stitching, and the optimal node is used to optimize the KD tree to accelerate the search for corresponding point pairs; all points P' of the registration point cloud P' are treated as i , search for the nearest corresponding point G in the target point cloud G i , forming corresponding point pairs; calculating the rotation matrix R and the translation vector T so that the root mean square error M between the corresponding point pairs k Minimum:

[0054]

[0055] Finally, set the threshold α (i.e., M k -M k+1 <α) and the maximum number of iterations N max .

[0056] Using the above solution, the computer stitches multiple point cloud images collected from different sampling points on the guide rail of the robot arm base, and obtains complete point cloud data within 1.6m of the coal mine tunnel arch surface; thereby realizing fast and efficient point cloud stitching technology.

[0057] As an improved solution of the present invention, the specific steps of implementing three-dimensional surface reconstruction using the greedy projection triangulation algorithm in step S6 are:

[0058] As an improved solution of the present invention, a greedy projection triangulation algorithm is used in step S6 to realize the three-dimensional reconstruction of the arch surface of the coal mine tunnel. The specific steps are:

[0059] S6-1: Using the recursive algorithm of dynamic programming problem, the normal vector of the three-dimensional point obtained in step 4 is The coordinates of the target point O determine the local tangent plane equation of O to obtain the two-dimensional tangent plane formed by adjacent points; if O = (x0, y0, z0), Then the equation of the tangent plane through point O is as follows:

[0060] A(x-x0)+B(y-y0)+C(c-c0)=0

[0061] S6-2: Project the 3D point onto the 2D tangent plane. The projection point is stored in the projection matrix for rotation transformation calculation. The calculation process is:

[0062]

[0063] Among them, T Matrix is the translation transformation matrix:

[0064] Where x0, y0, z0 are the translations of each coordinate axis;

[0065] R x Expressed as a rotation matrix of α degrees around the x-axis:

[0066] R y It is expressed as a rotation matrix of θ degrees around the y-axis:

[0067] Combining the tangent plane equation in step S6-1 with the above equation, we can get any point P(x i ,y i , z i ) on the tangent plane O П Projection on:

[0068]

[0069] S6-3: Use the Delaunay-based spatial region growing algorithm to triangulate the point cloud obtained by projection in the plane. The triangulation process should meet the empty circumscribed circle characteristics and maximum and minimum angle criteria of the Delaunay algorithm. The algorithm selects a sample triangle as the initial surface, and continuously expands the surface boundary to obtain the connection relationship between each point and form a complete triangular mesh surface. Finally, the topological connection between the original three-dimensional points is determined according to the connection relationship of the projected point cloud. The obtained triangular mesh is the reconstructed surface model.

[0070] With the above scheme, the computer triangulates the point cloud data using a greedy projection triangulation algorithm and splices the triangular meshes built on the projection plane, thereby achieving three-dimensional surface reconstruction of the coal mine tunnel arch surface.

[0071] As a preferred solution of the present invention, in step S7, the polygonal holes generated by diluting part of the point cloud on the arch surface are repaired online using a hole repair algorithm based on a triangular mesh, and the specific steps are as follows:

[0072] S7-1: Perform boundary point detection on the tangent plane. If a point P in the point cloud data is a boundary feature point, then the k-neighborhood points of P should fall on its side; if P is an internal point, then the k-neighborhood points of P should fall around point P; the boundary feature points are judged by measuring the distribution uniformity of the k-neighborhood in the point cloud, and the maximum angle difference is used as a measure of distribution uniformity:

[0073] Definition i = (i = 1, 2, ..., k) is the projection point of the P neighborhood point on the tangent plane, and the nearest neighbor O of P is taken. i Form a line segment with P by As a benchmark, calculate Rotate clockwise to Angle Make an angle sequence Sort the angles to get a new angle sequence Define the angle sequence difference as

[0074] From L = (L1, L2, ..., L k ) to find the maximum angle sequence difference L max , as the basis for judging the boundary feature points; setting a threshold, when L max When it is greater than the threshold, P is a boundary feature point, otherwise P is an internal point;

[0075] S7-2: After finding the boundary feature points, the disordered feature points are ordered, the nearest neighbor points of the feature points are selected as connection points, and the iteration is repeated until a closed boundary line is formed;

[0076] S7-3: After obtaining the boundary line, the polygonal hole is repaired online according to the angle α between two adjacent edges in the hole boundary;

[0077] (1) When α≤90°, only one triangular patch is constructed;

[0078] (2) When 90°<α≤135°, construct two triangular patches to fill point Q and satisfy P i Q divides ∠P equally i-1 P iP i+1 ;

[0079] (3) When 135°<α≤200°, three triangular patches are constructed to fill points Q and R, and P should be satisfied. i Q and P i R divides ∠P into three equal parts i-1 P i P i+1 ;

[0080] (4) When α>200°, three triangular facets are constructed to fill points Q and R and satisfy ΔQP i+1 P i and ΔPP i R i-1 is an equilateral triangle;

[0081] S7-4: For the above hole repair methods, the newly added triangular facets need to be checked for legality; this can be done by calculating arrive The trend of δ1 and arrive It can be judged by the direction δ2. If δ2>δ1, the newly constructed triangle patch is unreasonable. It can be judged by checking whether the newly constructed triangle patch intersects with the original hole polygon. If they intersect, the generated triangle patch is unreasonable.

[0082] By adopting the above scheme, the computer identified the holes in the reconstructed triangular mesh model, repaired the existing polygonal holes online and checked the legitimacy of the newly added triangular facets; finally, the accuracy of the three-dimensional surface reconstruction of the coal mine tunnel arch surface was improved, and the reconstruction effect was better; finally, the steps of the above three-dimensional surface reconstruction algorithm were repeated to solve the engineering problem of three-dimensional surface reconstruction within multiple 1.6m ranges of the entire tunnel arch surface.

[0083] As an improved solution of the present invention, a three-dimensional reconstruction is performed on the area of ​​the coal mine tunnel arch surface. The installation position of the laser radar is located on the guide rail. The hydraulic motor moving device is used to drive the laser radar to move at a constant speed on the guide rail to collect point cloud data of the coal mine tunnel arch surface. After obtaining the point cloud data of the complete coal mine tunnel arch surface, the three-dimensional surface reconstruction is performed.

[0084] The technical effects of the present invention are as follows: a computer can collect point cloud data of a coal mine tunnel arch surface, perform data preprocessing, extract normal vectors and perform point cloud splicing in real time; a Greedy PT algorithm can be used to quickly realize three-dimensional surface reconstruction of a coal mine tunnel arch surface and perform online repair of polygonal voids generated therefrom; the Greedy PT algorithm makes the reconstruction process more stable and more accurate, and the Greedy PT algorithm has a simple principle, strong real-time performance and fast calculation speed, and can meet the practical application requirements of engineering scenarios. BRIEF DESCRIPTION OF THE DRAWINGS

[0085] Figure 1 It is a flow chart of a three-dimensional surface reconstruction method of a coal mine tunnel arch surface;

[0086] Figure 2 Flowchart for preprocessing point cloud data;

[0087] Figure 3 is the characteristic area map;

[0088] Figure 4 It is a non-feature area map;

[0089] Figure 5 A flowchart to complete the stitching of multiple point cloud data;

[0090] Figure 6 is the target neighborhood point map;

[0091] Figure 7 This is a diagram of the triangulation process of Bowyer-Watson point cloud data;

[0092] Figure 8 This is the principle diagram of identifying point P as a boundary feature point;

[0093] Fig. 9 This is the schematic diagram of point P being identified as an internal point;

[0094] Fig.10 It is a schematic diagram of the angle sequence L;

[0095] Fig.11 This is a schematic diagram of cavity repair;

[0096] Fig.12 To destroy the feature map of the hollow polygon;

[0097] Fig.13 It is the intersection graph of the newly constructed triangle patch and the hole polygon;

[0098] Fig.14 Schematic diagram of the structure of the reconstruction component for the simulated coal mine tunnel arch surface. DETAILED DESCRIPTION

[0099] The technical solution of the present invention is further described below in conjunction with the accompanying drawings and embodiments.

[0100] A three-dimensional surface reconstruction method for a coal mine tunnel arch surface, the process of the method is as follows Figure 1 As shown in the figure, the method uses a 16-line laser radar collector to scan the arch surface of the coal mine tunnel in real time to obtain point cloud data. The collected point cloud data is reconstructed into a three-dimensional surface according to the following steps:

[0101] Preprocessing: Flip the laser radar 90 degrees and install it upside down on the 1.6m long guide rail at the base of the robot arm; set up multiple data collection stations to obtain complete tunnel arch surface data within a range of 1.6m at one time; drive the laser radar to move at a constant speed on the guide rail at a certain time interval through a hydraulic motor moving device to obtain complete tunnel arch surface point cloud data, and set the initial point and moving speed of the laser radar on the guide rail respectively; set the time interval for collecting point cloud data at each measuring station during the movement of the laser radar;

[0102] S1: Initialization: the laser radar on the guide rail of the robot base moves to the initial point, and the computer and the laser radar establish communication.

[0103] S2: The guide rail at the base of the robotic arm drives the laser radar collector to move to the pre-designed data sampling points on the guide rail at a certain speed. The laser collector collects real-time data on the arch surface of the coal mine tunnel at preset time intervals and transmits the collected point cloud data to the computer.

[0104] S3: The computer analyzes the point cloud data to determine whether there are noise, outliers and invalid points, and uses filters such as StatisticalOutlierRemoval, PassThrough, and VoxelGrid to pre-process the point cloud data.

[0105] S4: Based on the results of the data preprocessing in the previous step, the computer uses principal component analysis (PCA) to perform plane fitting on the N points in the K neighborhood of each point to estimate the normal vector and calculate the curvature.

[0106] S5: Multiple point cloud images collected from different sampling points on the guide rail are stitched in the order of key point extraction, SAC-IA coarse stitching based on FPFH feature descriptor, and ICP fine stitching, so that the intersecting areas between them overlap perfectly.

[0107] S6: The point cloud shape obtained from the arch surface of the coal mine tunnel is a non-closed point cloud data. According to its characteristics, the greedy projection triangulation algorithm is used to realize the three-dimensional surface reconstruction.

[0108] S7: After the three-dimensional surface of the coal mine tunnel arch is reconstructed, the polygonal voids generated by the dilution of part of the point cloud data on the tunnel arch are repaired online using a void repair algorithm based on a triangular mesh.

[0109] S8: The guide rail of the robot arm base drives the laser radar to move to 1.6m, completing the 3D surface reconstruction of the coal mine tunnel arch within the range of 1.6m, and returns to the initial point. After the grouting within this range is completed, it enters the next arch reconstruction within the range of 1.6m.

[0110] In this embodiment, the above steps all use the Linux-Ubuntu 18.04 operating system, the development platform is ROS, the laser radar uses PandarXT-16, and the open source library is the PCL library.

[0111] Combination Figure 1 It can be seen that in step S1 and step S2, the initial movement point of the robot arm base rail is (x0, y0), the movement speed is S, the time interval for the laser radar collector to collect point cloud data is T, and the collector establishes communication between the computer and the laser radar through the IP address.

[0112] Among them, the 1.6m-long guide rail at the base of the robotic arm drives the lidar collector through a hydraulic motor moving device to collect point cloud data at a certain speed and time interval, and obtain the distance d between different data collection stations and the tunnel arch scanning width w. Multiple point cloud images will be generated within the entire 1.6m range.

[0113] Combination Figure 1 and Figure 2 It can be seen that in step S3, filters such as PassThrough, VoxelGrid, and StatisticalOutlierRemoval are used to preprocess the point cloud data. Point cloud preprocessing includes noise and outlier removal, data capture, and data simplification, specifically:

[0114] S3-1: First, load the point cloud data of the pcd file collected by the collector as input, and use the removeNaNFromPointCloud function to remove NAN (invalid) points; because too much data will reduce the computer's running speed, use the StatisticalOutlierRemoval statistical filter to remove outliers, set the neighborhood points to 50, and the distance threshold to 1.0;

[0115] S3-2: Secondly, the point cloud is imported into the PassThrough filter, which mainly operates in the z direction and sets the threshold to (-8, 8), thereby filtering out points outside the parameter range and obtaining point cloud data within the specified range;

[0116] S3-3: Data simplification mainly uses the VoxelGrid grid method, setting the voxel grid threshold to 1.0*1.0*1.0, and calculating the centroid value of all data points in each grid to replace all points in the grid, thereby greatly reducing the amount of computer calculations.

[0117] Combination Figure 1It can be seen that in step S4, the computer creates a normal estimation object through the NormalEstimation function, establishes a kd-tree data structure and sets the number of K nearest search points. The specific steps of the PCA algorithm to estimate the surface normal vector of the point cloud are:

[0118] S4-1: For any point p in the point cloud i (1≤i≤n) query its k-domain point p ij (1≤j≤k), calculate p i The centroid of its k-neighborhood points

[0119]

[0120] S4-2: Construct the covariance matrix of local features of point cloud:

[0121]

[0122] S4-3: The computer calculates the eigenvalues ​​λ0, λ1, λ2 (λ0≤λ1≤λ2) and eigenvectors v0, v1, v2 of the covariance matrix C, where λ and v correspond one to one;

[0123] Among them, v0, v1, v2 are orthogonal, and the eigenvectors v1 and v2 determine the point p i An optimal tangent plane at, v0 is orthogonal to the tangent plane, so point p i Normal at n i It can be approximately represented by the eigenvector v0; λ0, λ1, λ2 are the degrees of change in the direction of their respective eigenvectors.

[0124] S4-4: Computer calculation point p i The curvature τ pi :

[0125]

[0126] S4-5: Normal vector feature description:

[0127] like Figure 3 As shown in , the angle of the point cloud normal vector in the local area is larger, indicating that the geometric features of this area are more prominent; Figure 4 If the normal vector angle does not change much, it means that the area is relatively smooth. Therefore, we define a point p in the point cloud as i The arithmetic mean of the angles between the normal vectors of its neighboring points:

[0128]

[0129] In the formula, θ ij is the point cloud p i The angle with the neighboring normal vector.

[0130] According to point p i The angle between the normal vector of its neighboring point is used to extract the feature point, and an appropriate threshold δ is selected. i >δ, point p i The curvature is larger, so p i is a feature point; when f i <δ, point p i The curvature is small, so p i is a non-feature point.

[0131] In order to ensure the consistency of the direction of the normal vector of the point cloud surface and extract feature points more accurately, it is necessary to adjust its direction to satisfy the above formula:

[0132] n i *n j <0, (i≠j)

[0133] Where n i is the source point cloud normal vector, n j is the target point cloud normal vector.

[0134] Combination Figure 1 and Figure 5 It can be seen that in step S5, the computer first selects the tunnel arch surface point cloud data of the two adjacent starting stations for splicing, and then splices them with the point cloud of the next measuring station in turn after completion. The specific steps are:

[0135] S5-1: Extract key points of point cloud image using internal morphological descriptor algorithm ISSKeypoint3D function: Let point cloud P = {p i (x i ,y i , z i )|1≤i≤n}, for each point p i Establish a local coordinate system and set the search radius p r , determine p i is the center of the sphere, p r For all points in the sphere of radius Figure 6 shown.

[0136] Calculate the weight ω ij :

[0137] |p i -p j |<p r

[0138] Calculate each point p i The covariance matrix of is:

[0139]

[0140] The eigenvalues ​​of the covariance matrix Arrange in descending order; set thresholds η1 and η2 to satisfy and The point is the key point, and iterate until all key points are found;

[0141] S5-2: Create an FPFH feature descriptor, input point cloud data, establish a kd tree data structure and set K to search for the number of nearest neighbor points. After calculating these descriptors, save them in the FPFHSignature33 function for preliminary matching in the registration process.

[0142] S5-3: In order to solve the problem that ICP is prone to fall into local optimality, the SAC-IA algorithm is first used to roughly stitch the source point cloud P and the target point cloud G. The algorithm steps are: (1) Select m feature points to be registered from P; (2) Find points in G that are similar to the FPFH features of the source point cloud, and select the points with the closest distance as the corresponding relationship points; (3) Calculate the rotation and translation matrix of the corresponding point pairs; Use the Huber function to represent the distance error and function after the corresponding point pairs are rotated and translated, denoted as

[0143]

[0144] Where r is the distance threshold, ||m i || represents the Euclidean distance between the i-th group of corresponding points after transformation, min H(m i ) corresponds to the optimal transformation matrix after rough splicing, thus completing the initial splicing.

[0145] S5-4: Use the roughly stitched source point cloud P′ and the target point cloud G as input point clouds for ICP fine stitching, use the IterativeClosestPoint function, create the Best Bin First (BBF) optimized KD tree data structure to accelerate the search for corresponding point pairs; treat all points P′ in the stitched point cloud P′ i , search for the nearest corresponding point G in the target point cloud G i , forming corresponding point pairs; calculating the rotation matrix R and the translation vector T so that the root mean square error M between the corresponding point pairs k Minimum:

[0146]

[0147] The rotation matrices R and T obtained in the first iteration are:

[0148]

[0149] T=[-4.44089e-16 6.66134e-16 -2.98023e-08]

[0150] Finally, set the threshold α (i.e., M k -M k+1 <α) and the maximum number of iterations N max To complete the stitching of multiple point cloud images, so that the overlapping areas between them are perfectly aligned.

[0151] Combination Figure 1 It can be seen that in step S6, since the point cloud shape obtained from the coal mine tunnel arch surface is non-closed point cloud data, a greedy projection triangulation algorithm is used according to its characteristics to realize the three-dimensional surface reconstruction of the coal mine tunnel arch surface. The specific steps are:

[0152] S6-1: Using the recursive algorithm of dynamic programming problem, the normal vector of the three-dimensional point obtained in step 4 is The coordinates of the positive target point O determine the local tangent plane equation of O to obtain the two-dimensional tangent plane formed by adjacent points. If O = (x0, y0, z0), Then the equation of the tangent plane through point O is as follows:

[0153] A(x-x0)+B(y-y0)+C(c-c0)=0

[0154] S6-2: Project the 3D point onto the 2D tangent plane. The projection point is stored in the projection matrix for rotation transformation calculation. The calculation process is:

[0155]

[0156] Among them, T Matrix is the translation transformation matrix:

[0157] Where x0, y0, z0 are the translation amounts of each coordinate axis.

[0158] R x Expressed as a rotation matrix of α degrees around the x-axis:

[0159] R y It is expressed as a rotation matrix of θ degrees around the y-axis:

[0160] Combining the tangent plane equation in step S6-1 with the above equation, we can get any point P(x i ,y i , z i ) on the tangent plane O П Projection on:

[0161]

[0162] S6-3: Figure 7 As shown, the point cloud obtained by projection is triangulated in the plane using the spatial region growing algorithm based on Delaunay. The triangulation process should meet the empty circumscribed circle characteristics and the maximum and minimum angle criteria of the Delaunay algorithm. The algorithm selects a sample triangle as the initial surface, and continuously expands the surface boundary to obtain the connection relationship of each point and form a complete triangular mesh surface. Finally, the topological connection between the original three-dimensional points is determined according to the connection relationship of the projected point cloud. The obtained triangular mesh is the reconstructed surface model.

[0163] Combination Figure 1 It can be seen that after the three-dimensional surface reconstruction is completed in step S6, the computer in step S7 performs online repair on the polygonal voids generated after the reconstruction of the coal mine tunnel arch surface using a void repair algorithm based on a triangular mesh. The specific steps are:

[0164] S7-1: Boundary point detection: The boundary feature points are judged by measuring the distribution uniformity of the k-neighborhood in the point cloud, and the maximum angle difference is used as a measure of distribution uniformity, such as Figure 8 and Fig. 9 shown.

[0165] Definition i = (i = 1, 2, ..., k) is the projection point of the neighborhood point P on the tangent plane, and the nearest neighbor point O of P is taken. i Form a line segment with P by As a benchmark, calculate Rotate clockwise to Angle Make an angle sequence Sort the angles to get Define the angle sequence difference as like Fig.10 shown.

[0166] From L = (L1, L2, ..., L k ) to find the maximum angle sequence difference L max , as the basis for judging the boundary feature points, set the threshold, when L max When it is greater than the threshold, P is a boundary feature point, otherwise P is an internal point;

[0167] S7-2: After finding the boundary feature points, the disordered feature points are ordered, the nearest neighbor points of the feature points are selected as connection points, and the iteration is repeated until a closed boundary line is formed;

[0168] S7-3: After obtaining the boundary line, the hole is repaired online according to the angle α between two adjacent edges in the hole boundary, such as Fig.11 As shown:

[0169] (1) When α≤90°, only one triangular patch is constructed, such as Fig.11 (a)

[0170] (2) When 90°<α≤135°, construct two triangular patches to fill point Q and satisfy P i Q divides ∠P equally i-1 P i P i+1 ,like Fig.11 (b)

[0171] (3) When 135°<α≤200°, three triangular patches are constructed to fill points Q and R, and P should be satisfied. i Q and P i R divides ∠P into three equal parts i-1 P i P i+1 ,like Fig.11 (c) as shown.

[0172] (4) When α>200°, three triangular facets are constructed to fill points Q and R and satisfy ΔQP i+1 P i and ΔPP i R i-1 is an equilateral triangle, such as Fig.11 (d) as shown.

[0173] S7-4: Added a new triangle patch validity check: Fig.12 The characteristics of the hollow polygon are destroyed. Fig.13 The newly constructed triangles intersect with the hole polygons.

[0174] for Fig.12 The test is done by calculating arrive The angle δ1 and arrive The angle δ2 is used to judge. If δ2>δ1, the newly constructed triangle is unreasonable. Fig.13 , by checking whether the newly constructed triangle face intersects with the original hole polygon. If so, the generated triangle face is unreasonable.

[0175] Combination Fig.14It can be seen that the area required for performing a three-dimensional reconstruction of the coal mine tunnel arch surface is 1.6m, so the design of the robot arm base rail length is 1.6m, the installation position of the laser radar is located on the rail, and the hydraulic motor moving device is used to drive the laser radar to move at a constant speed on the rail to collect the point cloud data of the coal mine tunnel arch surface; the experiment uses a 16-line explosion-proof laser radar. Due to the limitation of the laser radar polar line constraint, the vertical field of view angle can only collect point cloud data within the range of (-15°~+15°) at a time, and it is impossible to obtain point clouds within a range of 1.6m at a time. Therefore, multiple sampling points are designed on the rail. The point cloud data of each sampling point are spliced ​​to obtain the complete point cloud within 1.6m of the coal mine tunnel arch surface, and on this basis, the three-dimensional surface reconstruction of the coal mine tunnel arch surface is realized.

[0176] Finally, it should be noted that the above embodiments are only used to illustrate the technical solution of the present invention rather than to limit it. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solution of the present invention can be modified or replaced by equivalents without departing from the purpose and scope of the technical solution of the present invention, which should be included in the scope of the claims of the present invention.

Claims

1. A three-dimensional surface reconstruction method for a coal mine tunnel arch surface, characterized in that: The laser radar is used to scan the arch surface of the coal mine tunnel in real time to obtain point cloud data. The collected point cloud data is reconstructed into a three-dimensional surface according to the following steps: Preprocessing: Flip the laser radar 90° and install it upside down on the guide rail of the robot arm base; collect tunnel arch data and set up multiple data collection stations; use the hydraulic motor moving device to drive the laser radar to move at a constant speed on the guide rail at a certain time interval to obtain complete tunnel arch point cloud data, and set the initial point and moving speed of the laser radar on the guide rail respectively; set the time interval for collecting point cloud data at each measuring station during the movement of the laser radar; S1: Initialization, the laser radar on the guide rail of the robot base moves to the initial point, and the computer and the laser radar establish communication; S2: The guide rail of the robot arm base drives the laser radar collector to move to the data sampling points on the guide rail at a certain speed. The laser collector collects real-time data on the arch surface of the coal mine tunnel at a preset time interval and transmits the collected point cloud data to the computer; S3: The computer analyzes the point cloud data to determine whether there are noise, outliers and invalid points, and uses a filter to pre-process the point cloud data; S4: Based on the results of the data preprocessing in the previous step, the computer uses the principal component analysis method to perform plane fitting on the N points in the K neighborhood of each point to estimate the normal vector and calculate the curvature at the same time; S5: Multiple point cloud images collected from different sampling points on the guide rail are stitched according to the sequence of key point extraction, SAC-IA rough stitching based on FPFH feature descriptor, and ICP fine stitching, so that the intersecting areas between them overlap perfectly; S6: The point cloud shape obtained from the arch surface of the coal mine tunnel is non-closed point cloud data. According to its characteristics, the greedy triangulation algorithm is used to realize the 3D surface reconstruction; S7: After completing the 3D surface reconstruction of the tunnel arch, the polygonal holes generated by diluting part of the point cloud data on the tunnel arch are repaired online using a hole repair algorithm based on a triangular mesh; S8: The guide rail at the base of the robot arm drives the laser radar to move, complete the three-dimensional surface reconstruction of the arch surface of the coal mine tunnel, and return to the initial point. After the grouting is completed, it enters the next arch surface reconstruction.

2. A three-dimensional curved surface reconstruction method for a coal mine tunnel arch surface according to claim 1, characterized in that: In steps S1 and S2, the initial moving point of the guide rail of the robot base is (x0, y0), the moving speed is S, the time interval for the laser radar collector to collect point cloud data is T, and the collector establishes communication between the computer and the laser radar through the IP address; Among them, the guide rail of the base of the robotic arm drives the laser radar collector to collect point cloud data at a certain speed and time interval through a hydraulic motor moving device, and obtains the distance d between different data collection stations and the tunnel arch scanning width w. Multiple point cloud images will be generated within the entire guide rail range; the robotic arm is installed on a mobile platform trolley.

3. The method for three-dimensional surface reconstruction of a coal mine tunnel arch surface according to claim 1, characterized in that: In step S3, the specific steps of point cloud data preprocessing are: S3-1: Use statistical filters to perform statistical analysis on the K neighborhood points of each point in the point cloud, and calculate the average distance from it to all neighboring points; then calculate the average value μ and standard deviation σ of the average distance of each point to determine the distance threshold thresh_d; according to the point cloud density distribution, remove the point cloud whose average neighborhood distance of a point is lower or higher than its threshold range; Assume that the original point cloud dataset P = {p i (x i ,y i , z i )|1≤i≤n}, then the distance threshold is: thresh_d=μ+m·σ The point cloud dataset after outlier filtering is: P′=(um*σ,u+m*σ) Where m is the standard deviation multiple; S3-2: Import P′ into the pass filter, set multiple dimensional directions and different point cloud thresholds, filter out points outside the parameter range, and obtain point cloud data within the specified range; S3-3: Import the point cloud data obtained in S3-2 into a voxel filter or a conditional filter, set a variety of voxel grid thresholds and other conditions of different specifications to simplify the point cloud data, obtain the point cloud preprocessing results, and greatly reduce the amount of computer operations.

4. The method for three-dimensional surface reconstruction of a coal mine tunnel arch surface according to claim 1, characterized in that: In step S4, the specific steps of estimating the point cloud normal using the PCA algorithm are: S4-1: For any point p in the point cloud i (1≤i≤n) query its k-domain point p ij (1≤j≤k), calculate p i The centroid of its k-neighborhood points S4-2: Construct the covariance matrix of local features of point cloud: S4-3: Calculate the eigenvalues ​​λ0, λ1, λ2 of the covariance matrix C, λ0≤λ1≤λ2, and the eigenvectors v0, v1, v2, where λ and v correspond one to one; Among them, v0, v1, v2 are orthogonal, and the eigenvectors v1 and v2 determine the point p i An optimal tangent plane at, v0 is orthogonal to the tangent plane, so point p i Normal at n i It can be approximately represented by the eigenvector v0; λ0, λ1, λ2 are the degree of change in the direction of their respective eigenvectors; S4-4: Calculate point p i The curvature τ pi :

5. The method for three-dimensional surface reconstruction of a coal mine tunnel arch surface according to claim 1, characterized in that: In step S5, the computer completes the point cloud image stitching according to the key point extraction, SAC-IA rough stitching based on FPFH feature descriptor and ICP fine stitching sequence for multiple point cloud images collected from different sampling points on the guide rail. The specific steps are as follows: S5-1: Extract key points of point cloud images using internal morphological descriptor algorithm: Let point cloud P = {p i (x i ,y i , z i )|1≤i≤n}, for each point p i Establish a local coordinate system and set the search radius p r , determine p i is the center of the sphere, p r are all points within a sphere of radius ; Calculate the weight ω ij ; Calculate each point p i The covariance matrix of is: The eigenvalues ​​of the covariance matrix Arrange in descending order; set thresholds η1 and η2 to satisfy and The point is the key point, and iterate until all key points are found; S5-2: Calculate the fast point feature histogram features of the point cloud to be spliced ​​and the target point cloud, and obtain each calculation point M p The relative relationship between all its neighboring points is used to establish a simple point feature histogram; the FPFH feature is calculated based on the SPFH feature, denoted as F(M P ): Among them, d i is the Euclidean distance of corresponding point pairs; S5-3: Use the sampling consistency initial registration algorithm to stitch the model point cloud and the target point cloud. The algorithm steps are as follows: (1) Select m feature points to be registered in the model point cloud; (2) Find points in the target point cloud that are similar to the FPFH features of the model point cloud, and select the points with the closest distance as the corresponding relationship points; (3) Calculate the rotation and translation matrices of the corresponding point pairs; Use the Huber function to represent the distance error and function after the corresponding point pairs are rotated and translated, denoted as Where r is the distance threshold, ||m i || represents the Euclidean distance between the i-th group of corresponding points after transformation, min H(m i ) corresponds to the optimal transformation matrix after rough splicing; S5-4: The roughly stitched source point cloud P′ and the target point cloud G are used as input point clouds for ICP fine stitching, and the optimal node is used to optimize the KD tree first to accelerate the search for corresponding point pairs; all points P′ of the registration point cloud P′ are treated as i , search for the nearest corresponding point G in the target point cloud G i , forming corresponding point pairs; calculating the rotation matrix R and the translation vector T so that the root mean square error M between the corresponding point pairs k Minimum: Finally, set the threshold α and the maximum number of iterations N max , M k -M k+1 <α.

6. The method for three-dimensional surface reconstruction of a coal mine tunnel arch surface according to claim 1, characterized in that: In step S6, a greedy projection triangulation algorithm is used to realize the three-dimensional reconstruction of the arch surface of the coal mine tunnel. The specific steps are as follows: S6-1: Using the recursive algorithm of dynamic programming problem, the normal vector of the three-dimensional point obtained in step 4 is and the coordinates of the target point O to determine the local tangent plane equation of O, so as to obtain the two-dimensional tangent plane formed by adjacent points; if Then the equation of the tangent plane through point O is as follows: A(x-x0)+B(y-y0)+C(c-c0)=0 S6-2: Project the 3D point onto the 2D tangent plane. The projection point is stored in the projection matrix for rotation transformation calculation. The calculation process is: Among them, T Matrix is the translation transformation matrix: Where x0, y0, z0 are the translations of each coordinate axis; R x Expressed as a rotation matrix of α degrees around the x-axis: R y It is expressed as a rotation matrix of θ degrees around the y-axis: Combining the tangent plane equation in step S6-1 with the above equation, we can get any point P(x i ,y i , z i ) on the tangent plane O ∏ Projection on: S6-3: Use the Delaunay-based spatial region growing algorithm to triangulate the projected point cloud in the plane. The triangulation process should satisfy the empty circumscribed circle characteristics and the maximum and minimum angle criteria of the Delaunay algorithm. The algorithm selects a sample triangle as the initial surface and continuously expands the surface boundary to obtain the connection relationship between each point and form a complete triangular mesh surface. Finally, the topological connection between the original three-dimensional points is determined according to the connection relationship of the projected point cloud. The obtained triangular mesh is the reconstructed surface model.

7. The method for three-dimensional surface reconstruction of a coal mine tunnel arch surface according to claim 1, characterized in that: In step S7, the polygonal holes generated by the dilution of part of the point cloud by the arch surface are repaired online using a hole repair algorithm based on a triangular mesh. The specific steps are as follows: S7-1: Perform boundary point detection on the tangent plane. If a point P in the point cloud data is a boundary feature point, then the k-neighborhood points of P should fall on its side; if P is an internal point, then the k-neighborhood points of P should fall around point P; the boundary feature points are judged by measuring the distribution uniformity of the k-neighborhood in the point cloud, and the maximum angle difference is used as a measure of distribution uniformity: Definition i = (i = 1, 2, ..., k) is the projection point of the P neighborhood point on the tangent plane, and the nearest neighbor O of P is taken. i Form a line segment with P by As a benchmark, calculate Rotate clockwise to Angle Make an angle sequence Sort the angles to get a new angle sequence Define the angle sequence difference as From L = (L1, L2, ..., L k ) to find the maximum angle sequence difference L max , as the basis for judging the boundary feature points; setting a threshold, when L max When it is greater than the threshold, P is a boundary feature point, otherwise P is an internal point; S7-2: After finding the boundary feature points, the disordered feature points are ordered, the nearest neighbor points of the feature points are selected as connection points, and the iteration is repeated until a closed boundary line is formed; S7-3: After obtaining the boundary line, the polygonal hole is repaired online according to the angle α between two adjacent edges in the hole boundary; (1) When α≤90°, only one triangular patch is constructed; (2) When 90°<α≤135°, construct two triangular facets to fill point Q and satisfy P i Q divides ∠P equally i-1 P i P i+1 ; (3) When 135°<α≤200°, three triangular patches are constructed to fill points Q and R, and P should be satisfied. i Q and P i R divides ∠P into three equal parts i-1 P i P i+1 ; (4) When α>200°, three triangular facets are constructed to fill points Q and R and satisfy ΔQP i+1 P i and ΔPP i R i-1 is an equilateral triangle; S7-4: For the above hole repair methods, the newly added triangular facets need to be checked for legality; this can be done by calculating arrive The trend of δ1 and arrive It can be judged by the direction δ2. If δ2>δ1, the newly constructed triangle patch is unreasonable. It can be judged by checking whether the newly constructed triangle patch intersects with the original hole polygon. If they intersect, the generated triangle patch is unreasonable.

8. The method for three-dimensional surface reconstruction of a coal mine tunnel arch surface according to claim 1, characterized in that: For the area where a three-dimensional reconstruction of the coal mine tunnel arch surface is performed, the installation position of the laser radar is located on the guide rail. The hydraulic motor moving device is used to drive the laser radar to move at a constant speed on the guide rail to collect the point cloud data of the coal mine tunnel arch surface. After obtaining the point cloud data of the complete coal mine tunnel arch surface, the three-dimensional surface reconstruction is performed.

Images