Multi-arm drill jumbo self-adaptive drilling track planning method based on visual navigation

Abstract

The invention provides a multi-arm drill jumbo self-adaptive drilling track planning method based on visual navigation, which relates to the technical field of visual identification, and comprises the following steps: acquiring image data, identifying a working surface change area, generating a drilling grid with topological optimization, determining a target drilling position set and a drilling direction set, and performing drilling grouping. Generating a space-time conflict constraint for each group of drill holes, and iteratively optimizing to determine a multi-arm collaborative operation sequence; in the process of executing the multi-arm operation sequence, a structured light projection technology and a visual phase shift method are adopted to obtain a deformation field of the surface of a rock mass, rock mass response characteristics including elastic deformation, plastic deformation and abnormal deformation information are established, and according to the rock mass response characteristics, a progressive B-spline track adjustment algorithm is adopted to adjust the drilling track in real time.

Description

Technical Field

[0001] The present invention relates to vision recognition technology, and particularly to an adaptive drilling trajectory planning method for a multi-boom rock drilling jumbo based on vision navigation. Background Art

[0002] As an important equipment for tunnel and mine construction, the automation and intelligence level of the multi-boom rock drilling jumbo needs to be improved urgently. The following problems and deficiencies exist in the prior art: The traditional drilling layout method mainly relies on manual experience and lacks the ability to accurately perceive and analyze the change characteristics of the working face. The fixed drilling grid mode is difficult to adapt to complex geological conditions, which easily results in low utilization rate of blasting energy and affects the construction effect and efficiency.

[0003] The operation planning of the existing multi-boom rock drilling jumbo mostly adopts a simple nearest-neighbor principle or a fixed mode, without fully considering the spatio-temporal conflicts between the booms and the optimization of task allocation. This method easily leads to interference between the booms, too long waiting time, reduces the overall operation efficiency of the equipment, lacks real-time monitoring and adaptive control capabilities, and cannot timely perceive and respond to the rock mass deformation state during the drilling process. The preset fixed drilling trajectory is difficult to cope with the uncertainties during the construction process, which easily causes potential safety hazards and affects the drilling quality.

[0004] Therefore, there is an urgent need to develop an adaptive drilling trajectory planning method for a multi-boom rock drilling jumbo based on vision navigation to improve the intelligence level and operation efficiency of the equipment. Summary of the Invention

[0005] The embodiments of the present invention provide an adaptive drilling trajectory planning method for a multi-boom rock drilling jumbo based on vision navigation, which can solve the problems in the prior art.

[0006] In the first aspect of the embodiments of the present invention, an adaptive drilling trajectory planning method for a multi-arm rock drilling jumbo based on visual navigation is provided, including: acquiring multi-view image data of a working face, extracting features from the multi-view image data, and determining the matching relationship of feature points by combining temporal consistency constraints, identifying the changed area of the working face according to the matching relationship and the point cloud differential result, generating a drilling grid with topological optimization based on the changed area of the working face, and determining a set of target drilling positions and a set of drilling directions; calculating the drilling arm task allocation cost according to the set of target drilling positions and the set of drilling directions, where the task allocation cost includes the drilling arm movement distance and the drilling direction adjustment amount; grouping the drillings based on the task allocation cost, generating spatio-temporal conflict constraints for each group of drillings, and iteratively optimizing to determine the multi-arm collaborative operation sequence; during the execution of the multi-arm operation sequence, using the structured light projection technology to monitor the micro-deformation state of the rock mass during drilling, using the visual phase shift method to obtain the deformation field of the rock mass surface, establishing a rock mass response feature including elastic deformation, plastic deformation and abnormal deformation information, and according to the rock mass response feature, using the progressive B-spline trajectory adjustment algorithm to adjust the drilling trajectory in real time.

[0007] Optionally, generating a drilling grid with topological optimization based on the changed area of the working face and determining a set of target drilling positions and a set of drilling directions includes: constructing a grid density function based on the changed area of the working face, where the grid density function comprehensively considers local curvature, the changed area distance function and the structural importance measure, applying a topological optimization algorithm to solve the optimal drilling layout, and the topological optimization algorithm adopts a hybrid strategy combining an improved SIMP method and a genetic algorithm; performing local evaluation and global evaluation on the optimal drilling layout, where the local evaluation includes energy utilization rate, and the global evaluation includes grid uniformity, and iteratively optimizing the problem area based on the evaluation results to obtain a set of target drilling positions and a set of drilling directions.

[0008] Optionally, the hybrid strategy that the topological optimization algorithm adopts to combine the improved SIMP method and the genetic algorithm includes: the improved SIMP method introduces an adaptive penalty factor and local sensitivity correction, the adaptive penalty factor is dynamically adjusted with the number of iterations, and the local sensitivity correction refines the grid of the design domain based on the stress distribution gradient; converting the result optimized by the improved SIMP method into the initial population of the genetic algorithm, the genetic algorithm adopts a real number coding method, each chromosome contains both drilling position information and drilling direction information, a crossover operator and an adaptive mutation operator are designed, the crossover operator contains a local perturbation term, the mutation step size of the adaptive mutation operator decays with the number of generations, constructing a fitness objective function including energy utilization rate and grid uniformity, and performing population iterative optimization based on the fitness objective function.

[0009] Optionally, calculate the task assignment cost for the set of target drilling positions and the set of drilling directions. The task assignment cost includes the drill arm movement distance and the drilling direction adjustment amount. Based on the task assignment cost, perform drilling grouping, generate spatio-temporal conflict constraints for each group of drillings, and iteratively optimize to determine the multi-arm collaborative operation sequence, which includes: constructing a task assignment cost function for a single drilling. The task assignment cost function includes a drill arm movement distance cost and a drilling direction adjustment cost. The drill arm movement distance cost is calculated based on the Euclidean distance between the current position of the drill arm and the target drilling position. The drilling direction adjustment cost is calculated based on the Frobenius norm of the current attitude matrix of the drill arm and the target drilling direction matrix, and introduce a sequence continuity term based on the feature distance. Based on the task assignment cost function, construct a similarity matrix between drillings, use the spectral clustering algorithm to divide the drillings into multiple initial groups, generate spatio-temporal conflict constraints according to the drill arm positions, and the spatio-temporal conflict constraints are determined based on the safety distance threshold between drill arms. Construct a global optimization objective function including a task assignment cost term, a spatio-temporal conflict constraint term, and a sequence continuity term, use the simulated annealing algorithm to solve the global optimization objective function, generate neighborhood solutions by exchanging the order of drillings within a group and exchanging drillings across groups, calculate the acceptance probability based on the change of the global optimization objective function and the temperature parameter, and iteratively optimize to obtain the multi-arm collaborative operation sequence.

[0010] Optionally, using the simulated annealing algorithm to solve the global optimization objective function includes: designing an adaptive neighborhood generation strategy, dynamically adjusting the probabilities of within-group exchange and across-group exchange according to the current iteration stage. The within-group exchange probability is inversely proportional to the temperature, and the across-group exchange probability is proportional to the temperature. When the improvement rate of the global optimization objective function in N consecutive iterations is lower than the preset threshold, expand the neighborhood search range and introduce a trajectory memory mechanism. Construct a method for calculating the acceptance probability based on historical information, record the improvement degree of the objective function of each accepted solution, adaptively adjust the acceptance threshold according to the historical improvement rate, and introduce a second-order perturbation term based on the Taylor expansion. The amplitude of the second-order perturbation term is proportional to the curvature of the current solution space. Adopt a segmented cooling strategy, use different cooling coefficients in different temperature intervals. The cooling coefficient in the initial temperature interval is greater than that in the subsequent temperature intervals. When no better solution can be obtained after several consecutive iterations, adaptively determine the re-heating temperature parameter based on the distribution of the historical optimal solutions, and maintain the characteristics of the historical optimal solutions during the re-heating process.

[0011] Optionally, the structured light projection technology is used to monitor the micro-deformation state of the rock mass during the drilling process, the visual phase-shift method is used to obtain the deformation field on the surface of the rock mass, and the rock mass response characteristics including elastic deformation, plastic deformation and abnormal deformation information are established. According to the rock mass response characteristics, the progressive B-spline trajectory adjustment algorithm is used to adjust the drilling trajectory in real time, which includes: using the structured light projection technology to project a sinusoidal grating pattern onto the surface of the rock mass, using the visual phase-shift method to obtain four images with different phases and calculate the phase difference, performing an unwrapping calculation on the phase difference to obtain the absolute phase, and calculating the three-dimensional deformation field data of the rock mass surface based on the absolute phase and the system calibration parameters; extracting the macroscopic deformation field description characteristics, local deformation abnormality indexes and temporal deformation trend characteristics of the three-dimensional deformation field data, constructing a multi-scale deformation feature space, using the RANSAC algorithm to set up an evaluation unit to iteratively identify the abnormal deformation regions in the multi-scale deformation feature space, calculating the elastic deformation parameters, plastic deformation parameters and abnormal deformation parameters according to the abnormal deformation regions, and obtaining the rock mass response characteristics by using the fuzzy inference method; constructing a multi-constraint objective function including a deformation field gradient constraint function, a drilling process parameter constraint function and an equipment kinematics constraint function according to the rock mass response characteristics, using the sliding time window strategy to locally optimize the multi-constraint objective function, and adjusting the positions of the B-spline trajectory control points based on the local optimization results to optimize and adjust the drilling trajectory in real time.

[0012] Optionally, constructing a multi-constraint objective function according to the rock mass response characteristics, performing local optimization using the sliding time window strategy, and adjusting the positions of the B-spline trajectory control points based on the local optimization results to optimize and adjust the drilling trajectory in real time includes: the deformation field gradient constraint function is calculated based on the deformation gradient tensor, the drilling process parameter constraint function integrates energy consumption terms, wear terms and time terms, and the equipment kinematics constraint function includes curvature constraints and velocity continuity constraints; establishing a multi-level optimization trigger mechanism based on time trigger, event trigger and performance trigger according to the rock mass response characteristics, and the optimization priority of the multi-level optimization trigger mechanism is dynamically associated with the rock mass risk level; when optimization is triggered, a sliding time window is defined, and the length of the sliding time window changes adaptively according to the rock mass response characteristics. The trajectory is represented as a B-spline curve within the sliding time window, and the control points of the B-spline curve include fixed control points, transition control points and optimizable control points; designing a degradation optimization strategy that perceives the calculation time, using the sequential quadratic programming method to locally optimize the multi-constraint objective function, and constructing physical constraints through the geometric relationship of the control points by using the B-spline characteristics; adjusting the positions of the B-spline trajectory control points based on the local optimization results, and using a mixed coefficient calculation method based on the influence domain of the control points for the transition control points to achieve smooth trajectory transition; automatically adjusting the optimization degree according to the time margin calculated by the sequential quadratic programming method, and adopting an adaptive mixing strategy of partial optimization results and the old trajectory to complete the trajectory adjustment within a preset time.

[0013] Optionally, obtaining multi-view image data of the working face, extracting features from the multi-view image data, and determining the matching relationship of feature points in combination with temporal consistency constraints, and identifying the changed area of the working face according to the matching relationship and the point cloud difference result includes: obtaining multi-view image data of the working face, performing camera calibration based on Zhang Zhengyou calibration method and performing image enhancement preprocessing in combination with Retinex theory, extracting multi-scale feature points from the preprocessed image using the SIFT algorithm, calculating the local structure tensor to evaluate the quality of feature points, constructing a spatial distribution uniformity constraint according to the distance between feature points, and optimizing and screening feature points; calculating the feature appearance similarity and geometric similarity for the optimized and screened feature points, dynamically adjusting the similarity weight coefficient according to the stability of feature points, verifying the geometric consistency of feature points using the RANSAC algorithm, determining the matching relationship of feature points in combination with temporal motion continuity and multi-view epipolar geometry constraints, and calculating the confidence of feature point matching based on the matching relationship; performing three-dimensional reconstruction on the multi-view image data to obtain point cloud data, performing projection transformation registration on the point cloud data, segmenting the changed area using an adaptive threshold, calculating the feature matching probability according to the matching relationship, calculating the point cloud difference probability in combination with the point cloud difference, determining the fusion weight coefficient based on the feature point matching confidence and point cloud density, and fusing the feature matching probability and the point cloud difference probability according to the fusion weight coefficient to identify the changed area of the working face.

[0014] In a second aspect of the embodiments of the present invention, an electronic device is provided, including: a processor; a memory for storing instructions executable by the processor; wherein, the processor is configured to call the instructions stored in the memory to execute the method described above.

[0015] In a third aspect of the embodiments of the present invention, a computer-readable storage medium is provided, on which computer program instructions are stored, and when the computer program instructions are executed by a processor, the method described above is implemented.

[0016] The present invention realizes the accurate identification of the changed area of the working face through multi-view image acquisition and feature matching technologies, and generates an optimal drilling grid layout based on the topology optimization method, which can effectively improve the rationality of drilling layout and the energy utilization efficiency, and provide a scientific basis for subsequent construction.

[0017] The present invention uses an improved task assignment strategy and simulated annealing algorithm to optimize the drilling grouping and operation sequence, and improves the algorithm convergence through adaptive neighborhood generation and segmented cooling mechanism, effectively solving the spatio-temporal conflict problem in multi-arm collaborative operation and significantly improving the construction efficiency.

[0018] The present invention obtains the deformation information of the rock mass in real time based on structured light monitoring and visual phase shift method, dynamically optimizes the drilling trajectory through a progressive B-spline trajectory adjustment algorithm, and adopts a multi-level optimization trigger mechanism to ensure the real-time calculation, thereby improving the construction safety and operation accuracy of the rock drilling jumbo. Brief Description of the Drawings

[0019] Figure 1 It is a schematic flow chart of the adaptive drilling trajectory planning method for a multi-arm rock drilling jumbo based on visual navigation in an embodiment of the present invention; Figure 2 It is a comparison chart of the performance of the optimization strategy; Figure 3 It is an analysis chart of the performance parameters of the multi-constraint objective function of the present invention. <able> Detailed Embodiment

[0020] To make the objectives, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0021] The technical solutions of the present invention will be described in detail below with specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be repeated in some embodiments.

[0022] Figure 1 It is a schematic flow chart of the adaptive drilling trajectory planning method for a multi-arm rock drilling jumbo based on visual navigation in an embodiment of the present invention, as Figure 1As shown, the method includes: acquiring multi-view image data of the working face, extracting features from the multi-view image data, and determining the matching relationship of feature points in combination with temporal consistency constraints, identifying the changed area of the working face according to the matching relationship and the point cloud difference result, generating a drilling grid with topological optimization based on the changed area of the working face, and determining the set of target drilling positions and the set of drilling directions; calculating the task assignment cost of the drilling arms according to the set of target drilling positions and the set of drilling directions, the task assignment cost including the moving distance of the drilling arms and the amount of drilling direction adjustment; grouping the drillings based on the task assignment cost, generating spatio-temporal conflict constraints for each group of drillings, and iteratively optimizing to determine the multi-arm collaborative operation sequence; during the execution of the multi-arm operation sequence, using the structured light projection technology to monitor the micro-deformation state of the rock mass during the drilling process, using the visual phase shift method to obtain the deformation field of the rock mass surface, establishing the rock mass response characteristics including elastic deformation, plastic deformation and abnormal deformation information, and using the progressive B-spline trajectory adjustment algorithm to adjust the drilling trajectory in real time according to the rock mass response characteristics.

[0023] Optionally, generating a drilling grid with topological optimization based on the changed area of the working face and determining the set of target drilling positions and the set of drilling directions includes: constructing a grid density function based on the changed area of the working face, the grid density function comprehensively considering local curvature, changed area distance function and structural importance measure, applying a topological optimization algorithm to solve the optimal drilling layout, the topological optimization algorithm adopting a hybrid strategy combining an improved SIMP method and a genetic algorithm; performing local evaluation and global evaluation on the optimal drilling layout, the local evaluation including energy utilization rate, the global evaluation including grid uniformity, and iteratively optimizing the problem area based on the evaluation results to obtain the set of target drilling positions and the set of drilling directions.

[0024] Exemplarily, a grid density function is constructed for the working face change area, and the three-dimensional point cloud data of the working face is collected. The local curvature value of each point is calculated through principal curvature analysis. A higher density weight (weight value is 1.5) is assigned to the area with larger curvature (for example, the area with curvature value greater than 0.05), and a lower density weight (weight value is 0.8) is assigned to the flat area with smaller curvature (for example, the area with curvature value less than 0.01). At the same time, the shortest distance from the points on the working face to the change area is calculated to form a distance field. For the points within 50 cm from the boundary of the change area, weights are assigned according to the inverse distance function, and the closer the distance, the higher the weight. Specifically, when the distance is d, the weight is 1 / (d + 0.1). In addition, combined with on-site survey and rock mass structure analysis, the importance of the structure is measured. For example, the importance of fault and fracture-intensive areas is set to a high value (value is 1.8), and the importance of homogeneous areas is set to a medium value (value is 1.0). Finally, the above three indicators are weighted and summed to form a comprehensive density function, with weights of local curvature 0.3, distance function 0.4, and structural importance 0.3 respectively.

[0025] A hybrid strategy combining the improved SIMP method and the genetic algorithm is used to solve the optimal drilling layout. Local evaluation and global evaluation are carried out on the optimal drilling layout. Local evaluation mainly examines the energy utilization rate. By simulating the energy distribution during the drilling blasting process, the energy coverage rate and energy utilization efficiency are calculated. Specifically, an energy influence sphere with a radius of 2 meters is established for each drill hole, and the overlap rate between the spheres is calculated, and the overlap rate is required to be controlled between 15% - 25%. The area with an energy utilization rate lower than 85% is marked as a problem area and needs to be adjusted. Global evaluation includes grid uniformity analysis. By calculating the standard deviation and maximum deviation of the drill hole spacing, the standard deviation is required not to exceed 15% of the designed spacing, and the maximum deviation does not exceed 30% of the designed spacing. The area with non-compliant uniformity is also marked as a problem area. For the problem area, the system automatically adjusts the drill hole position and direction. For the area with low energy utilization rate, 1 - 3 auxiliary drill holes are added; for the area with poor grid uniformity, the uniformity is optimized by locally fine-tuning the drill hole position (the adjustment range does not exceed 20% of the designed spacing). The adjusted plan is evaluated again until the requirements are met or the maximum number of iterations of 5 times is reached. The finally determined set of drill hole positions and set of directions are used as the input for the subsequent drill hole task planning.

[0026] In this embodiment, by constructing a comprehensive grid density function, the accurate identification and quantification of the working face characteristics are realized. By using a hybrid strategy combining the SIMP method and the genetic algorithm, the optimization efficiency and quality of the drilling layout are significantly improved. Through the evaluation mechanism combining local evaluation and global evaluation, the finally generated drill hole grid not only meets the requirements of energy utilization rate but also ensures the grid uniformity, improves the blasting efficiency, and reduces energy waste.

[0027] Optionally, the topology optimization algorithm adopts a hybrid strategy combining the improved SIMP method and the genetic algorithm, which includes: the improved SIMP method introduces an adaptive penalty factor and local sensitivity correction. The adaptive penalty factor is dynamically adjusted according to the number of iterations. The local sensitivity correction refines the mesh of the design domain based on the stress distribution gradient. The result optimized by the improved SIMP method is converted into the initial population of the genetic algorithm. The genetic algorithm uses real-number encoding. Each chromosome contains both the drilling position information and the drilling direction information. A crossover operator and an adaptive mutation operator are designed. The crossover operator contains a local perturbation term. The mutation step size of the adaptive mutation operator decays with the number of generations. A fitness objective function including energy utilization rate and grid uniformity is constructed, and population iterative optimization is performed based on the fitness objective function.

[0028] Exemplarily, in the improved SIMP method, two technical points of an adaptive penalty factor and local sensitivity correction are introduced. The adaptive penalty factor is used to control the penalty degree of the intermediate density to ensure that the final result tends to a 0-1 distribution. The initial penalty factor is set to 3.0 and is dynamically adjusted according to the current iteration round during the optimization process. The specific adjustment strategy is as follows: in the first 30 iterations, the penalty factor increases by 0.2 every 10 iterations; in the 31-60 iterations, the penalty factor increases by 0.3 every 5 iterations; after the 61st iteration, the penalty factor increases by 0.4 every 3 iterations until the maximum value of 5.0 is reached. This adjustment strategy ensures that more intermediate density values are allowed to exist in the initial stage of optimization, and gradually forces it to tend to a 0-1 distribution in the later stage. For local sensitivity correction, first calculate the stress distribution of the entire design domain. Taking the drilling grid optimization as an example, the stress distribution is calculated by simulating the energy propagation of the drill hole in the rock mass. In each iteration, calculate the stress gradient and identify the regions where the stress gradient is greater than the threshold of 0.25. Refine the mesh of these regions, doubling the original mesh resolution. For example, the original regional mesh size is 10×10 cm, and after refinement, it becomes 5×5 cm. To ensure the calculation efficiency, the maximum proportion of the refined area is set not to exceed 30% of the overall design domain. Recalculate the local sensitivity on the refined mesh and map it back to the original mesh using linear interpolation to obtain the corrected sensitivity distribution.

[0029] After the SIMP method completes the optimization, the result is converted into the initial population of the genetic algorithm. During the conversion process, first identify the voxel positions with density values greater than 0.85 in the SIMP result as candidate drilling centers. Then, use the density gradient information to determine the drilling direction, and the direction vector is consistent with the density gradient direction. For each candidate drill hole, a 6D chromosome containing the position coordinates (x, y, z) and the direction vector (dx, dy, dz) is established. To increase the population diversity, 70% of the initial individuals come from the SIMP optimization result, and 30% are randomly generated. The population size is set to 200, which is large enough to cover the solution space without causing an excessive computational burden.

[0030] The genetic algorithm uses a real - number encoding method, and each chromosome contains the position and direction information of multiple drill holes. To improve the optimization efficiency, special crossover operators and adaptive mutation operators are designed. The crossover operator adopts the Simulated Binary Crossover (SBX) method, and the crossover probability is set to 0.85. A local perturbation term is introduced during the crossover process, and the perturbation range is 5% of the variable range to increase the local search ability. The initial mutation probability of the adaptive mutation operator is set to 0.1, and the initial mutation step size is 10% of the variable range. As the number of generations increases, the mutation step size decreases according to the exponential decay law. Specifically, it decays to 0.9 times the original value every 50 generations. For example, if the initial mutation step size is 10%, it becomes 9% after the 50th generation and 8.1% after the 100th generation, and so on. This strategy ensures that the algorithm has sufficient global search ability in the initial stage and gradually focuses on local fine - tuning search in the later stage.

[0031] The fitness objective function comprehensively considers two aspects: energy utilization rate and grid uniformity. The calculation of the energy utilization rate is based on the energy coverage model between drill holes. An energy - influence sphere is established for each drill hole, and the energy coverage ratio of the working face is calculated. The grid uniformity is evaluated by calculating the coefficient of variation of the adjacent drill - hole spacing. A coefficient of variation not exceeding 0.2 is considered excellent. The weights of the two indicators are 0.6 and 0.4 respectively. After the algorithm converges, the individual with the highest fitness is selected as the final drill - hole layout plan.

[0032] Figure 2For the optimized strategy performance comparison chart, the chart compares the performance of the hybrid optimization strategy of this application and the prior art in five key performance indicators. Specifically: in terms of the convergence speed, the method of this application reaches 0.95, which is significantly better than 0.72 of the traditional SIMP method and 0.68 of the simple GA, reflecting the obvious advantage of the hybrid strategy in optimization efficiency. In terms of the grid uniformity index, the method of this application reaches 0.92, which is about 18% higher than 0.78 of the traditional SIMP method and 0.74 of the simple GA, indicating that the hybrid strategy can effectively improve the grid quality. In terms of the energy utilization rate, the method of this application is 0.85, which is about 12% higher than 0.76 of the traditional SIMP method and 0.72 of the simple GA, showing that the practicality of the optimization result is better. In terms of the computational efficiency index, the method of this application reaches 0.82, which is higher than 0.75 of the traditional SIMP method and 0.65 of the simple GA, proving that the hybrid strategy has good computational performance while ensuring the optimization effect. In terms of the diversity of solutions, the method of this application is 0.86, which is significantly better than 0.65 of the traditional SIMP method and slightly higher than 0.78 of the simple GA, indicating that the hybrid strategy not only maintains the diversity characteristics of GA but also improves the quality of solutions. Generally speaking, the hybrid optimization strategy proposed in this application, which combines the improved SIMP method and the genetic algorithm, effectively solves the intermediate density problem and the local accuracy deficiency problem in the traditional SIMP method through the adaptive penalty factor and the local sensitivity correction technology. The method of initializing the genetic algorithm population based on the SIMP result, combined with the specially designed crossover operator and the adaptive mutation operator, significantly improves the convergence speed and the quality of solutions of the algorithm, achieving the organic combination of the global optimization and the local fine adjustment of the drilling layout.

[0033] Optionally, calculate the task assignment cost for the set of target drilling positions and the set of drilling directions, where the task assignment cost includes the drill arm movement distance and the drilling direction adjustment amount; group the drillings based on the task assignment cost, generate spatio-temporal conflict constraints for each group of drillings, and iteratively optimize to determine the multi-arm collaborative operation sequence, including: constructing a task assignment cost function for a single drilling, where the task assignment cost function includes a drill arm movement distance cost and a drilling direction adjustment cost, the drill arm movement distance cost is calculated based on the Euclidean distance between the current position of the drill arm and the target drilling position, the drilling direction adjustment cost is calculated based on the Frobenius norm of the current attitude matrix of the drill arm and the target drilling direction matrix, and a sequence continuity term based on the feature distance is introduced; constructing a similarity matrix between drillings based on the task assignment cost function, using the spectral clustering algorithm to divide the drillings into multiple initial groups, generating spatio-temporal conflict constraints based on the drill arm positions, where the spatio-temporal conflict constraints are determined based on the safety distance threshold between drill arms; constructing a global optimization objective function including a task assignment cost term, a spatio-temporal conflict constraint term, and a sequence continuity term, using the simulated annealing algorithm to solve the global optimization objective function, generating neighborhood solutions by exchanging the order of drillings within a group and exchanging drillings across groups, calculating the acceptance probability based on the change in the global optimization objective function and the temperature parameter, and iteratively optimizing to obtain the multi-arm collaborative operation sequence.

[0034] Exemplarily, when constructing the task assignment cost function for a single drilling, both the drill arm movement distance cost and the drilling direction adjustment cost are comprehensively considered. For the drill arm movement distance cost, the Euclidean distance between the current position of the drill arm and the target drilling position is used for calculation. Taking a certain actual application scenario as an example, the current position coordinates of the drill arm are (8.5m, 4.2m, 1.8m), and the target drilling position coordinates are (12.3m, 5.7m, 2.1m), then the Euclidean distance is 4.25 meters. After the distance is calculated, it is converted into a cost value through a piecewise function. For example, when the distance is less than 5 meters, the cost is the distance value multiplied by 1.0; when the distance is between 5 - 10 meters, the cost is the distance value multiplied by 1.2; when the distance is greater than 10 meters, the cost is the distance value multiplied by 1.5, reflecting the principle that the farther the distance, the higher the cost.

[0035] For the cost of adjusting the drilling direction, it is calculated based on the Frobenius norm of the current posture matrix and the target drilling direction matrix of the drill boom. The posture matrix is represented by a 3×3 rotation matrix, and the direction adjustment amount is obtained by calculating the Frobenius norm of the difference between the current posture matrix and the target posture matrix. For example, when the difference between the current posture and the target posture is small (the norm value is 0.3), the direction adjustment cost is set to 0.9; when the difference is medium (the norm value is 0.6), the direction adjustment cost is set to 1.8; when the difference is large (the norm value is 1.0), the direction adjustment cost is set to 3.0. In addition, a sequence continuity term based on the feature distance is introduced. When the feature distances (including position and direction features) between adjacent operation drill holes are small, an additional continuity reward is allocated, and the reward value is between 0.5 and 1.5 to encourage the system to select a drill hole sequence with good continuity and reduce unnecessary adjustment times of the equipment.

[0036] The moving distance cost of the drill boom for a single drill hole, the drilling direction adjustment cost, and the sequence continuity term are weighted and combined with weights of 0.5, 0.35, and 0.15 respectively to obtain the comprehensive task assignment cost. Based on this cost function, a similarity matrix between drill holes is further constructed. The matrix element value represents the task similarity between two drill holes, and the similarity value is between 0 and 1. The larger the value, the higher the similarity. For example, for two drill holes with similar positions and directions, the similarity value can reach above 0.85; for drill holes with similar positions but large direction differences, the similarity value is about 0.6; for drill holes with large differences in both position and direction, the similarity value is as low as below 0.3.

[0037] Based on the similarity matrix between drill holes, the spectral clustering algorithm is used to divide the drill holes into multiple initial groups. First, the Laplace transform is performed on the similarity matrix, then its eigenvectors are calculated, and the K most informative eigenvectors are selected to construct a new feature space. The K-means clustering method is applied in the new space to obtain the drill hole grouping result. The value of K is determined according to the number of drill booms and the size of the working face, and is usually set to 1.2 - 1.5 times the number of drill booms. For example, for a 3-arm drill jumbo, the value of K is set to 4, forming 4 initial drill hole groups, and each group contains about 25% of the drill holes.

[0038] After obtaining the initial grouping, spatio-temporal conflict constraints are generated according to the drill boom positions. Based on the operating range and safety distance requirements of the drill boom, the areas where conflicts may occur between drill booms are calculated. For example, if the extension limit of the drill boom is 5 meters and the safety distance threshold is set to 2 meters, then when the distance between the operating points of two drill booms is less than 7 meters, there may be a conflict risk. For such cases, the relevant drill hole pairs are marked and spatio-temporal conflict constraints are generated, and the conflict drill hole pairs that cannot operate simultaneously are recorded in the constraints. In practical applications, for a certain mine working face with a width of 15 meters, the generated spatio-temporal conflict constraints contain a total of 42 conflict drill hole pairs, accounting for about 18% of the total number of drill hole pairs.

[0039] Finally, a global optimization objective function including a task assignment cost term, a spatio-temporal conflict constraint term, and a sequence continuity term is constructed. The task assignment cost term is the sum of the total task costs of all drill arms. The spatio-temporal conflict constraint term is represented in the form of a penalty function, and the penalty value for violating the constraint once is 100. The sequence continuity term reflects the smoothness of the operation sequence. The weights of the three terms are 0.65, 0.25, and 0.1 respectively. The simulated annealing algorithm is applied to solve the global optimization objective function. In each iteration, neighborhood solutions are generated by swapping the within-group drilling order and cross-group drilling, and the acceptance probability is calculated based on the change in the objective function value and the current temperature to obtain the final multi-arm collaborative operation sequence.

[0040] In this embodiment, by comprehensively considering the task assignment cost function of the drill arm movement distance cost, the drilling direction adjustment cost, and the sequence continuity, a precise quantitative evaluation of the drilling task is achieved. The spectral clustering algorithm is used for drilling grouping and spatio-temporal conflict constraints are generated based on the drill arm positions, effectively avoiding the risk of equipment interference in multi-arm operations. The global optimization method combined with the simulated annealing algorithm significantly improves the multi-arm collaborative efficiency and achieves a double improvement in equipment utilization rate and operation safety.

[0041] Optionally, using the simulated annealing algorithm to solve the global optimization objective function includes: designing an adaptive neighborhood generation strategy, dynamically adjusting the probabilities of within-group swapping and cross-group swapping according to the current iteration stage. The within-group swapping probability is inversely proportional to the temperature, and the cross-group swapping probability is directly proportional to the temperature; when the improvement rate of the global optimization objective function in N consecutive iterations is lower than the preset threshold, expand the neighborhood search range and introduce a trajectory memory mechanism; construct an acceptance probability calculation method based on historical information, record the improvement degree of the objective function of each accepted solution, adaptively adjust the acceptance threshold according to the historical improvement rate, and introduce a second-order perturbation term based on Taylor expansion, and the amplitude of the second-order perturbation term is proportional to the curvature of the current solution space; adopt a segmented cooling strategy, using different cooling coefficients in different temperature intervals, and the cooling coefficient in the initial temperature interval is greater than that in the subsequent temperature intervals; when no better solution can be obtained after several consecutive iterations, adaptively determine the re-heating temperature parameter based on the historical optimal solution distribution and maintain the characteristics of the historical optimal solution during the re-heating process.

[0042] Exemplarily, in terms of the adaptive neighborhood generation strategy, the probabilities of in-group exchange and cross-group exchange are dynamically adjusted according to the current iteration stage. In-group exchange refers to adjusting the operation sequence of drilling within the drilling group responsible for the same drilling arm, and cross-group exchange refers to allocating a certain drilling to another drilling arm group. The system sets the in-group exchange probability to be inversely proportional to the temperature, and the cross-group exchange probability to be proportional to the temperature. Specifically, when the temperature is above 80% of the initial temperature, the in-group exchange probability is 0.3 and the cross-group exchange probability is 0.7; when the temperature drops to 30%-80% of the initial temperature, the in-group exchange probability is 0.5 and the cross-group exchange probability is 0.5; when the temperature is below 30% of the initial temperature, the in-group exchange probability is 0.8 and the cross-group exchange probability is 0.2. This design makes the algorithm tend to explore widely in the high-temperature stage and find a better drilling arm task allocation scheme through cross-group exchange; while in the low-temperature stage, it focuses on local optimization and finely adjusts the operation sequence through in-group exchange.

[0043] When the improvement rate of the global optimization objective function in N consecutive iterations is lower than the preset threshold, the neighborhood search range is automatically expanded and the trajectory memory mechanism is introduced. The value of N is usually set to 10% of the current iteration number, with a minimum of 20 times and a maximum of 100 times; the preset threshold is 0.5%. For example, when the improvement rate of the objective function in 30 consecutive iterations is less than 0.5%, the system expands the neighborhood search range to 1.5 times the original, and at the same time activates the trajectory memory mechanism. The trajectory memory mechanism maintains a memory bank containing historical search paths, recording the solutions and their evaluation values visited in the last 200 iterations. After expanding the search range, the system preferentially explores the areas in the historical memory bank that have not been fully explored, avoiding repeated searches and improving the algorithm efficiency.

[0044] In terms of acceptance probability calculation, a method based on historical information is constructed. In the traditional simulated annealing algorithm, the acceptance probability is only related to the change of the objective function and the temperature in the current iteration. In the improved method, the improvement degree of the objective function of each accepted solution is additionally recorded, and a queue of improvement rates with a length of 30 is maintained. According to the statistical characteristics (mean and standard deviation) of the historical improvement rates, the acceptance threshold is adaptively adjusted. For example, when the mean of the historical improvement rate is 2% and the standard deviation is 0.5%, the acceptance threshold is adjusted to 1.5% accordingly; when the mean is 0.8% and the standard deviation is 0.2%, the acceptance threshold is adjusted to 0.65%. This mechanism enables the algorithm to dynamically adjust the acceptance strategy according to the actual search situation at different stages, avoiding premature convergence or excessive exploration. In addition, a second-order perturbation term based on Taylor expansion is introduced into the acceptance probability calculation, and the amplitude of this perturbation term is proportional to the curvature of the current solution space. Specifically, the curvature of the objective function near the current solution is estimated through multiple neighborhood samplings. The larger the curvature value, the larger the amplitude of the perturbation term. For example, when the estimated curvature value is 0.05, the amplitude of the perturbation term is set to 0.02; when the curvature value is 0.15, the amplitude increases to 0.06. At low temperatures, this second-order perturbation can help the algorithm jump out of local optima, especially in regions where the solution space is highly non-linear.

[0045] In terms of the cooling strategy, a segmented cooling method is adopted, and different cooling coefficients are used in different temperature intervals. The entire temperature range is divided into three segments: the high-temperature segment (70% - 100% of the initial temperature), the medium-temperature segment (30% - 70% of the initial temperature), and the low-temperature segment (below 30% of the initial temperature). The cooling coefficient in the high-temperature segment is 0.95, in the medium-temperature segment is 0.98, and in the low-temperature segment is 0.99. This segmented strategy enables the algorithm to cool down quickly in the initial stage, narrowing the search range; while cooling down slowly in the later stage, finely exploring promising regions.

[0046] When no better solution can be obtained after 50 consecutive iterations, the re-heating temperature parameter is adaptively determined based on the distribution of historical optimal solutions. First, the distribution characteristics of historical optimal solutions are analyzed, including the aggregation degree and diversity of the solutions. If the historical optimal solutions are highly aggregated in a few regions, it indicates that the current situation may be trapped in local optima, and at this time, the re-heating temperature is set relatively high (about 60% of the initial temperature); if the distribution of historical optimal solutions is relatively dispersed, it indicates that the search space is relatively flat, and at this time, the re-heating temperature is set relatively low (about 30% of the initial temperature). During the re-heating process, the key features of the historical optimal solutions are maintained, such as the core allocation mode of the drilling group and the efficient operation sequence segments, and only the non-critical parts are perturbed, so that it can jump out of local optima without completely losing the existing optimization results.

[0047] The adaptive neighborhood generation strategy proposed in this embodiment realizes a smooth transition of the search process from global exploration to local fine-tuning, avoiding the algorithm from falling into local optimality. The acceptance probability calculation method based on historical information and the introduction of the second-order perturbation term enhance the adaptability of the algorithm and its ability to jump out of local optimality. The combined action of the segmented cooling strategy and the reheating mechanism realizes a wider exploration of the solution space while maintaining the key features of the historical optimal solution. Compared with the traditional simulated annealing algorithm, both the optimization effect and the convergence speed are improved.

[0048] Optionally, the structured light projection technology is used to monitor the micro-deformation state of the rock mass during the drilling process, the visual phase-shifting method is used to obtain the deformation field on the surface of the rock mass, and the rock mass response characteristics including elastic deformation, plastic deformation, and abnormal deformation information are established. According to the rock mass response characteristics, the progressive B-spline trajectory adjustment algorithm is used to adjust the drilling trajectory in real time, including: using the structured light projection technology to project a sine grating pattern onto the surface of the rock mass, using the visual phase-shifting method to obtain four images with different phases and calculate the phase difference, performing an unfolding calculation on the phase difference to obtain the absolute phase, and calculating the three-dimensional deformation field data of the rock mass surface based on the absolute phase and the system calibration parameters; extracting the macroscopic deformation field description features, local deformation abnormality indicators, and temporal deformation trend features of the three-dimensional deformation field data, constructing a multi-scale deformation feature space, using the RANSAC algorithm to set up an evaluation unit to iteratively identify the abnormal deformation regions in the multi-scale deformation feature space, calculating the elastic deformation parameters, plastic deformation parameters, and abnormal deformation parameters according to the abnormal deformation regions, and obtaining the rock mass response characteristics by using the fuzzy inference method; constructing a multi-constraint objective function including a deformation field gradient constraint function, a drilling process parameter constraint function, and an equipment kinematics constraint function according to the rock mass response characteristics, using the sliding time window strategy to locally optimize the multi-constraint objective function, and adjusting the positions of the B-spline trajectory control points based on the local optimization results to optimize and adjust the drilling trajectory in real time.

[0049] Exemplarily, the structured light projection system consists of a projection unit and an imaging unit. The projection unit includes a DLP projector (resolution 1920×1080) and an optical focusing device. The imaging unit includes two industrial cameras (resolution 2448×2048), with a camera acquisition frequency of up to 30Hz. The system is installed at the front end of the drill boom, maintaining a fixed position relative to the drill bit, with a projection distance of 1.5-3 meters. The projected pattern is a sinusoidal grating with a spatial frequency of 32 lines / meter and a projection light intensity of 2000 lumens, ensuring clarity even in harsh environments. During monitoring, the system projects a sinusoidal grating pattern onto the rock surface. The projected sinusoidal grating consists of four phase versions with a phase difference of 90 degrees (i.e., 0°, 90°, 180°, and 270°). By controlling the projector to rapidly switch between grating patterns of different phases and simultaneously triggering the cameras to capture the corresponding images, four sets of phase images are acquired. The two cameras acquire data synchronously, with an acquisition time difference of less than 2 milliseconds, ensuring consistent phase measurements during rock micro-deformation. After acquiring four phase images, the visual phase shift method is used to calculate the phase difference. Specifically, for each pixel position, the wrapped phase is calculated using a four-step phase shift method based on the grayscale values of the four phase images. Because the wrapped phase range is limited to -π to π, phase unwrapping is required to obtain the absolute phase. Phase unwrapping uses an improved quality-guided algorithm with a quality threshold of 0.75. The unwrapping path begins in high-quality areas and gradually expands to low-quality areas, effectively avoiding unwrapping errors in discontinuous areas.

[0050] After phase unwrapping is complete, the system's calibration parameters are combined to calculate the 3D deformation field data of the rock mass surface. System calibration is based on Zhang's calibration method, using a 9×12 grid calibration plate, with a calibration error of less than 0.05 pixels. A mapping relationship between phase and 3D coordinates is established based on the calibration parameters, and the precise 3D coordinates of the rock mass surface are calculated. By comparing the changes in 3D coordinates before and after drilling, the deformation field of the rock mass surface is determined. In practical applications, the system's 3D reconstruction accuracy reaches 0.1 mm, enabling detection of minute rock mass deformations.

[0051] When extracting the features of three-dimensional deformation field data, it is done at three levels: macroscopic, local, and temporal. The characteristics describing the macroscopic deformation field include the distribution of deformation amplitude, the main deformation direction, and the deformation field gradient. By statistically analyzing the deformation field data, the average deformation (typical value is 0.2-5 mm), standard deviation (generally 30%-60% of the average value), and maximum deformation (usually 3-5 times the average value) are calculated. The main deformation direction is determined by principal component analysis and expressed as a three-dimensional unit vector. The deformation field gradient reflects the spatial rate of change of deformation. Areas with large gradient values (greater than 0.05) usually indicate the presence of local stress concentration.

[0052] The local deformation abnormality index mainly detects the discontinuities and singular points in the deformation field. Using the local similarity index calculation method, the rock mass surface is divided into grid cells of 10×10 cm, and the deformation similarity between each cell and its adjacent cells is calculated. When the similarity is lower than the threshold of 0.75, it is marked as an abnormal area. At the same time, the local curvature is calculated, and the area where the curvature value exceeds the threshold of 0.2 is also determined as an abnormal area. The time-series deformation trend characteristics are obtained through the time-series analysis of continuous monitoring data. The system collects data at a frequency of 5 Hz and constructs deformation time-series data with a duration of 10 seconds. The deformation rate, deformation acceleration, and deformation stability index are calculated through linear regression analysis.

[0053] Integrate the above features to form a multi-scale deformation feature space, and use the RANSAC (Random Sample Consensus) algorithm to identify abnormal deformation areas. The RANSAC algorithm sets the inlier threshold to 85% of the total data volume, the maximum number of iterations to 200, and each evaluation unit contains 50 randomly sampled points. Through iterative calculations, the abnormal points that do not conform to the normal deformation model are removed, and finally the abnormal deformation areas are obtained.

[0054] According to the abnormal deformation areas and overall deformation characteristics, calculate the elastic deformation parameters, plastic deformation parameters, and abnormal deformation parameters. The elastic deformation parameters include the estimated elastic modulus (e.g., about 50 - 80 GPa for hard rock) and the proportion of elastic deformation (usually accounting for 60% - 80% of the total deformation). The plastic deformation parameters include the area ratio of the plastic deformation area (generally less than 20%) and the cumulative amount of plastic deformation. The abnormal deformation parameters include the number of abnormal areas, location distribution, and abnormal degree grading (divided into three levels: mild, moderate, and severe).

[0055] Use the fuzzy inference method to combine the above parameters to obtain the rock mass response characteristics. The fuzzy inference system contains 15 core rules, with the above parameters as inputs and the rock mass stability rating (four levels: safe, attention, warning, dangerous) and the recommended adjustment plan for drilling parameters as outputs. For example, when the proportion of elastic deformation is higher than 75%, the area ratio of the plastic deformation area is less than 10%, and there is no severe abnormal deformation, the rock mass stability rating is "safe"; when the proportion of elastic deformation is between 60% - 75%, the area ratio of the plastic deformation area is between 10% - 20%, and there are 1 - 2 moderate abnormal deformations, the rock mass stability rating is "attention", and it is recommended to appropriately reduce the feed rate. The rock mass response characteristics obtained through fuzzy inference provide a scientific basis for subsequent drilling trajectory adjustment.

[0056] In this embodiment, through structured light projection technology and visual phase shift method, high-precision monitoring of the micro-deformation state of rock masses is achieved, overcoming the limitations of traditional contact measurement methods. The multi-scale deformation feature extraction method comprehensively captures the deformation features of rock masses from three levels: macroscopic, local, and temporal. The RANSAC algorithm effectively identifies abnormal deformation areas. The fuzzy inference system converts complex deformation parameters into intuitive rock mass response features and stability ratings, providing a scientific basis for real-time adjustment of borehole trajectories, significantly improving the perception ability of rock mass state changes and the warning accuracy.

[0057] Optionally, a multi-constraint objective function is constructed according to the rock mass response features, and a sliding time window strategy is adopted for local optimization. Based on the local optimization results, the positions of the B-spline trajectory control points are adjusted. The real-time optimization and adjustment of the borehole trajectory include: the deformation field gradient constraint function is calculated based on the deformation gradient tensor, the drilling process parameter constraint function integrates energy consumption terms, wear terms, and time terms, and the equipment kinematics constraint function includes curvature constraints and velocity continuity constraints; a multi-level optimization trigger mechanism based on time trigger, event trigger, and performance trigger is established according to the rock mass response features, and the optimization priority of the multi-level optimization trigger mechanism is dynamically associated with the rock mass risk level; when optimization is triggered, a sliding time window is defined, and the length of the sliding time window changes adaptively according to the rock mass response features. Within the sliding time window, the trajectory is represented as a B-spline curve, and the control points of the B-spline curve include fixed control points, transition control points, and optimizable control points; a degradation optimization strategy considering calculation time is designed, and the sequential quadratic programming method is used to locally optimize the multi-constraint objective function. Physical constraints are constructed through the geometric relationship of the control points using the B-spline characteristics; based on the local optimization results, the positions of the B-spline trajectory control points are adjusted, and for the transition control points, a hybrid coefficient calculation method based on the influence domain of the control points is adopted to achieve smooth trajectory transition; the optimization degree is automatically adjusted according to the time margin calculated by the sequential quadratic programming method, and an adaptive mixing strategy of partial optimization results and the old trajectory is adopted to complete the trajectory adjustment within a preset time.

[0058] Exemplarily, the multi-constraint objective function includes three main parts: the deformation field gradient constraint function, the drilling process parameter constraint function, and the equipment kinematics constraint function. The deformation field gradient constraint function is calculated based on the deformation gradient tensor and reflects the spatial distribution characteristics of rock mass deformation. In specific implementation, the three-dimensional space is discretized into a grid of 20×20×20 cm, and the deformation gradient tensor of each grid point is calculated. The eigenvalues of the deformation gradient tensor represent the main deformation directions and amplitudes. When the maximum eigenvalue exceeds the threshold of 0.08, it indicates that the deformation gradient in this area is too large, and the constraint function gives a high penalty value; when all eigenvalues are less than the threshold of 0.03, the constraint function gives a low penalty value. In this way, the drilling trajectory is guided to avoid areas with excessive deformation gradients, reducing the damage to the rock mass structure. The drilling process parameter constraint function integrates the energy consumption term, the wear term, and the time term. The energy consumption term is calculated based on the drilling power, and the ratio of drilling depth to energy consumption is used as an index, with a standard value of 0.8 - 1.2 kWh of electricity consumption per meter of drilling depth. The wear term evaluates the degree of bit wear by monitoring the torque fluctuation and vibration characteristics of the bit, and a wear accumulation threshold of 100 standard units is set. The time term evaluates the drilling speed, and the standard drilling speed in normal rock mass is 2 - 4 m / min. The three sub-terms are combined with weights to form the drilling process parameter constraint function, with weights of 0.4 for the energy consumption term, 0.35 for the wear term, and 0.25 for the time term. The equipment kinematics constraint function includes curvature constraint and velocity continuity constraint. The curvature constraint limits the maximum curvature of the drilling trajectory to prevent the drill pipe from getting stuck due to an overly curved trajectory. The velocity continuity constraint ensures that the drilling speed changes smoothly during the trajectory adjustment process, with a maximum allowable acceleration of 0.5 m / s² and a maximum allowable deceleration of 0.8 m / s². These constraint parameters are determined according to the actual operating characteristics of the equipment to ensure that the trajectory adjustment does not exceed the equipment's capacity range.

[0059] Based on the rock mass response characteristics, a multi-level optimization trigger mechanism is established, including three modes: time trigger, event trigger, and performance trigger. The time trigger performs trajectory optimization at a fixed time interval (usually 2 seconds); the event trigger immediately triggers optimization when a significant change (change amplitude exceeds 25%) in the rock mass response characteristics is detected; the performance trigger triggers optimization when the drilling parameters (such as drilling speed, torque) deviate from the normal range (deviation exceeds 30%). The priority of the optimization trigger mechanism is dynamically associated with the rock mass risk level. When the rock mass stability rating is "warning" or "dangerous", the event trigger mode is preferentially adopted, and the optimization frequency is increased.

[0060] After triggering the optimization, a sliding time window is defined for local optimization. The length of the sliding time window adaptively changes according to the rock mass response characteristics. The window length is longer (5 - 8 seconds) in stable rock masses and shorter (2 - 3 seconds) in unstable rock masses. Within the sliding time window, the trajectory is represented as a cubic B-spline curve, and the spacing between control points is 1 / 4 of the window length. B-spline control points are divided into three categories: fixed control points (the starting point of the trajectory, the end point of the executed part), transition control points (connecting the fixed segment and the optimizable segment), and optimizable control points (the unexecuted part). For example, within a 5-second time window, 9 control points are set, including 2 fixed control points, 2 transition control points, and 5 optimizable control points.

[0061] A degradation optimization strategy considering computing time is also designed to automatically adjust the optimization complexity according to the available computing time (usually 100 - 200 milliseconds). When the computing time is sufficient, the complete sequential quadratic programming method is used to optimize the multi-constrained objective function. When the computing time is limited, the optimization accuracy requirement is reduced (for example, the convergence threshold is relaxed from 0.001 to 0.01) or the number of optimization variables is reduced (for example, some adjacent control points are grouped together). The iteration step size of the sequential quadratic programming method is set to 0.05, the maximum number of iterations is 50, and the convergence threshold is 0.005.

[0062] Physical constraints are constructed using the B-spline characteristics through the geometric relationships of control points. The spacing between adjacent control points is not less than the minimum safety distance (usually 10 centimeters), and the angle between the tangent vector of the control point and the drill bit direction does not exceed the maximum allowable angle (usually 5 degrees).

[0063] Based on the local optimization results, the positions of the B-spline trajectory control points are adjusted. For transition control points, a hybrid coefficient calculation method based on the influence domain of control points is used to achieve smooth trajectory transition. Specifically, for each transition control point, the influence weights from fixed control points and optimizable control points are calculated, and the weights are inversely proportional to the distance. For example, a transition control point closer to a fixed control point has a smaller position adjustment amplitude (adjustment ratio is 10% - 30%); a transition control point closer to an optimizable control point has a larger position adjustment amplitude (adjustment ratio is 50% - 80%). The optimization degree is automatically adjusted according to the time margin calculated by the sequential quadratic programming method. When the computing time is sufficient (the remaining time exceeds 50% of the allocated time), the system performs more refined trajectory optimization, increasing the number of iterations or refining the control points. When the computing time is tight (the remaining time is less than 20% of the allocated time), the system terminates the optimization process in advance and uses the current optimal result.

[0064] An adaptive hybrid strategy of using partial optimization results and the old trajectory is adopted to complete the trajectory adjustment within a preset time. The mixing ratio is dynamically determined according to the improvement degree of the optimization results. When the improvement is significant (the objective function value decreases by more than 30%), the proportion of the new trajectory is high (80%-100%); when the improvement is slight (the objective function value decreases by less than 10%), the proportion of the new trajectory is low (20%-50%). Through this adaptive hybrid strategy, it is ensured that the trajectory adjustment is stable and effective, avoiding equipment shocks caused by drastic changes.

[0065] Figure 3 This is an analysis diagram of the performance parameters of the multi-constraint objective function of the present invention. The diagram shows the changing trends of three key parameters in the multi-constraint objective function during the iterative optimization process. The eigenvalue of the deformation gradient gradually decreases from the initial 0.10 to below the safety threshold of 0.03, indicating that the optimized trajectory effectively avoids the high-deformation area; the drilling energy consumption ratio is optimized from 0.11 to 0.08; the trajectory curvature decreases from 0.06 to 0.035, meeting the requirements of the equipment kinematic constraints. The convergence processes of the three parameters show that the optimization strategy of this application can effectively coordinate multiple constraints and achieve the comprehensive optimization of the drilling trajectory. The multi-constraint objective function constructed in this embodiment comprehensively considers the deformation field gradient, drilling process parameters, and equipment kinematic constraints, realizing the scientific optimization of the drilling trajectory. The multi-level optimization trigger mechanism based on the rock mass risk level ensures that the system can respond to rock mass changes in a timely manner. The sliding time window and the B-spline trajectory representation method make the optimization process efficient and controllable, and the degradation optimization strategy and the adaptive hybrid strategy considering the calculation time ensure the real-time performance and stability of the trajectory adjustment.

[0066] Optionally, obtaining multi-view image data of the working face, extracting features from the multi-view image data, and determining the matching relationship of feature points in combination with temporal consistency constraints. Identifying the changed area of the working face according to the matching relationship and the point cloud difference result includes: obtaining multi-view image data of the working face, performing camera calibration based on the Zhang Zheng-you calibration method and performing image enhancement preprocessing in combination with the Retinex theory, extracting multi-scale feature points from the preprocessed image using the SIFT algorithm, calculating the local structure tensor to evaluate the quality of feature points, constructing a spatial distribution uniformity constraint according to the distance between feature points, and performing optimized screening on feature points; calculating the feature appearance similarity and geometric similarity for the optimized and screened feature points, dynamically adjusting the similarity weight coefficient according to the stability of feature points, verifying the geometric consistency of feature points using the RANSAC algorithm, determining the matching relationship of feature points in combination with temporal motion continuity and multi-view epipolar geometry constraints, and calculating the confidence of feature point matching based on the matching relationship; performing three-dimensional reconstruction on the multi-view image data to obtain point cloud data, performing projection transformation registration on the point cloud data, using adaptive threshold segmentation to segment the changed area, calculating the feature matching probability according to the matching relationship, calculating the point cloud difference probability in combination with the point cloud difference, determining the fusion weight coefficient based on the confidence of feature point matching and the point cloud density, and fusing the feature matching probability and the point cloud difference probability according to the fusion weight coefficient to identify the changed area of the working face.

[0067] Exemplarily, in the image acquisition and preprocessing stage, first, a collection device is arranged at different positions around the working face using multiple industrial cameras. The cameras use CMOS sensors with a resolution of 1280×1024 and low-distortion lenses with a focal length of 12 mm, and the field of view angle is 61°×48°. The baseline distance between the cameras is 0.8 meters, and the difference in the camera coverage angles is about 30°-45° to ensure full coverage and multi-view observation of the working face. The camera acquisition frequency is set to 5 frames per second, and the exposure time is automatically adjusted within the range of 5-15 milliseconds according to the lighting conditions of the working face.

[0068] After obtaining the images, camera calibration is performed using the Zhang Zheng-you calibration method. A black and white checkerboard calibration board with a 9×11 grid is used in the calibration process, and the size of each grid is 30×30 mm. 20 groups of calibration images are collected at different distances and angles, and the camera internal parameter matrix and distortion coefficient are calculated through the calibration algorithm. In the evaluation of camera calibration accuracy, the reprojection error is controlled within 0.15 pixels, and the relative error of three-dimensional measurement is less than 0.5%.

[0069] After image acquisition, there are problems such as uneven illumination and insufficient contrast. The Retinex theory is used for image enhancement preprocessing. Specifically, the multi-scale Retinex algorithm is used, and three Gaussian kernel sizes are set to 15, 80, and 200 pixels respectively, with corresponding weights of 0.3, 0.4, and 0.3. By calculating the illumination component and reflection component of the image, the illumination component is adjusted to improve the image quality, while the detailed information in the reflection component is retained. The dynamic range of the processed image is expanded by 30%, and the visibility of details in low-illumination areas is increased by 65%, which is crucial for subsequent feature extraction.

[0070] The SIFT algorithm is used to extract multi-scale feature points from the preprocessed image. During the SIFT feature extraction process, the number of Gaussian pyramid layers is set to 4, the number of scales within each layer is set to 5, the feature point detection threshold is set to 0.04, and the edge response threshold is set to 10. For an image with a resolution of 1280×1024, usually 2000 - 3000 SIFT feature points can be extracted. The SIFT feature descriptor is represented by an 8-direction gradient histogram within a 4×4 grid, generating a 128-dimensional feature vector.

[0071] To evaluate the quality of feature points, the local structure tensor of the area near each feature point is calculated. The specific method is to calculate the average of the outer product of the image gradients within a window of 31×31 pixels around the feature point to obtain the structure tensor. The corner response value, anisotropy index, and stability index are calculated based on the eigenvalues of the structure tensor. The corner response value determines the significance of the feature point, the anisotropy index evaluates the feature discrimination ability, and the stability index evaluates the stability of the feature point under illumination and perspective changes.

[0072] Based on the results of feature point quality evaluation and spatial distribution, feature point optimization and screening are carried out. First, the feature points are sorted in descending order according to the quality score, and the top 60% of high-quality feature points are selected. Then, a spatial distribution uniformity constraint is constructed based on the distances between feature points. Specifically, the image is divided into a grid of 80×64, with each grid size of 16×16 pixels, and the number of feature points within each grid is calculated. For grids with overly dense feature points (more than 5 feature points), the top 3 feature points with the highest quality are retained; for grids without feature points, the feature point detection threshold is appropriately reduced for re-detection. In this way, the screened feature points not only ensure quality but also maintain spatial distribution uniformity, which is beneficial for subsequent accurate matching.

[0073] Multi-view matching is performed on the optimized and screened feature points. The matching process considers two aspects: feature appearance similarity and geometric similarity. Appearance similarity is measured by calculating the Euclidean distance between feature descriptors, and the ratio test method is used to determine reliable matches, with the ratio threshold set to 0.75. Geometric similarity is calculated by considering the similarity degree of the neighborhood structures of feature points, including local shape and scale ratio relationships.

[0074] Dynamically adjust the weight coefficients of appearance similarity and geometric similarity according to the stability index of feature points. For feature points with high stability (stability index greater than 0.8), the appearance similarity weight is set to 0.7 and the geometric similarity weight is set to 0.3; for feature points with medium stability (stability index 0.5 - 0.8), both weights are 0.5; for feature points with low stability (stability index less than 0.5), the appearance similarity weight is set to 0.3 and the geometric similarity weight is set to 0.7. This dynamic weight adjustment strategy improves the accuracy of feature point matching.

[0075] Use the RANSAC algorithm to verify the geometric consistency of feature points. The number of initial matching pairs is usually 500 - 1500 pairs, the number of iterations of the RANSAC algorithm is set to 2000, the inlier threshold is set to 2.0 pixels, and the confidence level is set to 0.99. Feature point matching pairs that meet the geometric consistency are obtained through screening by the RANSAC algorithm, and the number is approximately 60% - 80% of the initial matching pairs.

[0076] To further improve the matching accuracy, optimize the matching relationship by combining temporal motion continuity and multi-view epipolar geometry constraints. The temporal motion continuity constraint requires that the displacement change of the same feature point on consecutive time frames conforms to smoothness, and the change rate threshold is set to 20 pixels / second. The multi-view epipolar geometry constraint is calculated based on the fundamental matrix and requires that the matching point pairs satisfy the epipolar line constraint, and the deviation threshold is set to 1.5 pixels.

[0077] Calculate the confidence of feature point matching based on the optimized matching relationship. The confidence calculation considers the appearance similarity score, geometric consistency score, and temporal continuity score, and the weights of the three are 0.3, 0.4, and 0.3 respectively. The confidence score ranges from 0 to 1, and a high score indicates high matching reliability.

[0078] Perform three-dimensional reconstruction on multi-view image data to obtain point cloud data. The reconstruction uses a structured light-assisted multi-view stereo vision method. First, sparse point clouds are obtained through feature matching and triangulation, and then dense point clouds are obtained through a local region growth strategy.

[0079] Perform projection transformation registration on the point cloud data to analyze the changes in the point cloud data collected at different time points. In the registration process, first extract the plane features and edge features in the point cloud as registration markers, and then use the iterative closest point (ICP) algorithm for fine registration, and control the registration error within 2 mm.

[0080] After registration is completed, an adaptive threshold is used to segment the changed area. The adaptive threshold is dynamically calculated based on the local distribution characteristics of the point cloud, considering the point cloud density, noise level, and change amplitude. Specifically, the point cloud space is divided into grid cells of 50×50×50 cm. The mean and standard deviation of the height difference of the point cloud are calculated within each cell, and the mean plus twice the standard deviation is used as the change threshold for that cell. For areas with low point cloud density, the threshold tolerance is appropriately increased; for areas with high point cloud density, the threshold is decreased to improve sensitivity.

[0081] Calculate the feature matching probability based on the feature matching results, which represents the probability that the area has changed based on image features. The calculation of the feature matching probability considers three factors: the matching status of feature points, the confidence of the match, and the density of feature points. Areas where the feature point matching rate within the region is less than 50% or the match confidence is less than 0.6 are determined as high-change-probability regions.

[0082] Calculate the point cloud differential probability in combination with the point cloud differential results, which represents the probability that the area has changed based on the point cloud data. The point cloud differential probability is positively correlated with the height difference of the point cloud. Areas where the height difference exceeds the adaptive threshold have a high change probability.

[0083] Determine the fusion weight coefficient based on the feature match confidence and the point cloud density. When the feature match confidence is high (greater than 0.8) and the point cloud density is low, the weight of the feature matching probability is relatively large (0.7 - 0.9); when the feature match confidence is low and the point cloud density is high, the weight of the point cloud differential probability is relatively large (0.7 - 0.9); when the two indicators are balanced, the weights of the two probabilities are similar (0.4 - 0.6).

[0084] Finally, the feature matching probability and the point cloud differential probability are weighted and fused according to the fusion weight coefficient to obtain the comprehensive determination result of the changed area of the working face. The threshold of the fusion result is set to 0.65, and areas exceeding the threshold are marked as changed areas. Through this method of multi-source information fusion, the system can accurately identify the changed area of the working face, providing a reliable basis for subsequent feature extraction of the working face and optimization of the drilling layout.

[0085] The feature extraction of this method based on the SIFT algorithm combined with the spatial distribution uniformity constraint ensures the high quality and comprehensive coverage of feature points. The feature matching and dynamic weight adjustment strategies under multiple constraints improve the matching accuracy. The method for identifying the changed area by combining 3D reconstruction and point cloud analysis integrates image features and geometric information, enhancing the accuracy and robustness of the identification of the changed area of the working face.

[0086] In the second aspect of the embodiments of the present invention, an electronic device is provided, including: a processor; a memory for storing instructions executable by the processor; wherein, the processor is configured to call the instructions stored in the memory to execute the method described above.

[0087] In a third aspect of the embodiments of the present invention, a computer-readable storage medium is provided, on which computer program instructions are stored, and when the computer program instructions are executed by a processor, the foregoing method is implemented.

[0088] The present invention may be a method, an apparatus, a system, and / or a computer program product. The computer program product may include a computer-readable storage medium on which computer-readable program instructions for executing various aspects of the present invention are carried.

[0089] Finally, it should be noted that: the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it; although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that: they can still modify the technical solutions recorded in the foregoing embodiments, or perform equivalent replacements on some or all of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A multi-arm drilling rig adaptive drilling trajectory planning method based on vision navigation, characterized in that: include: Acquire multi-view image data of the working face, perform feature extraction on the multi-view image data, and determine the matching relationship of feature points in combination with temporal consistency constraints. Identify the working face change area based on the matching relationship and point cloud difference results, generate a topologically optimized drilling grid based on the working face change area, and determine the target drilling position set and drilling direction set; calculate the drill arm task allocation cost based on the target drilling position set and drilling direction set, and the task allocation cost includes the drill arm movement distance and the drilling direction adjustment amount; group the drilling holes based on the task allocation cost, generate spatiotemporal conflict constraints for each group of drilling holes, and iteratively optimize to determine the multi-arm collaborative operation sequence; in the process of executing the multi-arm operation sequence, use structured light projection technology to monitor the micro-deformation state of the rock mass during the drilling process, use the visual phase shift method to obtain the deformation field of the rock mass surface, establish the rock mass response characteristics containing elastic deformation, plastic deformation and abnormal deformation information, and use the progressive B-spline trajectory adjustment algorithm to adjust the drilling trajectory in real time based on the rock mass response characteristics.

2. The method according to claim 1, characterized in that A drilling grid with topology optimization is generated based on the working surface change area, and a target drilling position set and a drilling direction set are determined, including: constructing a grid density function based on the working surface change area, the grid density function comprehensively considering local curvature, change area distance function and structural importance metric, applying a topology optimization algorithm to solve the optimal drilling layout, the topology optimization algorithm adopts a hybrid strategy combining an improved SIMP method and a genetic algorithm; performing local and global evaluations on the optimal drilling layout, the local evaluation including energy utilization, and the global evaluation including grid uniformity, iteratively optimizing the problem area based on the evaluation results to obtain a target drilling position set and a drilling direction set.

3. The method according to claim 2, characterized in that The topology optimization algorithm adopts a hybrid strategy combining an improved SIMP method with a genetic algorithm, including: the improved SIMP method introduces an adaptive penalty factor and a local sensitivity correction, the adaptive penalty factor is dynamically adjusted with the number of iterations, and the local sensitivity correction refines the grid of the design domain based on the stress distribution gradient; the optimization result of the improved SIMP method is converted into the initial population of the genetic algorithm, the genetic algorithm adopts a real number encoding method, each chromosome contains both drilling position information and drilling direction information, and is designed with a crossover operator and an adaptive mutation operator, the crossover operator contains a local perturbation term, and the variable asynchronous length of the adaptive mutation operator decays with algebra, a fitness objective function including energy utilization and grid uniformity is constructed, and population iterative optimization is performed based on the fitness objective function.

4. The method according to claim 1, wherein The drill arm task allocation cost is calculated based on the target drilling position set and the drilling direction set, and the task allocation cost includes the drill arm movement distance and the drilling direction adjustment amount; drilling is grouped based on the task allocation cost, and spatiotemporal conflict constraints are generated for each group of drilling. The iterative optimization determines the multi-arm collaborative operation sequence, including: constructing a task allocation cost function for a single drilling hole, and the task allocation cost function includes the drill arm movement distance cost and the drilling direction adjustment cost. The drill arm movement distance cost is calculated based on the Euclidean distance between the current position of the drill arm and the target drilling position, and the drilling direction adjustment cost is calculated based on the Frobenius matrix of the current posture matrix of the drill arm and the target drilling direction matrix. norm calculation, and introduce a sequence continuity term based on feature distance; construct an inter-borehole similarity matrix based on the task allocation cost function, divide the boreholes into multiple initial groups using a spectral clustering algorithm, generate spatiotemporal conflict constraints according to the drill arm positions, and the spatiotemporal conflict constraints are determined based on the safety distance threshold between the drill arms; construct a global optimization objective function including the task allocation cost term, the spatiotemporal conflict constraint term and the sequence continuity term, and solve the global optimization objective function using a simulated annealing algorithm. Generate a neighborhood solution by exchanging the drilling sequence within a group and across groups, calculate the acceptance probability based on the change of the global optimization objective function and the temperature parameter, and obtain the multi-arm collaborative operation sequence through iterative optimization.

5. The method according to claim 4, characterized in that The simulated annealing algorithm is used to solve the global optimization objective function, including: designing an adaptive neighborhood generation strategy, dynamically adjusting the probabilities of intra-group exchange and cross-group exchange according to the current iteration stage, the intra-group exchange probability is inversely proportional to the temperature, and the cross-group exchange probability is proportional to the temperature; when the improvement rate of the global optimization objective function of N consecutive iterations is lower than a preset threshold, expanding the neighborhood search range and introducing a trajectory memory mechanism; constructing an acceptance probability calculation method based on historical information, recording the improvement degree of the objective function of each accepted solution, adaptively adjusting the acceptance threshold according to the historical improvement rate, and introducing a second-order perturbation term based on Taylor expansion, the amplitude of the second-order perturbation term is proportional to the curvature of the current solution space; adopting a segmented cooling strategy, using different cooling coefficients in different temperature intervals, and the cooling coefficient of the initial temperature interval is greater than the cooling coefficient of the subsequent temperature intervals; when a better solution cannot be obtained after several consecutive iterations, the reheating temperature parameters are adaptively determined based on the distribution of historical optimal solutions, and the characteristics of the historical optimal solution are maintained during the reheating process.

6. The method according to claim 1, characterized in that The structured light projection technology is used to monitor the micro-deformation state of the rock mass during the drilling process, the visual phase shift method is used to obtain the deformation field of the rock mass surface, and the rock mass response characteristics including elastic deformation, plastic deformation and abnormal deformation information are established. According to the rock mass response characteristics, the progressive B-spline trajectory adjustment algorithm is used to adjust the drilling trajectory in real time, including: using structured light projection technology to project a sinusoidal grating pattern onto the rock mass surface, using the visual phase shift method to obtain four groups of images with different phases and calculate the phase difference, expanding the phase difference to obtain the absolute phase, and calculating the three-dimensional deformation field data of the rock mass surface based on the absolute phase and system calibration parameters; extracting the macroscopic deformation field description features and local deformation field description features of the three-dimensional deformation field data. The deformation anomaly index and time-series deformation trend characteristics are used to construct a multi-scale deformation feature space. The RANSAC algorithm is used to set an evaluation unit to iteratively identify abnormal deformation areas in the multi-scale deformation feature space. The elastic deformation parameters, plastic deformation parameters and abnormal deformation parameters are calculated according to the abnormal deformation areas, and the fuzzy reasoning method is used to obtain the rock mass response characteristics. According to the rock mass response characteristics, a multi-constraint objective function including the deformation field gradient constraint function, the drilling process parameter constraint function and the equipment kinematic constraint function is constructed. The sliding time window strategy is used to locally optimize the multi-constraint objective function. The position of the B-spline trajectory control point is adjusted based on the local optimization results, and the drilling trajectory is optimized and adjusted in real time.

7. The method according to claim 6, characterized in that According to the rock mass response characteristics, a multi-constraint objective function is constructed, a sliding time window strategy is used for local optimization, and the position of the B-spline trajectory control point is adjusted based on the local optimization result. The real-time optimization and adjustment of the drilling trajectory include: the deformation field gradient constraint function is calculated based on the deformation gradient tensor, the drilling process parameter constraint function integrates energy consumption terms, wear terms and time terms, and the equipment kinematic constraint function includes curvature constraints and speed continuity constraints; according to the rock mass response characteristics, a multi-level optimization trigger mechanism based on time triggering, event triggering and performance triggering is established, and the optimization priority of the multi-level optimization trigger mechanism is dynamically associated with the rock mass risk level; when the optimization is triggered, a sliding time window is defined, and the length of the sliding time window is determined according to the rock mass response characteristics. The rock mass response characteristics change adaptively, and the trajectory is represented as a B-spline curve within the sliding time window. The control points of the B-spline curve include fixed control points, transition control points and optimizable control points. A time-aware degradation optimization strategy is designed and calculated, and the multi-constraint objective function is locally optimized using a sequential quadratic programming method. Physical constraints are constructed through the geometric relationship of control points using the B-spline characteristics. The positions of the B-spline trajectory control points are adjusted based on the local optimization results, and a mixing coefficient calculation method based on the influence domain of the control points is used for the transition control points to achieve smooth transition of the trajectory. The optimization degree is automatically adjusted according to the time margin calculated by the sequential quadratic programming method, and an adaptive mixing strategy of partial optimization results and the old trajectory is adopted to complete the trajectory adjustment within a preset time.

8. The method according to claim 1, characterized in that Acquire multi-view image data of the working surface, perform feature extraction on the multi-view image data, determine the matching relationship of feature points in combination with temporal consistency constraints, and identify the changed area of the working surface according to the matching relationship and point cloud difference results, including: acquiring multi-view image data of the working surface, performing camera calibration based on Zhang Zhengyou calibration method and performing image enhancement preprocessing in combination with Retinex theory, extracting multi-scale feature points from the preprocessed image using SIFT algorithm, calculating local structure tensor to evaluate the quality of feature points, constructing spatial distribution uniformity constraints according to the distance between feature points, and optimizing and screening feature points; calculating feature appearance similarity and geometric similarity of the optimized and screened feature points, and performing image enhancement preprocessing according to the stability of feature points. Qualitatively and dynamically adjust the similarity weight coefficient, use the RANSAC algorithm to verify the geometric consistency of the feature points, combine the temporal motion continuity and multi-view epipolar geometric constraints to determine the matching relationship of the feature points, and calculate the feature point matching confidence based on the matching relationship; perform three-dimensional reconstruction of the multi-view image data to obtain point cloud data, perform projection transformation and registration on the point cloud data, use adaptive thresholding to segment the changed area, calculate the feature matching probability based on the matching relationship, calculate the point cloud difference probability based on the point cloud difference, determine the fusion weight coefficient based on the feature matching confidence and point cloud density, fuse the feature matching probability and point cloud difference probability according to the fusion weight coefficient, and identify the changed area of the working surface.

9. An electronic device, characterized in that: include: processor; A memory for storing processor-executable instructions; wherein the processor is configured to call the instructions stored in the memory to execute the method according to any one of claims 1 to 8.

10. A computer-readable storage medium having computer program instructions stored thereon, characterized in that: When the computer program instructions are executed by a processor, the method according to any one of claims 1 to 8 is implemented.

Images