Blasting area surface morphology inversion method based on unmanned aerial vehicle

Abstract

The present invention belongs to the technical field of digital mine safety, and particularly relates to a blasting area surface morphology inversion method based on an unmanned aerial vehicle. The method comprises: S01, on-site data collection; S02, blasting area morphology inversion; S03, a muck pile throw distance; and S04, a muck pile surface fragmentation distribution. In the present invention, oblique photography by an unmanned aerial vehicle is used, and three-dimensional model reconstruction is performed by capturing a blasting area image, so that a feasible aerial survey scheme is formulated, and the integrity of collected data is good; and reverse modeling of the blasting area is performed by means of the steps of feature point extraction, spatial information conversion, point cloud generation, grid generation, etc.; voxel grid downsampling, point cloud matching and error detection are used to perform registration on point clouds of a target area before and after blasting, so that a good effect is achieved; and a color-based region growing algorithm is used to perform coarse segmentation of point cloud features on an ore-rock block on the surface of a muck pile, and a PointNet++ algorithm is used to perform fine segmentation of point cloud features on the ore-rock block on the surface of the muck pile, so that the muck pile throw distance and the muck pile surface fragmentation distribution are calculated.

Description

Unmanned aerial vehicle-based blasting area surface morphology inversion method TECHNICAL FIELD

[0001] The present application belongs to the technical field of digital mine safety, in particular to an unmanned aerial vehicle-based blasting area surface morphology inversion method. BACKGROUND

[0002] Open-pit mine bench blasting effect analysis is an important link for evaluating and optimizing bench blasting construction, mainly including block size distribution, blast pile height, blast flyrock distance, blast vibration, etc. Among them, the analysis of blast pile morphology is crucial, reflecting the rationality of blasting parameters and charge structure, and directly affecting the shovel loading efficiency and mining area layout. In the case of determined mine excavation equipment selection, to improve the shovel loading efficiency and full bucket coefficient of the excavator, the blast pile morphology is very important. The blast pile should not be too high or too scattered. The upper part of the blast pile will collapse if it exceeds the maximum excavation height of the excavator, and the scattered blast pile will reduce the full bucket coefficient of the excavator. At the same time, analyzing the blast pile morphology and surface characteristics is of great significance for evaluating the blasting effect and fine blasting design. At present, unmanned aerial vehicle aerial survey technology is increasingly widely used in open-pit mines, especially in three-dimensional model reconstruction, safety hazard monitoring and investigation, and engineering measurement, and has good application effect. However, the unmanned aerial vehicle carries three-dimensional laser scanning equipment, and the collected data is easily affected by the environment, resulting in many point cloud noise points, inaccurate collected data, and inaccurate analysis of the surface morphology of the blasting area. The use of unmanned aerial vehicle oblique photography to collect blast pile images and invert the blast pile morphology is relatively weak in research and application.

[0003] Chinese invention patent application No. CN202310597675.8 discloses an open-pit bench blasting blast pile morphology prediction method and system, relating to the field of engineering blasting, which comprises: obtaining the two-dimensional contour curve and morphology parameters of the test blast pile; randomly selecting the blast pile height at the farthest throw distance and the slope of the two-dimensional contour curve to determine the initial Weibull curve shape control parameter; comparing the curve drawn using the above shape parameter with the two-dimensional contour curve to determine the optimal blast pile height at the farthest throw distance and the slope of the two-dimensional contour curve; calculating the throw speed and throw angle of each broken block for all spherical cartridges; calculating the predicted throw distance of each broken block and comparing to obtain the farthest predicted throw distance; estimating the blast pile loose coefficient; obtaining the optimal curve shape control parameter; drawing the optimal Weibull distribution model curve; and using the visible depth of the blasting crater to correct the above curve to obtain the open-pit bench blasting blast pile morphology prediction curve.

[0004] SUMMARY

[0005] The present application aims to provide an unmanned aerial vehicle-based blasting area surface morphology inversion method, which includes blast pile throw distance and surface block size distribution.

[0006] The object of the application is achieved by the technical solutions as follows:

[0007] The unmanned aerial vehicle-based blasting area surface shape inversion method of the application comprises the following steps:

[0008] S01: field data acquisition: in order to obtain the oblique photography data before and after the blasting of the open pit, the application selects an unmanned aerial vehicle as a flight platform, uses an airborne camera to collect image data of the target, and photographs the same area vertically: the overlapping area of every adjacent two pictures must reach more than 60%; the inclination angle: the inclination angle of the camera for photographing the same object is not less than ±15°.

[0009] S0101: specific scheme: when collecting the field image data, the unmanned aerial vehicle is used to perform three-layered nested flight route photographing above the open pit, and the interval height of each layer is 5 m. The flight route of the unmanned aerial vehicle is shown in Fig. 1, and the overlapping degree of the photographed images is generally more than 60% in the lateral direction, and it is more appropriate to be more than 70% in the route overlapping degree. The flight route of the unmanned aerial vehicle should be slightly larger than the aerial survey area to ensure the integrity of the reconstructed three-dimensional model.

[0010] S02: blasting area shape inversion: based on the oblique photography, the image of the blasting area is obtained, a group of corresponding images of the characteristic points are selected from the multiple images in different directions, the corresponding two-dimensional characteristic points are obtained, and the position of the modeling area in the real world is obtained by using the space information conversion method. As shown in Fig. 2, then the characteristic points of each image are obtained by using the image index position, and the three-dimensional point cloud data is restored. The steps are as follows:

[0011] S0201: characteristic point extraction. The characteristic points such as edges and corner points in the image are analyzed, the key features of the image are extracted, and the extracted characteristic points are used for subsequent image matching and three-dimensional reconstruction;

[0012] S020101: feature matching. The SIFT algorithm (Scale-invariant Feature Transform) is used for feature point detection and matching.

[0013] S020102: feature point detection. The target pixel points are compared with the pixel points in the corresponding pixel regions on the left, right, top and bottom of the target pixel points, and the extreme points in the regions are taken as the characteristic points.

[0014] S020103: Feature point direction judgment. According to the coordinates of each feature point, one or more directions are assigned to the feature point. And with the feature point as the center, all directions within the surrounding range are compared to construct a gradient histogram. The maximum value in each gradient is taken as the main direction, and when the remaining direction values exceed 80% of the main direction, they are taken as auxiliary directions, as shown in FIG. 2.

[0015] S0202: Spatial information conversion: Using a spatial position transformation method, the basic matrix and matching formula are calculated according to the matching principle of each feature point and feature point pair, and the position information of the feature points in the real space is obtained, as shown in FIGS. 3 and 4.

[0016] S020201: Pixel coordinate to image coordinate: The pixel coordinate system is a two-dimensional rectangular coordinate system that reflects the arrangement of pixels in the image captured by the camera.

[0017] S020202: Image coordinate to camera coordinate: The blast pile in the blasting area is converted from the image coordinate system to the camera coordinate system.

[0018] S020203: Camera coordinate to world coordinate: Camera calibration is to convert the coordinate information of the blast pile in the real world to a coordinate system based on the camera's perspective.

[0019] S020204: Error detection: To avoid insufficient reprojection accuracy between images before determining the reference coordinate system.

[0020] S0203: Point cloud generation: After completing image matching, a dense point cloud is generated based on the spatial position of the feature points in the image, as shown in FIG. 5.

[0021] S0204: Mesh generation: A mesh model is generated by triangulating the point cloud. The mesh model is formed by many small triangles, each with its vertex and normal direction.

[0022] S0205: Texture mapping: To make the generated model more realistic and highlight the surface rock block distribution characteristics of the blast pile, the texture in the image is mapped to the surface of the mesh model, as shown in FIG. 6.

[0023] S03: Blast pile throw distance: The blast pile throw distance is calculated by point cloud registration, blast pile point cloud segmentation, and blast pile throw distance calculation method.

[0024] S0301: Point cloud registration: Before and after the blasting of the blast area is modeled, affected by the shooting environment, the flight state of the unmanned aerial vehicle and the picture collection quality, the established before and after the blasting of the blast area model often has different coordinate systems. The calculation of the throwing distance of the blast pile needs to integrate the before and after blasting point cloud models in different coordinate systems in the same coordinate system, and the core is to calculate the rotation matrix R and the translation matrix T of the two pieces of point cloud by using the downsampling and registration method. The before and after blasting point clouds are shown in Figure 7-a.

[0025] S030101: Downsampling: The before and after blasting point clouds are processed by downsampling, and ten thousand points are downsampled to several hundred to improve the speed of the algorithm. The downsampling method of the application: the point cloud data is divided into a plurality of cubic voxels with the same size, and the nearest point to the center of the voxel is selected as the sampling result.

[0026] S030102: Point cloud matching: Assuming that the point cloud {Q} is the target point cloud (reference point cloud), {P} is the source point cloud (to-be-registered point cloud), p i (i∈1,2,...N) is a point in {S}, q i is the nearest point in {E} to p i . The transformation matrix from {P} to {Q}, i.e., the rotation matrix R and the translation matrix T, is calculated. If the transformation parameters are accurate, each point p i in the point cloud {P} should coincide completely with the point q i in the point cloud {Q} after transformation, i.e.: q i =Rp i +T. However, due to the existence of noise, it is impossible for all points to coincide completely, so the objective function is defined as:

[0027] The R and T that minimize the objective function E are the transformation parameters to be solved. F is actually the average distance between the reference point cloud {Q} and {P'} that has been subjected to the R and T matrix space transformation.

[0028] S030103: Error detection: For each point p i in {P}, find its nearest point q i in {Q} to form a one-to-one point pair. Calculate the centroids of the two groups of point clouds, denoted as u p and u q :

[0029] The centroids of the two groups of point clouds are removed to obtain:

[0030] p′ i =p i -u p , q′ i =qi -u q (13)

[0031] Construct the matrix H:

[0032] SVD decomposition is performed on the H matrix:

[0033] H=UΣV T (15)

[0034] R and T are obtained:

[0035] R=VU T , T=u q -Ru p (16)

[0036] After obtaining the R and T matrices, a new point set is obtained by performing a spatial transformation on the to-be-registered space, and is substituted into the objective function:

[0037] If the new transformed point set and the reference point set satisfy the average distance between the two point sets being less than a given threshold, the iterative calculation is stopped, otherwise the new transformed point set is taken as a new {P i} to continue iteration until the requirements of the objective function are met. The registration result is shown in Figure 7-b.

[0038] S0302: Blasting pile point cloud segmentation: Point cloud segmentation is to divide points according to spatial, geometric and texture features, and points in the same division have similar features. The purpose of point cloud segmentation is to block, thereby facilitating individual processing.

[0039] The point cloud segmentation of the present application adopts a point cloud segmentation method based on geometric features, sets the slice thickness and the table thickness, and cuts the point cloud. The present application performs equidistant slicing processing on the point cloud data, which is convenient for subsequent calculation of blasting pile morphology throwing distance, and the point cloud segmentation effect is shown in Figure 8.

[0040] S0303: Throwing distance calculation: Through the above S0302, the registered point cloud is segmented into multiple segments, and for each segmented table body, the blasting pile throwing distance is calculated using the step feature. By counting the coordinates of the extreme points of the point cloud before and after blasting, the blasting pile throwing distance is calculated, as shown in Figure 9, the step height is H, the blasting pile height is h, the blasting step width is B, the mine rock blasting pile extension distance is b, and the blasting pile throwing distance is b0,

[0041] b0=b-(B0-B) (18)

[0042] S04: Blast pile surface block size distribution: The surface block size calculation uses the point cloud model of the blast area after blasting with texture and color features. The target is detected through two steps of point cloud feature coarse segmentation and block size fine segmentation, and the surface block size distribution of the blast pile is calculated.

[0043] S0401: Point cloud feature coarse segmentation: A color region growing algorithm is used for preliminary screening to determine the suspicious area. According to the same color (similar color) and close distance, it is highly likely that they belong to the same target, and the continuous scene point cloud is changed into different objects.

[0044] The algorithm mainly consists of two steps:

[0045] (1) Segmentation: The color difference between the current seed point and the field point is less than the color difference threshold, which is considered as a cluster.

[0046] (2) Merge: The color difference between clusters is less than the color difference threshold, and the number of points in the current cluster is less than the number of points in the nearest cluster, which are merged together.

[0047] Distance calculation of RGB:

[0048] The preliminary screening results are shown in Figure 10.

[0049] S0402: Block size fine segmentation: PointNet++ algorithm is used for blast pile block size segmentation to identify and segment the surface large blocks of the blast pile.

[0050] Combined with the production situation of the mine, the volume threshold is set, and the blocks larger than the threshold are considered as large blocks. The PointNet++ point cloud segmentation process is shown in Figure 11.

[0051] PointNet++ uses the concept of neighborhood to effectively model local features, extract local features at different scales, and obtain deep features through a multi-layer network structure. The PointNet++ point cloud segmentation process is shown in Figure 11, and the block size segmentation results are shown in Figure 12.

[0052] S0403: Block size distribution calculation: Based on the above blast pile surface block size segmentation results, it can be known that the left surface of the blast area has the most rock blocks, the right area of the blast area has the second most rock blocks, and the middle area of the blast area has the least rock blocks. From the above results, it can be known that the left side of the blast area needs to be optimized, as shown in Figure 13.

[0053] Large block quantity statistics method:

[0054] if x 大块i <x 左 num 左 +1

[0055] Else if x 左 <x 大块i <x 右 num 中 +1

[0056] Else x 右 <x 大块i num 右 +1

[0057] (i=1、2、3、…n)

[0058] Wherein, x 大块i : x-axis coordinate of the i-th large block, x 左 : left region boundary, x 右 : right region boundary, num: number of different region large blocks.

[0059] The key point of the present application

[0060] The key point of the present application is:

[0061] 1. The present application aims at the problems of low efficiency, poor accuracy, weak real-time performance and insufficient analysis capability in the process of on-site data collection in mine blasting production, which still adopts manual measurement. A three-layer nested cultivation flight line shooting scheme is designed above the stope by using unmanned aerial oblique photography, a feasible photogrammetry scheme is formulated, the data collected has good integrity, and the steps of feature point extraction, spatial information conversion, point cloud generation and grid generation are performed for reverse modeling of the blasting area; voxel grid downsampling, point cloud matching and error detection are used to register the point clouds of the target area before and after blasting, and then the blasting heap point cloud is segmented to calculate the blasting throw distance, which has good effect; the color region growing algorithm is used for coarse segmentation of the point cloud features of the blasting heap surface rock mass, and the PointNet++ algorithm is used for fine segmentation of the point cloud features of the blasting heap surface rock mass to realize the calculation of the blasting heap surface block size distribution. The present application uses two indexes of blasting throw distance and blasting heap surface block size distribution to analyze the blasting heap form, the blasting heap form is the external manifestation of the blasting effect of the open-pit mine, and comprehensively reflects the rationality of blasting design and blasting construction, which provides help for blasting design optimization.

[0062] Advantages of the present application:

[0063] (1) The blasting area surface form inversion method proposed in the present application uses unmanned aerial oblique photography to perform blasting heap three-dimensional model inversion through shooting of blasting area images. It can effectively solve the problems of many noise points, poor effect and high cost of three-dimensional laser scanning of blasting heap, thereby greatly improving the efficiency of blasting heap form inversion;

[0064] (2) The point cloud registration, segmentation and throwing distance calculation method suitable for open-air blast pile form analysis provided by the present application can effectively solve the key problem that the throwing distance of the blast pile is difficult to quickly, safely and accurately collect, thereby greatly improving the measurement efficiency of the key indicators of the blasting effect, and providing important technical support for realizing the blasting effect acceptance of the open-pit mine based on the unmanned aerial vehicle oblique photography;

[0065] (3) The "double feature" extraction method of the blast pile surface block degree based on the point cloud data provided by the present application solves the key problem that the point cloud feature is difficult to quickly extract, greatly improves the feature extraction speed of the blast pile surface block degree, simultaneously, solves the statistical precision problem of the surface block degree based on image recognition by using the point cloud to statistically count the surface block degree of the blast pile, improves the efficiency of the blasting effect acceptance measurement, and provides important technical support for realizing the blasting effect acceptance of the open-pit mine based on the unmanned aerial vehicle oblique photography. BRIEF DESCRIPTION OF DRAWINGS

[0066] Fig. 1 is a flight scheme of an unmanned aerial vehicle based on a blasting area surface form inversion method provided by an embodiment of the present application.

[0067] Fig. 2 is a feature point direction descriptor diagram of an unmanned aerial vehicle based on a blasting area surface form inversion method provided by an embodiment of the present application.

[0068] Fig. 3 is feature point position information of an unmanned aerial vehicle based on a blasting area surface form inversion method provided by an embodiment of the present application.

[0069] Fig. 4 is a camera imaging principle of an unmanned aerial vehicle based on a blasting area surface form inversion method provided by an embodiment of the present application.

[0070] Fig. 5 is a pre-blasting and post-blasting point cloud diagram of an unmanned aerial vehicle based on a blasting area surface form inversion method provided by an embodiment of the present application.

[0071] Fig. 6 is a point cloud diagram with texture of an unmanned aerial vehicle based on a blasting area surface form inversion method provided by an embodiment of the present application.

[0072] Fig. 7 is a pre-registration and post-registration diagram of an unmanned aerial vehicle based on a blasting area surface form inversion method provided by an embodiment of the present application.

[0073] Fig. 8 is a point cloud segmentation effect and segmentation body schematic diagram of an unmanned aerial vehicle based on a blasting area surface form inversion method provided by an embodiment of the present application.

[0074] Fig. 9 is a blast pile throwing profile diagram of an unmanned aerial vehicle based on a blasting area surface form inversion method provided by an embodiment of the present application.

[0075] Fig. 10 is a preliminary screening result of an unmanned aerial vehicle based on a blasting area surface form inversion method provided by an embodiment of the present application.

[0076] Fig. 11 is a PointNet++ point cloud segmentation process of a blasting area surface morphology inversion method based on a UAV according to an embodiment of the present application.

[0077] Fig. 12 is a block segmentation result diagram of a blasting area surface morphology inversion method based on a UAV according to an embodiment of the present application.

[0078] Fig. 13 is a blasting heap surface morphology analysis diagram of a blasting area surface morphology inversion method based on a UAV according to an embodiment of the present application.

[0079] Fig. 14 is a flowchart of a blasting area surface morphology inversion method based on a UAV according to an embodiment of the present application. DETAILED DESCRIPTION

[0080] The specific embodiments of the present application will be further described below in conjunction with the accompanying drawings.

[0081] The UAV-based blasting heap morphology analysis system and method of the present application comprises the following steps:

[0082] S01: On-site data acquisition: In order to obtain the oblique photography data before and after the blasting of the open pit, the present application selects a UAV as the flight platform and uses the on-board camera to collect the image data of the target. The same area is photographed vertically: the overlapping area of every two adjacent pictures must reach more than 60%; the inclination angle: the camera inclination angle for photographing the same object is not less than ±15°.

[0083] S0101: Specific scheme: when collecting on-site image data, the UAV is used to fly three layers of nested plowing routes above the pit, and the interval height of each layer is 5m. The flight route is shown in Fig. 1. The overlapping degree of the photographed images is generally more than 60% in the lateral direction, and it is more appropriate to be more than 70% in the route overlapping degree. The route of the UAV flight should be slightly larger than the surveying area to ensure the integrity of the reconstructed three-dimensional model

[0084] S02: Blasting area morphology inversion: based on the oblique photography to obtain the blasting area image, a group of feature corresponding images are selected from several different direction images, the corresponding two-dimensional feature points are obtained, and the spatial information conversion method is used to obtain the position of the modeling area in the real world. As shown in Fig. 2, then the feature points of each image are obtained by using the image index position, and the three-dimensional point cloud data is restored. The steps are as follows:

[0085] S0201: Feature extraction: analyze the feature points in the image, such as edges, corner points, etc., and extract the key features of the image. The extracted feature points are used for subsequent image matching and three-dimensional reconstruction;

[0086] S020101: Feature matching. SIFT algorithm (Scale-invariant Feature Transform) is used for feature point detection and matching. The theoretical basis of SIFT algorithm is as follows: first, a Gaussian difference pyramid is established, a Gaussian function G(x, y, σ) with a variance of σ is convolved with an image L(x, y) which has been calculated by various methods in advance to obtain the scale space L(x, y, σ) of the image, as shown in formula (1):

[0087] L(x, y, σ) = G(x, y, σ) * L(x, y) (1)

[0088] wherein

[0089] In order to efficiently extract stable feature points of the scale space, the Gaussian difference function is convolved with the preprocessed image to construct a Gaussian difference pyramid, and the formula is as follows:

[0090] D(x, y, σ) = [G(x, y, kσ) - G(x, y, σ)] * L(x, y) = L(x, y, kσ) - L(x, y, σ) (2)

[0091] S020102: Feature point detection: each pixel point in each layer of the Gaussian difference pyramid is compared with the points in the corresponding pixel region on the left, right, top and bottom, respectively, and the extreme points in the region are found as feature points.

[0092] S020103: Feature point direction judgment: according to the coordinates of each feature point, one or more directions are assigned to the feature point. For an image with pixels L(x, y), gradient modulus m(x, y) and direction angle θ(x, y) are used to distinguish different feature vectors:

[0093] S0202: Spatial information conversion: using the principle of matching of each feature point and feature point pair, the basic matrix and matching formula are calculated to obtain the spatial position information of the feature points, as shown in FIG. 3, wherein (x w ,y w ,z w ) is the world coordinate system, (x c ,y c ,z c ) is the camera coordinate system, (x, y) is the image coordinate system, and (u, v) is the pixel coordinate system, as shown in FIG. 4. The position information of the feature points in space is obtained in combination with the camera calibration principle, and the position information conversion process is as follows:

[0094] S020201: Pixel coordinate to image coordinate: Pixel coordinate system is a two-dimensional rectangular coordinate system, which reflects the arrangement of pixels in the image taken by the camera. The conversion formula between pixel coordinate system and image coordinate is as shown in formula (5):

[0095] The above formula is arranged, that is, the target object is converted from the three-dimensional world coordinate system to the pixel coordinate system:

[0096] Through formula (6), the internal and external parameters of the camera can be obtained, and the real blast pile in the modeling area can be uniquely projected into the image.

[0097] S020202: Image coordinate to camera coordinate: Convert the blast pile in the blasting area from the image coordinate system to the camera coordinate system.

[0098] Further, the coordinates of the target object in the pixel coordinate system are obtained.

[0099] Where f is the focal length, x, y is the image coordinate of the target object.

[0100] S020203: Camera coordinate to world coordinate: Camera calibration is to convert the coordinate information of the blast pile in the real world into a coordinate system with the camera view as the reference. This conversion method is divided into translation and rotation. The translation is represented by the translation vector T, and the rotation is represented by the rotation matrix R. The translation vector only needs to add the corresponding translation amount to the object coordinates in the world coordinate system. The rotation matrix needs to consider the rotation around the x, y, and z axes respectively. The public demonstration is shown in formula (8):

[0101] Where R is a 3*3 rotation matrix, T is a 3*1 translation vector, is the homogeneous coordinate of the camera coordinate system, is the homogeneous coordinate of the world coordinate system. In this way, the conversion between different coordinate systems is realized.

[0102] S020204: Error detection: In order to avoid the insufficient reprojection accuracy between images before determining the reference coordinate system, a basic structure is established using two images as reference. If the measured noise satisfies the Gaussian distribution, the distance between the projection point and the image is estimated using the maximum likelihood function, so as to minimize the reprojection error. The objective function is shown in formula (9):

[0103] Where, P i is the projection matrix of the i-th image, X i is the j-th matched feature point, x ij is the j-th point coordinate of the i-th camera, x ijis the homogeneous coordinate of the jth point of the ith image in the projection space. m is the number of images, and n is the number of feature points in the images.

[0104] S0203: Point cloud generation: After completing image matching, a dense point cloud is generated based on the feature points in the images. The generated point cloud is composed of a large number of three-dimensional points, each of which has its position in the camera coordinate system, as shown in FIG. 5.

[0105] S0204: Mesh generation: A mesh model is generated by triangulating the point cloud. The mesh model is formed by a large number of small triangles, each of which has its vertex and normal direction.

[0106] S0205: Texture mapping: In order to make the generated model more realistic and highlight the surface rock block distribution characteristics of the blast pile, the texture in the image is mapped onto the surface of the mesh model, as shown in FIG. 6.

[0107] S03: Blast pile throw distance: The blast pile throw distance is calculated by point cloud registration, blast pile point cloud segmentation, and blast pile throw distance calculation method.

[0108] S0301: Point cloud registration: Based on the pre-blasting and post-blasting modeling of the blast area in the open-pit mine by unmanned aerial photography, the pre-blasting and post-blasting models of the blast area are often established with different coordinate systems due to the influence of the shooting environment, the flight state of the unmanned aerial vehicle, and the quality of the picture collection. The calculation of the throw distance of the blast pile requires the integration of the pre-blasting and post-blasting point cloud models in different coordinate systems into the same coordinate system, and the core calculation is to solve the rotation matrix R and the translation matrix T between the two pieces of point cloud. The pre-blasting and post-blasting point clouds are shown in FIG. 7-a.

[0109] S030101: Down-sampling: The pre-blasting and post-blasting point clouds are down-sampled to several hundred points from tens of thousands of points to improve the speed of the algorithm. The down-sampling method of the present application is to divide the point cloud data into a plurality of cubic voxels of the same size, and select the nearest point to the center of the voxel as the sampling result.

[0110] S030102: Point cloud matching: Assuming that the point cloud {Q} is the target point cloud (reference point cloud), {P} is the source point cloud (point cloud to be registered), p i (i∈1,2,...N) is a point in {S}, q i is the nearest point in {E} to p i . The RT transformation matrix from {P} to {Q}, i.e., the rotation matrix R and the translation matrix T, is calculated. If the transformation parameters are accurate, each point p i in the point cloud {P} should coincide completely with the point q i in the point cloud {Q} after transformation, i.e.: q i =Rp i+T. But due to the existence of noise, it is impossible for all points to be completely coincident, so the objective function is defined as:

[0111] R, T that makes the objective function F minimum is the transformation parameter we want. F is actually the average distance between the reference point cloud {Q} and {P'} that has been spatially transformed by R, T matrix.

[0112] S030103: Error detection: for each point p in {P} i , find its nearest point q in {Q} i , form a one-to-one point pair p , u q :

[0113] De-center the two point clouds, get:

[0114] p′ i =p i -u p , q′ i =q i -u q (13)

[0115] Construct matrix H:

[0116] SVD decompose H matrix:

[0117] H=UΣV T (15)

[0118] Get R and T:

[0119] R=VU T , T=u q -Ru p (16)

[0120] After getting R and T matrix, use it to spatially transform the space to be registered to get a new point set, and substitute it into the objective function:

[0121] If the new transformed point set and the reference point set satisfy the average distance between the two point sets less than a given threshold, stop the iterative calculation, otherwise the new transformed point set is used as the new {P i} to continue iteration until the requirement of the objective function is met. The registration result is shown in Figure 7-b.

[0122] S0302: Point cloud segmentation of the blasting pile: point cloud segmentation is to divide points according to spatial, geometric and texture features, and points in the same division have similar features. The purpose of point cloud segmentation is to divide blocks, thereby facilitating individual processing.

[0123] The point cloud segmentation of the present application adopts a point cloud segmentation method based on geometric features, sets the slice thickness and the table body thickness, and cuts the point cloud. The present application performs equidistant slicing processing on the point cloud data, which facilitates subsequent calculation of the blasting pile morphology throwing distance, and the point cloud segmentation effect is shown in FIG. 8.

[0124] S0303: Throwing distance calculation: through the above S0302, the registered point cloud is segmented into multiple segments, and for each segmented table body, the maximum throwing distance of the blasting pile is calculated using the step feature; the point cloud after segmentation is used for calculation, and for the horizontal of the upper and lower step table surfaces, the position where the two point clouds start to separate is the starting position of the blasting area, which is recorded at this time and intercepted, and the points obtained when calculating the maximum throwing distance of the blasting pile are all within this range; by counting the extreme point coordinates of the point cloud before and after blasting, the throwing distance of the blasting pile is calculated, as shown in FIG. 9, the step height is H, the blasting pile height is h, the blasting step width is B, the mine rock blasting pile extension distance is b, and the blasting pile throwing distance is b0, wherein the throwing distance:

[0125] b0=b-(B0-B) (18)

[0126] S04: Surface block size distribution of the blasting pile: the surface block size calculation adopts a point cloud model of the blasting area after blasting with texture and color features. Through point cloud feature rough division and block size fine segmentation, the block size distribution is calculated based on the surface block size segmentation result of the blasting pile.

[0127] S0401: Point cloud feature rough division: a color region growing algorithm is used for preliminary screening to determine the suspicious area (this method runs fast, but the segmentation effect is not ideal, and only the color features are used to circumscribe the area). The same color (similar color) and close distance have a high possibility of being a class of targets, and the continuous scene point cloud is changed into different objects.

[0128] The algorithm mainly includes two steps:

[0129] (1) Segmentation: the color difference between the current seed point and the field point is less than the color difference threshold, which is regarded as a cluster.

[0130] (2) Merge: the color difference between the clusters is less than the color difference threshold, and the number of points in the current cluster is less than the number of cluster points, and the nearest cluster is merged together.

[0131] Distance calculation of RGB:

[0132] The preliminary screening of the blast pile surface block size distribution characteristics results, using PointNet++ algorithm for block size detailed division. The preliminary screening results are shown in Figure 10.

[0133] S0402: Block size fine segmentation: For the blast pile block size segmentation, PointNet++ algorithm is used to identify and segment the large blocks on the surface of the blast pile.

[0134] Combined with the mine production situation, set the volume threshold, which is greater than the threshold is a large block. The PointNet++ point cloud segmentation process is shown in Figure 11.

[0135] PointNet++ uses the concept of neighborhood to effectively model local features, extract local features at different scales, and obtain deep features through multi-layer network structure. The PointNet++ point cloud segmentation process is shown in Figure 11, and the block size segmentation results are shown in Figure 12.

[0136] S0403: Block size distribution calculation: Based on the above blast pile surface block size segmentation results, it can be known that the left surface of the blast area has the most rock blocks, the right area of the blast area has the second most rock blocks, and the middle area of the blast area has the least rock blocks. From the above results, it can be known that the left side of the blast area needs to be optimized, as shown in Figure 13.

[0137] Large block quantity statistics method:

[0138] if x 大块i <x 左 num 左 +1

[0139] Else if x 左 <x 大块i <x 右 num 中 +1

[0140] Else x 右 <x 大块i num 右 +1

[0141] (i=1、2、3、…n)

[0142] Where, x 大块i : the x-axis coordinate of the i-th large block, x 左 : the left side of the region limit, x 右 : the right side of the region limit, num: the number of large blocks in different regions.

Claims

1. A method for inversion of surface morphology of a blast area based on a drone, characterized in that, It comprises the following steps: S01: field data acquisition: in order to obtain the tilt photography data before and after the open pit blasting, the unmanned aerial vehicle is used as a carrying platform and a flight platform, and the image data of the target is collected by using the airborne camera, the same area is shot, the vertical shooting: the overlapping area of every adjacent two pictures reaches more than 60%; the tilt angle: the camera tilt shooting angle of the same object is not less than +15°; S02: blast area shape inversion: based on the tilt photography, the blast area image is obtained, a group of feature corresponding images are selected from multiple images in different directions, the corresponding two-dimensional feature points are obtained, the spatial information conversion method is used to obtain the position of the modeling area in the real world, the feature points of each image are obtained by using the image index position, and the three-dimensional point cloud data is restored; S03: blast pile throwing distance: the blast pile throwing distance is calculated by point cloud registration, blast pile point cloud segmentation and blast pile throwing distance calculation method; S04: blast pile surface block degree distribution: the surface block degree calculation adopts the blast area point cloud model with texture and color characteristics; the block degree distribution is calculated based on the blast pile surface block degree segmentation result through point cloud feature rough segmentation and block degree fine segmentation.

2. The UAV-based blast zone surface morphology inversion method of claim 1, wherein In the S01, the steps are as follows: S0101: when collecting the image data on site, the unmanned aerial vehicle is used to shoot the three-layer ploughing flight path above the mining field, and the interval height of each layer is 5m; the ploughing flight path usually divides the given flight area into two areas, and the first flight path in the first area is entered from the first area, after flying the first flight path, the first flight path in the second area is entered by turning, and the second flight path in the first area is returned, and the flight path in the second area is entered, and the flight path in the two areas is flown in this way, and the overlapping degree of the shot image is generally more than 60% in the lateral direction, and the overlapping degree of the flight path is more than 70%; when the unmanned aerial vehicle flies, in order to ensure that the complete three-dimensional model is reconstructed, the flight path of the unmanned aerial vehicle is greater than the flight survey area.

3. The UAV-based blast zone surface morphology inversion method of claim 1, wherein In the S02, the steps are as follows: S0201: feature extraction: the edge, corner feature point and image key feature extraction are analyzed, and the extracted feature points are used for subsequent image matching and three-dimensional reconstruction; S020101: feature matching: the SIFT algorithm is used for feature point detection and matching; the theoretical basis of the SIFT algorithm is as follows: first, a Gaussian difference pyramid is established, a Gaussian function G(x,y,σ) with a variance σ is convolved with an image L(x,y) which has been calculated by various methods in advance, and a scale space L(x,y,σ) corresponding to the image is obtained, as shown in formula (1): L(x,y,σ)=G(x,y,σ)*L(x,y) (1) wherein In order to efficiently extract the stable feature points of the scale space, the Gaussian difference function and the preprocessed image are convolved to construct a Gaussian difference pyramid, and the formula is as follows: D(x,y,σ)=[G(x,y,kσ)-G(x,y,σ)]*L(x,y)=L(x,y,kσ)-L(x,y,σ)(2) S020102: Feature point detection: compare each pixel in each layer of the Gaussian difference pyramid with the corresponding pixels in the left, right, top and bottom regions, respectively, to find the extreme points in the region as feature points; S020103: Feature point direction judgment: according to the coordinates of each feature point, one or more directions are assigned to the feature point; for an image with pixels L(x, y), the gradient modulus m(x, y) and the direction angle θ(x, y) are taken to distinguish different feature vectors: S0202: Spatial information conversion: using triangulation method, according to the matching principle of each feature point and feature point pair, the basic matrix and matching formula are calculated to obtain the spatial position information of the feature points, wherein (x w ,y w ,z w ) is the world coordinate system, (x c ,y c ,z c ) is the camera coordinate system, (x,y) is the image coordinate system, and (u,v) is the pixel coordinate system; the position information of the feature points in space is obtained by combining the camera calibration principle, and the position information conversion process is as follows: S020201: Pixel coordinate to image coordinate: The pixel coordinate system is a two-dimensional rectangular coordinate system, which reflects the arrangement of pixels in the image taken by the camera; the conversion formula of the pixel coordinate system and the image coordinate is as formula (5): The above formula is arranged, i.e. the target object is converted from a three-dimensional world coordinate system to a pixel coordinate system: The internal and external parameters of the camera can be obtained through formula (6), and the real blast pile in the modeling region can be uniquely projected into the image; S020202: Image to camera coordinate conversion: convert the blast pile of the blast area from the image coordinate system to the camera coordinate system; Where f is the focal length, x and y are the image coordinates of the target object; S020203: Camera coordinate to world coordinate: camera calibration is to convert the coordinate information of the blast pile in the blasting area in reality into a coordinate system based on the camera view angle. This conversion method includes translation and rotation. Translation is represented by translation vector T, and rotation is represented by rotation matrix R; The translation vector only needs to add the corresponding translation to the object coordinate in the world coordinate system; the rotation matrix needs to consider the rotation around the x, y, and z axes respectively, as shown in formula (8): where R is a 3*3 rotation matrix and T is a 3*1 translation vector, homogeneous coordinates of the camera coordinate system, The world coordinate system is homogeneous coordinate; thus, the conversion between different coordinate systems is realized; S020204: error detection: a basic structure is established using two images as a reference, if the measured noise satisfies Gaussian distribution, the distance between the projection point and the image is estimated using maximum likelihood function, so as to minimize the re-projection error, the objective function is shown in equation (9): where P i is the projection matrix of the i-th image, X i is the j-th matched feature point, x ij is the coordinate of the j-th point of the i-th camera, x ij is the homogeneous coordinate of the j-th point of the i-th image in the projection space. M is the number of images, and n is the number of feature points in the image; S0203: Point cloud generation: after image matching, a dense point cloud is generated according to the feature points in the image; the generated point cloud is composed of a large number of three-dimensional points, and each point has its position in the camera coordinate system; S0204: Mesh generation: a mesh model is generated by triangulating the point cloud; the mesh model is formed by many small triangles, and each triangle has its vertex and normal direction; S0205: Texture mapping: the texture in the image is mapped to the surface of the mesh model.

4. [Amended according to Rule 26 06.08.2024] The unmanned aerial vehicle-based blast zone surface morphology inversion method according to claim 1, characterized in that In S03, the steps are as follows: S0301: Point cloud registration: based on the pre-blasting and post-blasting modeling of the open-pit mine blast area by unmanned aerial photography, the pre-blasting and post-blasting models of the blast area are often established with different coordinate systems due to the influence of the shooting environment, the flight state of the unmanned aerial vehicle and the quality of the picture collection. The calculation of the throw distance of the blast pile requires integrating the pre-blasting and post-blasting point cloud models in different coordinate systems into the same coordinate system, and the core calculation is to solve the rotation matrix R and the translation matrix T between the two point clouds; S030101: Down-sampling: the pre-blasting and post-blasting point clouds are down-sampled to several hundred points from tens of thousands of points to improve the speed of the algorithm; the down-sampling method of the present application is to divide the point cloud data into a plurality of cubic voxels with the same size, and select the nearest point to the center of the voxel as the sampling result; S030102: Point cloud matching: Assume point cloud {Q} is the target point cloud (reference point cloud), {P} is the source point cloud (point cloud to be registered), p i is a point in {S}, q i is the point in {E} closest to p i ; Compute the RT transformation matrix from {P} to {Q}, i.e. the rotation matrix R and the translation matrix T; if the transformation parameters are accurate, every point p in the point cloud {P} i , after transformation, should coincide exactly with a point q i in the point cloud {Q}, i.e.: i q i = Rp + T, but due to the presence of noise, it is not possible for all points to coincide exactly, so a target function is defined: R and T that minimize the target function F are the transformation parameters to be solved; F is actually the average distance between the reference point cloud {Q} and {P'} that has been spatially transformed by R and T matrix; S030103: error detection: for each point p in {P} i , find its nearest point q in {Q} i , form a one-to-one point pair p , compute the centroids of the two point clouds, respectively denoted as u q : The two sets of point clouds are decentered to obtain: p' i = p i - u p , q' i = q i - u q (13) Constructing the matrix H: SVD decomposition of H matrix: H = U∑V T (15) R and T are obtained: R = VU T , T = u q -Ru p (16) After obtaining the R, T matrix, the new point set is obtained by spatial transformation of the to-be-registered space using the R, T matrix, and substituted into the objective function: If the average distance between the new transformed point set and the reference point set satisfies a given threshold, the iterative calculation is stopped, otherwise the new transformed point set is taken as the new {Pi} for continuous iteration until the requirement of the target function is met; S0302: Blast pile point cloud segmentation: point cloud segmentation is to divide the point cloud according to spatial, geometric and texture features, and the point clouds in the same division have similar features; The purpose of point cloud segmentation is to divide the point cloud into blocks for separate processing; The point cloud segmentation method based on geometric features is adopted, the slice thickness and the table thickness are set, and the point cloud is cut; The point cloud data is processed by equidistant slicing, which is convenient for subsequent analysis of the shape of the blast pile; In S03, the steps are as follows: S0301: Point cloud registration: based on the pre-blasting and post-blasting modeling of the open-pit mine blast area by unmanned aerial photography, the pre-blasting and post-blasting models of the blast area are often established with different coordinate systems due to the influence of the shooting environment, the flight state of the unmanned aerial vehicle and the quality of the picture collection. The calculation of the throw distance of the blast pile requires integrating the pre-blasting and post-blasting point cloud models in different coordinate systems into the same coordinate system, and the core calculation is to solve the rotation matrix R and the translation matrix T between the two point clouds; S030101: Down-sampling: the pre-blasting and post-blasting point clouds are down-sampled to several hundred points from tens of thousands of points to improve the speed of the algorithm; the down-sampling method of the present application is to divide the point cloud data into a plurality of cubic voxels with the same size, and select the nearest point to the center of the voxel as the sampling result; S0303: Calculate the throwing distance: through the above S0302, the point cloud after registration is segmented into multiple segments, and the maximum throwing distance of the blast pile is calculated for each segmented bench. The point cloud after segmentation is used for calculation. For the horizontal plane where the upper and lower bench surfaces are located, the position where the two point clouds start to separate is the starting position of the blast area. Record the point cloud data at this time and intercept. When calculating the maximum throwing distance of the blast pile, the points obtained are within this range. By counting the extreme point coordinates of the point clouds before and after blasting, the throwing distance of the blast pile is calculated. The step height is H, the blast pile height is h, the blast step width is B, the mine rock blast pile extension distance is b, and the blast pile throwing distance is b0, wherein the throwing distance: b0=b-(B0-B) (18).

5. The UAV-based blast zone surface morphology inversion method of claim 1, wherein In the S04, the steps are as follows: S0401: Point cloud feature rough division: a color region growing algorithm is used for preliminary screening to determine the suspicious area; the same color and close distance are likely to be a type of target, and continuous scene point clouds are changed into different objects; The algorithm mainly includes two steps: (1) Segmentation: the color difference between the current seed point and the field point is less than the color difference threshold, which is considered as a cluster; (2) Merge: the color difference between clusters is less than the color difference threshold, and the number of points in the current cluster is less than the number of cluster points, which is merged with the nearest cluster; Distance calculation for RGB: After preliminary screening of the blast pile surface block degree distribution feature results, the PointNet++ algorithm is used for detailed division of the block degree; S0402: Block fine segmentation: for blast pile block segmentation, the PointNet++ algorithm is used for blast pile surface block identification and segmentation; Combined with the production situation of the mine, a volume threshold is set, and the volume exceeding the threshold is a large block; PointNet++ uses the concept of neighborhood to effectively model local features, extract local features at different scales, and obtain deep features through a multi-layer network structure; S0403: Block degree distribution calculation: based on the blast pile surface block degree segmentation results in the above, the left surface of the blast area has the most rock blocks, the right area of the blast area has the second most rock blocks, and the middle area of the blast area has the least rock blocks. From the above results, it can be seen that the left side of the blast area needs to be optimized. Bulk quantity statistical method: if x 大块i <x 左 num 左 +1Else if x 左 <x 大块i <x 右 num 中 +1Else x 右 <x 大块i num 右 +1(i=1, 2, 3, …n) wherein x 大块i : x-axis coordinate of the i-th macroblock, x 左 : left region boundary, x 右 : right region boundary, num: number of macroblocks in different regions.

Images