A blasting parameter optimization method and system for blasting fragmentation control

Abstract

The present application relates to the technical field of blasting control, in particular to a blasting parameter optimization method and system for blasting block size control. The method comprises the following steps: extracting the rock mass structure plane network characteristics through structure plane digital identification and topological analysis; analyzing the rock mass wave impedance partition characteristics by using the measurement while drilling inversion and wave impedance clustering; obtaining the aperture grouping characteristics based on the cluster of blast hole diameter, combining the rock mass structure plane network characteristics, the rock mass wave impedance partition characteristics and the aperture grouping characteristics, calculating the rock mass explosibility comprehensive index through the attention enhanced deep neural network; expanding the block size distribution sample through the generative adversarial network, combining the simulated block size distribution characteristics and the actual block size distribution characteristics to construct the blasting area space block size correction field; according to the blast hole grid executability constraint characteristics, the rock mass explosibility comprehensive index and the blasting area space block size correction field, solving the optimal blasting parameter Pareto solution set. The present application realizes the accurate design of blasting parameters, which is conducive to improving the blasting effect.

Description

Technical Field

[0001] This invention relates to the field of blasting control technology, specifically to a method and system for optimizing blasting parameters for controlling the size of blasted blocks. Background Technology

[0002] Blasting is one of the core procedures in engineering construction. The rationality of blasting parameters, such as borehole spacing, charge quantity, and detonation time difference, directly determines the blasting effect, construction cost, and feasibility. Currently, blasting parameter optimization methods generally suffer from the following technical shortcomings, making it difficult to meet the demands of modern blasting engineering for high precision, low cost, and high feasibility: The simulation accuracy is low due to incomplete consideration of rock mass characteristics: Existing methods mostly use a single rock mass parameter (such as compressive strength) to describe rock mass characteristics, without fully integrating the characteristics of rock mass structural surface network and the spatial distribution characteristics of wave impedance, and ignoring the influence of rock mass heterogeneity on blasting energy transfer and rock mass fragmentation; at the same time, the simulation of explosive detonation energy release mostly uses a fixed rate model, which cannot accurately capture the spatiotemporal energy release law of the explosive from ignition, detonation to the end of the reaction, resulting in a large deviation between the block size simulation and reality.

[0003] Insufficient fidelity in block size prediction and lack of error correction: Existing block size simulations mostly use finite element method (FEM) or discrete element method (DEM), which cannot achieve cross-scale organic connection between macroscopic stress propagation, microscopic particle motion and block size distribution; and due to the influence of small sample and uneven sample distribution, the deviation between simulated block size and actual block size cannot be effectively corrected, resulting in parameter optimization based on simulation results deviating from engineering reality.

[0004] The lack of construction feasibility constraints makes it difficult to implement optimization solutions: Existing optimization methods focus on optimizing blasting effects and costs, without fully considering the kinematic characteristics of on-site construction equipment and drill bit inventory, or the accessibility of borehole grid layout and the efficiency of drilling rig coordination. As a result, the optimized blasting parameters often have problems such as equipment incompatibility, unreachable construction paths, and frequent drill bit replacements, making them difficult to apply directly to on-site construction.

[0005] Multi-objective optimization has limitations, and the weight setting is subjective: existing multi-objective optimization methods mostly use a single algorithm (such as genetic algorithm) to solve the problem, which cannot effectively handle the conflicting relationships between multiple objectives such as blasting effect, cost, and construction feasibility; moreover, the weight setting in the scheme selection process relies heavily on human subjective experience and lacks objective basis, resulting in insufficient comprehensive applicability of the optimal scheme.

[0006] Therefore, there is an urgent need for a blasting parameter optimization method that can integrate multi-source rock mass and construction characteristics, achieve accurate simulation across scales, balance multi-objective optimization and construction feasibility, and set objective weights, in order to overcome the shortcomings of existing technologies. Summary of the Invention

[0007] To address the shortcomings of existing methods and the needs of practical applications, and to solve the aforementioned problems, this invention provides a method for optimizing blasting parameters for controlling the size of blasted blocks, comprising the following steps: The rock mass structural surface network features are extracted through structural surface digital recognition and topology analysis. The rock mass wave impedance zoning features are analyzed using measurement-while-drilling inversion and wave impedance clustering. Hole diameter grouping features are obtained based on borehole diameter clustering. Combining the rock mass structural surface network features, the rock mass wave impedance zoning features, and the hole diameter grouping features, an attention-enhanced deep neural network is used to calculate the comprehensive rock mass blastability index. The block size distribution samples are expanded using a generative adversarial network, and a spatial block size correction field for the blasting area is constructed by combining simulated block size distribution features and actual block size distribution features. Based on the executability constraint features of the borehole grid, the comprehensive rock mass blastability index, and the spatial block size correction field for the blasting area, the Pareto solution set of the optimal blasting parameters is solved.

[0008] Optionally, the step of extracting the network features of rock mass structural surfaces through structural surface digital recognition and topology analysis includes the following steps: Based on the point cloud dataset of the rock mass surface, the geometric parameters of the structural surface are obtained by segmenting the point cloud using the region growing method; a topological analysis network of the rock mass structural surface is constructed according to the geometric parameters of the structural surface to obtain the network features of the rock mass structural surface.

[0009] Optionally, the step of analyzing the rock mass impedance zoning characteristics using measurement-while-drilling inversion and impedance clustering includes the following steps: A wave impedance inversion model is constructed using segmented drilling parameter datasets to obtain wave impedance profile data for all boreholes. Wave impedance clustering and partitioning are performed on the wave impedance profile data to generate a wave impedance partition map of the blast area, thereby obtaining the wave impedance partition characteristics of the rock mass.

[0010] Optionally, obtaining aperture grouping features based on borehole diameter clustering includes the following steps: Hierarchical clustering algorithm is used to cluster the designed borehole diameters; based on the objective function of the number of drill bit replacements, a globally optimal scheduling scheme is obtained through multi-drilling rig collaborative scheduling optimization; combining the clustering results and the globally optimal scheduling scheme, the borehole diameter grouping characteristics are obtained.

[0011] Optionally, the step of combining the rock mass structural surface network characteristics, the rock mass impedance zoning characteristics, and the aperture grouping characteristics to calculate the comprehensive rock mass blastability index using an attention-enhanced deep neural network includes the following steps: Construct an attention-enhanced deep neural network; and obtain a comprehensive rock mass blastability index by weightedly fusing the rock mass structural surface network features, the rock mass wave impedance zoning features, and the aperture grouping features through the attention-enhanced deep neural network.

[0012] Optionally, the step of generating adversarial networks to expand block size distribution samples and constructing a burst zone spatial block size correction field by combining simulated block size distribution features and actual block size distribution features includes the following steps: By generating adversarial networks to expand the block size distribution samples, a block size distribution sample set is obtained. Based on the block size distribution sample set, the block size error between the simulated block size distribution characteristics and the actual block size distribution characteristics is calculated. Using the two-dimensional spatial coordinates of the blast zone as the independent variable, the Kriging interpolation algorithm is used to interpolate the calculated block size error to each grid position in the entire blast zone, thereby obtaining the error value corresponding to each spatial coordinate point and forming the blast zone spatial block size correction field.

[0013] Optionally, extracting the simulated block size distribution features includes the following steps: By combining chemical reaction kinetics to reconstruct the time-space curve of explosive energy release, and through cross-scale coupling of population balance model-finite element-discrete element, the blasting process is simulated from the overall blasting zone to rock particles to obtain the explosive energy release law and block size distribution characteristics.

[0014] Optionally, extracting the actual block size distribution features includes the following steps: A U-Net segmentation model is constructed, and the equivalent grain size data of each rock block is obtained using the U-Net segmentation model; block size gradation statistics are performed based on the equivalent grain size data to obtain the actual block size characteristic parameters.

[0015] Optionally, the step of solving for the optimal Pareto solution set of blasting parameters based on the executability constraint characteristics of the borehole grid, the comprehensive rock mass blastability index, and the spatial block size correction field of the blasting zone includes the following steps: Combining the executability constraint characteristics of the borehole grid, the comprehensive index of rock mass blastability, and the spatial block size correction field of the blasting zone, a multi-objective optimization function is constructed. The multi-objective optimization function is solved by a non-dominated sorting genetic algorithm with an elite strategy to obtain the Pareto front solution set. Then, the optimal solution is selected from the solution set by TOPSIS multi-attribute decision-making.

[0016] This invention reconstructs the spatiotemporal release curve of explosive detonation energy by integrating the rock mass structural surface network and the spatial partitioning characteristics of wave impedance. It uses PBM-FEM-DEM cross-scale coupling to simulate the rock mass fracturing process, and combines GAN data enhancement and spatial correction field to correct block size prediction errors. It incorporates borehole grid construction constraints, borehole diameter grouping, and multi-drilling rig collaborative scheduling, and solves the Pareto front solution set based on the NSGA-II algorithm. It also combines entropy weight-TOPSIS to screen the optimal blasting parameters. This invention solves the shortcomings of existing methods, such as incomplete consideration of rock mass characteristics, low block size prediction fidelity, difficulty in implementing optimization schemes, and subjective weight setting. It achieves accurate optimization of blasting parameters, taking into account blasting effect, construction cost, and feasibility, effectively improving the comprehensive benefits of blasting projects, and is applicable to various blasting construction scenarios.

[0017] Secondly, to efficiently execute the blasting parameter optimization method for blasting block size control provided by this invention, this invention also provides a blasting parameter optimization system for blasting block size control, comprising: an input device, an output device, a processor, and a memory, wherein the input device, output device, processor, and memory are interconnected, and the memory stores program instructions for the blasting parameter optimization method for blasting block size control. The blasting parameter optimization system for blasting block size control of this invention has a compact structure and stable performance, and can stably execute the blasting parameter optimization method for blasting block size control provided by this invention, further improving the overall applicability and practical application capability of this invention. Attached Figure Description

[0018] Figure 1 A flowchart of a blasting parameter optimization method for controlling blasting block size is provided in an embodiment of the present invention; Figure 2 This is a framework diagram of a blasting parameter optimization system for controlling blasting block size, provided in an embodiment of the present invention. Detailed Implementation

[0019] Specific embodiments of the present invention will now be described in detail. It should be noted that the embodiments described herein are for illustrative purposes only and are not intended to limit the invention. In the following description, numerous specific details are set forth in order to provide a thorough understanding of the invention. However, it will be apparent to those skilled in the art that these specific details are not necessary to practice the invention. In other instances, well-known circuits, software, or methods have not been specifically described to avoid obscuring the invention.

[0020] Throughout this specification, references to an embodiment, example, or illustration mean that a particular feature, structure, or characteristic described in connection with that embodiment or example is included in at least one embodiment of the invention. Therefore, phrases appearing in various places throughout the specification, such as "in one embodiment," "in an embodiment," "an example," or "an illustration," do not necessarily refer to the same embodiment or example. Furthermore, specific features, structures, or characteristics can be combined in any suitable combination and / or sub-combination in one or more embodiments or examples. Moreover, those skilled in the art will understand that the illustrations provided herein are for illustrative purposes and are not necessarily drawn to scale.

[0021] Please see Figure 1 To address the above problems, this invention provides a method for optimizing blasting parameters for controlling the size of blasting blocks, such as... Figure 1 As shown, in one embodiment, the method includes the following steps: S1. Extract the network features of rock mass structural surfaces through digital identification and topology analysis of structural surfaces.

[0022] In this embodiment, the extraction of rock mass structural surface network features through structural surface digital recognition and topology analysis includes the following steps: S11. Based on the point cloud dataset of the rock mass surface, the geometric parameters of the structural surface are obtained by segmenting the point cloud using the region growing method.

[0023] The raw point cloud data acquired by 3D laser scanning is susceptible to environmental interference, containing a large number of outliers, duplicates, and uneven density. First, denoising is required to effectively remove scanning noise and environmental interference points. Second, multi-station point cloud registration is necessary. Since a single-station scan cannot cover the entire blasting area, the ICP (Iterative Closest Point) algorithm is used to align the multi-station scan data. Using the first station's point cloud as a reference, feature point pairs in overlapping areas between stations are selected. By iteratively optimizing the rotation matrix and translation vector, the Euclidean distance between feature point pairs is minimized, achieving seamless stitching of multi-station point clouds. Finally, resampling is performed using a uniform voxel resampling method. The point cloud space is divided into voxels of fixed size, with one most representative point retained within each voxel, homogenizing the point cloud density and reducing data redundancy. The final result is a well-organized, clean, and uniformly dense point cloud dataset of the rock mass surface.

[0024] Furthermore, since rock mass structural surfaces often exhibit continuous planar characteristics, planar clustering is achieved based on the similarity of point cloud normal vectors: First, key parameters are set, including the normal vector angle threshold, the minimum number of points to grow in the region, and the maximum number of points. The normal vector angle threshold is set to 5° (which can be adjusted according to the flatness of the rock mass structural surface; the flatter the structural surface, the smaller the threshold; for rough structural surfaces, it can be relaxed to 8°). The minimum number of points is set to 50 (to ensure that the segmented region is an effective structural surface rather than an isolated cluster of points), and the maximum number of points is reasonably set according to the scale of the blast zone structural surfaces. Then, seed points are selected. Points with stable point cloud normal vectors and uniform surrounding point cloud density are selected as initial seed points to avoid selecting edge noise points. Finally, region growth is initiated. Starting from the seed point, all points in its neighborhood are searched, and the angle between the normal vector of the neighboring point and the seed point is calculated. If the angle is less than a set threshold, the point is included in the current growth region and used as a new seed point to continue growth. This process is repeated until no new points can be included, and finally multiple continuous planar regions are segmented, which are potential rock mass structure surfaces. At the same time, isolated planar clusters that are too small or have irregular shapes are removed.

[0025] The potential structural surfaces obtained from point cloud segmentation may include pseudo-structural surfaces formed by uneven areas on the rock mass surface and locally damaged areas. Accurate verification using field-recorded data is necessary to ensure the accuracy of structural surface identification. Field-recorded data mainly includes manually measured structural surface attitudes (dip α, dip β), structural surface spacing, infill type, and distribution range. First, the attitude data of each field-recorded structural surface is converted into corresponding plane normal vector parameters, which are then compared one by one with the normal vectors of the potential structural surfaces obtained from point cloud segmentation. A normal vector matching threshold is set. If the angle between the normal vector of a potential structural surface and the normal vector of a field-recorded structural surface is less than 3°, and their spatial overlap is greater than 70%, then the potential structural surface is determined to be a real rock mass structural surface. If a potential structural surface lacks corresponding field-recorded data support, or the normal vector angle exceeds the threshold or the spatial overlap is insufficient, it is determined to be a pseudo-structural surface (such as non-penetrating fractures with uneven surfaces or irregular planes formed by rock weathering and erosion) and is discarded, ultimately resulting in a set of real and reliable rock mass structural surfaces.

[0026] Furthermore, the geometric parameters of the structural surfaces, including trace length L, aperture d, and attitude (α, β), are extracted from the verified real rock mass structural surfaces. The core of calculating the trace length L is to obtain the actual extension length of the structural plane on the rock mass surface. First, the boundary point cloud of each structural plane is extracted using a point cloud contour extraction algorithm to construct a two-dimensional contour graphic of the structural plane. The minimum bounding rectangle method is used to fit the two-dimensional contour of the structural plane, and the length of the longer side of the minimum bounding rectangle is taken as the trace length L of the structural plane. If the contour of the structural plane is irregular, it needs to be corrected by combining the structural plane spacing data recorded on site to ensure that the trace length calculation error is less than 5%. The aperture d is the gap width of the structural plane, reflecting the degree of development of the structural plane. It is calculated by statistically analyzing the distance difference between the point clouds on both sides of the structural plane in the normal direction of the structural plane. 10-15 measurement points at different positions of the structural plane are selected, and the normal distance difference of each measurement point is calculated. The average value is taken as the aperture d of the structural plane, and abnormal measurement values ​​(such as excessively large or small distance differences caused by point cloud noise) are eliminated. The attitude (α, β) is obtained by plane normal vector transformation. Let the plane normal vector of the structural plane be... ,in α and β are the components of the normal vector on the x, y, and z axes, respectively. The inclination α is the angle between the projection of the normal vector onto the horizontal plane and the positive x-axis, and the tilt angle β is the angle between the normal vector and the horizontal plane.

[0027] S12. Construct a topological analysis network of the rock mass structural surface based on the geometric parameters of the structural surface to obtain the network characteristics of the rock mass structural surface.

[0028] Each verified real rock mass structural surface is defined as a node in the network, and each node corresponds to a set of geometric parameters (L, d, α, β) of a structural surface. The intersection relationship between two structural surfaces is defined as an edge in the network. If two structural surfaces have an intersection point in space, and the intersection point is located within the effective range of the two structural surfaces (not the edge endpoint), then it is determined that the two have a connection relationship, and there is an edge between the two nodes in the network; otherwise, there is no edge.

[0029] Then, the core topology parameters, including connectivity C and line density, are calculated. Dough density The connectivity rate C reflects the degree of connectivity of the structural surface network. The formula is C = number of connected structural surfaces / total number of structural surfaces. A connected structural surface is one that intersects with at least one other structural surface. A larger C value indicates stronger connectivity of the rock mass's structural surfaces, making it easier to fracture along these surfaces during blasting. Linear density... The formula reflects the linear distribution density of structural planes within the rock mass space and is calculated as follows: =Total trace length / Rock mass volume, where the total trace length is the sum of the trace lengths L of all structural surfaces, and the rock mass volume is calculated based on the three-dimensional geometric model of the blasting area. A higher value indicates a denser development of structural planes; surface density The formula reflects the area distribution density of structural planes within the rock mass space and is calculated as follows: = Total structural surface area / rock mass volume, where the total structural surface area is estimated using the trace length L and average aperture d of each structural surface. The value directly reflects the integrity of the rock mass; the larger the value, the worse the integrity of the rock mass.

[0030] Furthermore, an adjacency matrix is ​​used to mathematically represent network relationships, and core topological features are extracted through matrix analysis. First, an adjacency matrix A is constructed, where the number of rows and columns equals the total number of structural surfaces. The matrix elements... The definition is: if structural plane i and structural plane j have an intersecting and connected relationship, then If structural plane i and structural plane j do not have an intersecting or connected relationship, then At the same time, set the diagonal elements (The structural planes are not connected to themselves), and finally a binary adjacency matrix is ​​obtained.

[0031] Subsequently, through eigenvalue analysis of the adjacency matrix, the core topological features of the network were extracted: the clustering coefficient reflects the degree of local clustering in the structural surface network. The calculation formula for each node's clustering coefficient is: the actual number of edges between adjacent nodes / the maximum possible number of edges between adjacent nodes. The average clustering coefficient of all nodes is taken as the clustering coefficient of the entire network. The larger the clustering coefficient, the easier it is for the structural surfaces to form dense clusters in local areas, and the easier it is to form small rock blocks during blasting. The average path length reflects the average connectivity distance between nodes in the structural surface network. The calculation formula is the average of the shortest path lengths between all pairs of nodes. The shortest path length is solved by exponentiation of the adjacency matrix. The shorter the average path length, the stronger the connectivity between structural surfaces and the more uniform the transfer of explosive energy in the rock mass.

[0032] Furthermore, the characteristic parameters are categorized into two main types: geometric features and topological features. Geometric features include the trace length L, aperture d, dip direction α, and dip angle β of each structural surface. To reflect the overall geometric characteristics of the rock mass's structural surfaces, statistical values ​​(mean, standard deviation, maximum, and minimum) of each geometric parameter need to be calculated, such as the average trace length, aperture standard deviation, and dip direction distribution range. Topological features include connectivity C and linear density. areal density Parameters such as network clustering coefficient, average path length, number of connected components extracted from adjacency matrix, and proportion of isolated structural surfaces.

[0033] Subsequently, a structured feature vector is constructed, and dimension normalization is performed (normalizing all feature parameters to the 0-1 range to eliminate dimensional differences). Each structure corresponds to a basic feature sub-vector (containing its own geometric parameters), and the entire blast zone rock mass corresponds to a global feature vector (containing all geometric parameter statistics and topological feature parameters).

[0034] S2. Analyze the wave impedance zoning characteristics of the rock mass using measurement-while-drilling inversion and wave impedance clustering.

[0035] In this embodiment, the analysis of rock mass acoustic impedance zoning characteristics using measurement-while-drilling inversion and acoustic impedance clustering includes the following steps: S21. Construct a wave impedance inversion model using segmented drilling parameter datasets to obtain wave impedance profile data for all boreholes.

[0036] During down-the-hole drilling (DWT) measurements, the raw DWT data (ROP, WOB, torque T, and impact frequency f) will contain a large number of outliers due to factors such as rig start-up and shutdown, drill rod vibration, rock powder blockage in the borehole, and measurement sensor errors. Outliers must be removed first, and then, in conjunction with the rig operation log, data with sudden changes under special conditions such as rig start-up and shutdown, drill rod replacement, and stuck drill in the borehole (e.g., ROP suddenly rises to 0 or increases to more than 3 times the normal range, or WOB suddenly exceeds the rig's rated drilling pressure range) should be manually removed to avoid misjudging valid data. The data was then segmented. Considering the variation of rock mass mechanical properties with drilling depth, and taking into account the drilling efficiency and measurement accuracy of the down-the-hole drilling rig, each data segment was set to 0.5m (if the degree of rock weathering varies greatly, it can be adjusted to 0.3m / segment to ensure relatively uniform rock mass properties within each segment). Statistical analysis was performed on the ROP, WOB, and T data within each depth segment, calculating the mean and variance of each parameter. The mean is used for subsequent wave impedance inversion modeling, and the variance is used to evaluate the stability of the rock mass properties at that depth segment (the larger the variance, the more heterogeneous the rock mass properties are, and this needs to be closely monitored during the inversion process). Finally, a well-structured segmented drilling parameter dataset was output. Wave impedance (ρv, where ρ is rock mass density and v is longitudinal wave velocity) is a core parameter characterizing the density and elastic characteristics of rock mass, and directly determines the energy transfer efficiency of explosives. However, measuring wave impedance through on-site core sampling is costly and time-consuming, and the spatial representativeness of the core samples is limited. Based on the strong correlation between drilling parameters and wave impedance, an inversion model is constructed to achieve rapid and accurate calculation of wave impedance.

[0037] Specifically, based on field coring data, a multiple linear regression model of wave impedance ρv and drilling parameters is established: First, representative boreholes at different locations and depths in the blast zone are selected for field coring tests. The density ρ and P-wave velocity v of each core sample are obtained through laboratory testing, and the wave impedance ρv at the corresponding depth is calculated (ρv=ρ×v). At the same time, the mean values ​​of drilling parameters (ROP, WOB, T) at the corresponding depths of these coring locations are extracted, and a sample dataset is constructed (each sample contains three independent variables: ROP, WOB, and T, and one dependent variable: ρv). The number of samples should be no less than 50 to ensure the reliability of the model fitting.

[0038] Next, the ρv value of each borehole depth segment is inverted using a multiple linear regression model to obtain the wave impedance profile on the borehole trajectory. The mean values ​​of ROP, WOB, and T for each borehole and each depth segment are substituted into the fitted multiple linear regression model to calculate the predicted wave impedance value for each depth segment. Combined with the borehole depth coordinates, the predicted wave impedance values ​​for each depth segment are connected in depth order to form a wave impedance profile on a single borehole trajectory, clearly showing the variation law of wave impedance with borehole depth (e.g., the wave impedance value is lower in weathered rock strata and higher in intact rock strata). Finally, the wave impedance profile data of all boreholes are output.

[0039] S22. Perform wave impedance clustering and partitioning on the wave impedance profile data to generate a wave impedance partitioning map of the blast area and obtain the wave impedance partitioning characteristics of the rock mass.

[0040] The wave impedance profile of a single borehole can only reflect the changes in rock mass properties along the borehole trajectory and cannot reflect the spatial distribution differences of the rock mass in the entire blasting area. Therefore, it is necessary to classify the wave impedance data of all boreholes through clustering algorithms and combine spatial interpolation to realize the spatial partitioning of wave impedance in the blasting area, so as to provide a partitioning basis for subsequent explosive energy matching.

[0041] Specifically, the K-means clustering algorithm is used to cluster the ρv values ​​of all boreholes in the blasting area: First, the wave impedance inversion values ​​of all boreholes and all depth segments in the blasting area are collected to construct a wave impedance sample set for the blasting area; The core of K-means clustering is to determine the optimal number of clusters K, and the maximum silhouette coefficient criterion is used for screening. The silhouette coefficient is used to evaluate the rationality of the clustering results (the closer the silhouette coefficient is to 1, the better the clustering effect, the higher the similarity of samples of the same type and the lower the similarity of samples of different types). Based on engineering experience, the number of clusters K is usually taken as 3~5 (corresponding to weak integrity, medium integrity and strong integrity rock masses respectively. If the rock mass properties of the blasting area are very different, it can be adjusted to 5~6 categories).

[0042] During the clustering process, K wave impedance values ​​are randomly selected as initial cluster centers. The Euclidean distance from each sample to each cluster center is calculated, and the sample is assigned to the category of the nearest cluster center. Then, each cluster center is updated (the mean wave impedance of all samples in that category is taken as the new cluster center). The process of distance calculation, sample assignment, and cluster center update is repeated until the cluster centers no longer change or the change is less than 0.01 (convergence condition). After the clustering is completed, each category is verified by calculating the mean and standard deviation of wave impedance for each category. Abnormal clusters (such as categories with too few samples or abnormal wave impedance ranges) are removed. Finally, the clustering results of the wave impedance of the rock mass in the blasting area are obtained, and the rock mass integrity level corresponding to each wave impedance category is determined.

[0043] Furthermore, spatial interpolation (Kriging interpolation) is performed on the clustering results to generate a wave impedance zoning map of the blast area. First, a two-dimensional spatial coordinate system for the blast area is constructed (with a corner point of the blast area as the origin, and the x-axis and y-axis corresponding to two orthogonal directions in the horizontal direction, respectively), and the wave impedance clustering results of all boreholes are associated with the corresponding spatial coordinates. Using the Kriging interpolation algorithm, based on the wave impedance clustering data of known boreholes, the wave impedance category of all unknown locations in the blast area is predicted to generate a wave impedance zoning map of the blast area. In the zoning map, each zoning is labeled with the corresponding clustering category, the mean wave impedance (the average wave impedance of all samples in the zoning), the zoning area (calculated using the grid counting method, with the grid size consistent with the grid used for modeling the feasible region of the blast area), and the width of the transition zone with adjacent zoning (the transition zone is the gradual area between two zoning categories, and the width is usually 1~2m, adjusted according to the spatial continuity of the clustering results). The wave impedance value of the transition zone is processed as a linear gradual change to avoid abrupt changes at the zoning boundary and ensure that the zoning results conform to the actual distribution law of the rock mass.

[0044] After clustering and partitioning, a standardized feature dataset is formed. Specific output features include: ① Partition number: Each partition is numbered from high to low rock mass integrity level (or from large to small wave impedance value) (e.g., Z1, Z2, Z3...) for easy subsequent identification and retrieval; ② ρv range of each partition: The minimum and maximum values ​​of all wave impedance samples within each partition are statistically analyzed to clarify the wave impedance value range of that partition, intuitively reflecting the differences in density and elastic characteristics of the rock mass within each partition; ③ Spatial distribution coordinates: Using raster coordinates or Cartesian coordinates, the boundary coordinates and center coordinates of each partition are marked, clarifying the spatial location and distribution range of each partition within the blasting area, as well as its positional relationship with the blasting area boundary, obstacles, and borehole grid; ④ Integrity coefficient: Calculated based on wave impedance values, used to quantitatively characterize rock mass integrity, satisfying: ,in The measured wave impedance of the intact rock mass in the blasting area is the average wave impedance of the intact rock core obtained through on-site core sampling tests. The value range is 0~1. The closer the value is to 1, the stronger the integrity of the rock mass, the closer the wave impedance is to that of an intact rock mass, and the higher the energy transfer efficiency of the explosive.

[0045] S3. Based on the borehole diameter clustering, obtain the borehole diameter grouping features, and combine the rock mass structural surface network features, the rock mass wave impedance zoning features, and the borehole diameter grouping features to calculate the comprehensive rock mass blastability index through attention-enhanced deep neural network.

[0046] In another embodiment, obtaining aperture grouping features based on borehole diameter clustering includes the following steps: S311. Use hierarchical clustering algorithm to cluster the designed aperture.

[0047] The absolute value of the difference in borehole diameters is used as the distance metric. The smaller the distance, the higher the similarity between the two borehole diameters, and the more suitable they are to be grouped together. A reasonable distance threshold is then set. The determination of the distance threshold needs to consider the drill bit model and specifications, the allowable range of borehole diameter design error, and engineering construction experience. It is usually set to 5mm (if there are significant differences in drill bit model and specifications on site, it can be adjusted to 3~8mm). The core function of the distance threshold is to define the range of borehole diameters that can be shared with drill bits. When the distance between two borehole diameters is ≤ the distance threshold, it can be considered that the error is within the allowable range, and they can be used with the same model of drill bit. When the distance is > the distance threshold, they need to belong to different groups and be used with different models of drill bits.

[0048] The clustering process employs agglomerative hierarchical clustering (bottom-up clustering): Initially, each effective design aperture is considered as an independent cluster; then, the distance between all clusters is calculated (the average distance between all apertures within a cluster and all apertures within another cluster is taken as the inter-cluster distance), and the two clusters with the smallest distance ≤ the distance threshold are merged into a new cluster; the inter-cluster distance calculation and merging operation are repeated until the distance between all clusters is > the distance threshold, at which point the clustering stops, ultimately resulting in several aperture cluster groups G1, G2, ..., Gk (k is the number of cluster groups, determined by the aperture distribution characteristics and the distance threshold). The apertures within each group have high similarity and can share the same type of drill bit.

[0049] After clustering, there are slight differences in the designed borehole diameter within each group. A unified target borehole diameter needs to be determined as the actual construction borehole diameter and drill bit compatibility benchmark for all boreholes in that group. The determination of the target borehole diameter needs to take into account both the borehole diameter design accuracy and drill bit compatibility. In this embodiment, the arithmetic mean method is used for calculation. After the calculation is completed, the target borehole diameter is rounded (keeping the integer to match the integer specification of the drill bit model; for example, if the calculated target borehole diameter is 102.3mm, it is rounded to 102mm; if there is no 102mm drill bit on site, it can be adjusted to the closest standard drill bit diameter, such as 100mm or 105mm, within the allowable error range). At the same time, the matching of the target borehole diameter with the drill bit inventory on site is checked. If there is no corresponding drill bit for a certain group of target borehole diameters, the distance threshold of that group needs to be readjusted, and the cluster groups need to be split or merged until the target borehole diameter of all groups can match the existing drill bit models. Finally, the target borehole diameter of each group is determined to ensure that the corresponding drill bit can be obtained smoothly in subsequent construction.

[0050] S312. Based on the objective function of the number of drill bit replacements, the globally optimal scheduling scheme is obtained through multi-drilling rig collaborative scheduling optimization.

[0051] The core objective of multi-rig collaborative scheduling is to minimize the number of drill bit changes, thereby reducing construction time and drill bit wear costs. Therefore, it is necessary to construct an accurate objective function for the number of drill bit changes to quantify the drill bit change costs under different scheduling schemes. Specifically, when a single drilling rig is responsible for drilling holes of a certain diameter, the corresponding drill bit does not need to be changed initially. However, for each subsequent new drilling task with a different diameter, the drill bit needs to be changed once. Therefore, the number of drill bit changes for a single drilling rig in a certain group is the number of group-related changes corresponding to the number of holes in that group that the drilling rig is responsible for.

[0052] Furthermore, the objective function for the number of drill bit replacements satisfies: Where N represents the total number of drill bit changes for all drilling rigs (which needs to be minimized). The number of pore size clusters. This refers to the number of drilling rigs used in on-site construction. The j-th drilling rig is responsible for the number of boreholes of the i-th diameter group. The number of times the drill bit of the j-th drilling rig is replaced during the construction of the i-th borehole diameter (when... At that time, the drilling rig was not responsible for the construction of this group, and the number of replacements was -1. This number needs to be discarded during the summation process; that is, only the number of replacements is calculated. (In cases where...) By rationally allocating the borehole diameter group and the number of blast holes to each drilling rig, the total number of drill bit changes for all drilling rigs is minimized, thereby maximizing construction efficiency.

[0053] Meanwhile, to ensure the feasibility and applicability of the scheduling plan, two types of core constraints need to be set based on the actual on-site construction conditions to prevent the scheduling plan from deviating from the scope of construction capacity: ① Drill bit compatibility constraint: The target borehole diameter for each group of boreholes must match the existing drill bit model, and each drilling rig responsible for the construction of a certain group of boreholes must be equipped with the corresponding model of drill bit (which can be allocated from the on-site drill bit inventory). If a drilling rig does not have a drill bit for the corresponding group of target boreholes and cannot be allocated, then the drilling rig shall not be responsible for the construction of that group of boreholes; the number of drill bits corresponding to each group of target boreholes must meet the construction needs of that group (i.e., the number of drill bits ≥ the number of drilling rigs responsible for the construction of that group) to avoid construction stagnation due to insufficient drill bits. ② Drilling rig workload constraints: The total workload of a single drilling rig (i.e., the total number of boreholes of all groups handled by the drilling rig × the drilling depth of a single hole) shall not exceed its rated operating efficiency (unit: m / day, determined by the drilling rig technical parameter manual and on-site construction conditions, such as a down-the-hole drilling rig with a rated operating efficiency of 100m / day); the boreholes handled by a single drilling rig should be concentrated in a continuous area of ​​the blast zone to avoid frequent movement of the drilling rig within the blast zone, which would reduce construction efficiency; in addition, the number of borehole groups handled by a single drilling rig should not be too many (usually no more than 3 groups) to avoid disrupting the construction rhythm and reducing efficiency due to frequent drill bit changes.

[0054] Next, a greedy algorithm is used to solve the scheduling scheme. First, the scheduling allocation is initialized by counting the model and quantity of drill bits available for each drilling rig. Priority is given to assigning boreholes of the same diameter to rigs that already have drill bits corresponding to the target diameter, and to rigs with sufficient drill bits and high rated operating efficiency. This avoids assigning the same borehole group to rigs without corresponding drill bits (requiring additional drill bit replacement). Second, boreholes of the same diameter are allocated in batches. If multiple drilling rigs can accommodate a particular borehole group (with corresponding drill bits), the boreholes in that group are allocated to one of them (prioritizing rigs with less workload and fewer assigned groups), avoiding splitting the same borehole group among multiple rigs. The first step is to adjust and optimize the scheduling scheme. After allocation, it is verified whether the total workload of each drilling rig meets the constraints (not exceeding the rated working efficiency) and whether the number of groups it is responsible for is reasonable. If the workload of a certain drilling rig exceeds the rated value, the excess holes will be adjusted to drilling rigs with insufficient workload and suitable for the corresponding hole diameter. If a certain drilling rig is responsible for too many groups (more than 3 groups), some groups will be adjusted to drilling rigs with sufficient corresponding drill bits and fewer groups, until all constraints are met and the number of drill bit replacements is minimized, and finally the globally optimal scheduling scheme is obtained.

[0055] S313. Combining the clustering results with the global optimal scheduling scheme, the aperture grouping characteristics are obtained.

[0056] After completing the aperture clustering and multi-rig collaborative scheduling optimization, standardized feature parameters are compiled and output, including: ① Aperture grouping results: presented in a standardized table format, including the group number (G1~Gk) of each group, all borehole numbers and corresponding design apertures in the group, the target aperture of the group (rounded and matched with the drill bit model), the total number of boreholes in the group, and the spatial distribution range of the boreholes in the blasting area, clearly presenting the core information of each aperture group, facilitating identification and operation by on-site construction personnel; ② Minimum number of drill bit replacements. Output the total number of drill bit replacements (integer) after scheduling optimization, and mark the details of the number of drill bit replacements for each drilling rig (the number of replacements for the j-th drilling rig = the number of groups that the drilling rig is responsible for -1), clarifying the drill bit replacement requirements for each drilling rig, so that construction personnel can prepare drill bits in advance and plan replacement time; ③ The compatible drilling rig models and scheduling allocation results for each group: clarify the compatible drilling rig models (all compatible drilling rigs) corresponding to the borehole diameters of each group, the final assigned drilling rig number, mark the borehole diameters and corresponding number of blast holes for each drilling rig, and link the rated operating efficiency of the drilling rig and the number of existing drill bits to explain the rationality of the scheduling allocation. If there is a need for drilling rig relocation (such as drill bit transfer), add specific relocation suggestions to ensure that on-site construction can be directly executed according to the scheduling plan, improving construction efficiency and reducing construction costs.

[0057] In this embodiment, the step of calculating the comprehensive rock mass blastability index by combining the rock mass structural surface network characteristics, the rock mass wave impedance zoning characteristics, and the aperture grouping characteristics using an attention-enhanced deep neural network includes the following steps: S321. Construct an attention-enhanced deep neural network.

[0058] The architecture employs a five-level structure: input layer - hidden layer 1 - attention layer - hidden layer 2 - output layer, detailed as follows: ① Input Layer: The input dimension is set to 11, corresponding to 11 core input features, specifically 7 rock mass structural surface network features (statistical mean of trace length L, statistical mean of aperture d, standard deviation of dip α, standard deviation of dip angle β, connectivity C, linear density). areal density ), 2 wave impedance partition characteristics (mean ρv, integrity coefficient) ), 2 aperture grouping features (target aperture) Minimum number of drill bit replacements The number of neurons in the input layer is consistent with the feature dimension, which is 11. Linear activation is used, which does not change the numerical distribution of the input features, but only completes the input transmission of features. ② Hidden Layer 1: 64 neurons are set, using the ReLU activation function (rectified linear unit). The ReLU function expression is f(x)=max(0,x). Its core advantage is that it can effectively alleviate the gradient vanishing problem, accelerate model training convergence, and reduce computational complexity. The number of 64 neurons has been verified in engineering. It can fully extract the deep semantic associations of input features (such as the synergistic effect of structural surface connectivity and wave impedance) and avoid overfitting caused by too many neurons. ③ Attention layer: As the core optimization module of the network, it is embedded between hidden layer 1 and hidden layer 2. Its core function is to dynamically allocate the weights of each feature, highlight the key features that affect the blastability of the rock mass (such as the connectivity C of the structural surface and the mean wave impedance), suppress the interference of secondary features, and realize the adaptive weighted fusion of features. ④ Hidden layer 2: Set 32 ​​neurons, also using the ReLU activation function, to receive the weighted features output from the attention layer, further extract higher-order correlations between features, compress and optimize the weighted features, reduce feature dimensionality, and provide accurate feature input for the output layer. The number of 32 neurons forms a reasonable gradient with hidden layer 1 to avoid information loss during feature transmission. ⑤ Output layer: Set up 1 neuron, use sigmoid activation function, match the value range of rock mass blastability comprehensive index I (0≤I≤1), and can directly output a single blastability comprehensive index, where 0 represents rock mass that is extremely difficult to blast (intact structure, high wave impedance) and 1 represents rock mass that is extremely easy to blast (fragmented structure, low wave impedance), so as to realize the quantitative evaluation of blastability.

[0059] Furthermore, the attention layer dynamically assigns weights based on feature importance. Through learnable weight parameters, it quantifies the contribution of each feature to the explosiveness index, achieving adaptive weighted fusion of features. The specific calculation process consists of three steps: The first step is feature input and weight initialization. The 64-dimensional feature vector output from hidden layer 1 is used as the input to the attention layer, denoted as . ,in The output features of the i-th neuron in hidden layer 1 are used to initialize the learnable weight vector. The dimension of the weight vector is the same as the dimension of the output feature of hidden layer 1, and the initial values ​​adopt a random normal distribution (mean is 0, standard deviation is 0.1) to ensure the randomness and diversity of the weights; The second step is to calculate the attention weights for each input feature. Calculate the corresponding attention weights The softmax function is used to normalize the weights, ensuring that the sum of the weights of all features is 1, satisfying the following: (j ranges from 1 to 64), where For weight With features The dot product reflects the matching degree between the weights and features. The exp function is used to amplify the differences between features, and the softmax function is used to normalize the weights to make the weight distribution reasonable. The third step is to output the weighted features, which are the input features. With the corresponding attention weights Multiply to obtain the weighted features All weighted features are concatenated to form the output feature vector of the attention layer. And input it into hidden layer 2.

[0060] Furthermore, historical blasting project data is collected, and historical projects similar to the current blasting area in terms of lithology, construction conditions, and blasting scale are selected. Simultaneously, the actual explosiveness performance of each historical project after blasting (such as rock fragmentation, explosive energy utilization rate, and large block ratio) is collected. Secondly, the explosiveness performance of historical projects is quantitatively labeled to generate explosiveness tags (i.e., a comprehensive explosiveness index). The labels range from 0 to 1, and are labeled using a combination of expert scoring and actual fragmentation effects: 3-5 experts in blasting engineering are organized to give explosiveness scores (0-1) based on indicators such as rock fragmentation degree, explosive consumption per unit volume, and large block rate of historical projects. The average of the expert scores is taken as the explosiveness label for the project, where 0 represents extremely difficult to blast (intact rock mass, large block rate after blasting >20%, explosive energy utilization rate <60%), 1 represents extremely easy to blast (fragmented rock mass, large block rate after blasting <5%, explosive energy utilization rate >90%), and 0.3-0.7 represents moderate explosiveness, ensuring the accuracy and objectivity of the labels. Finally, the collected historical data is screened and divided, and invalid data with missing features or abnormal labels are removed to ensure the integrity of the dataset. The dataset is then divided into training and validation sets in an 8:2 ratio, with the training set (80%) used for model parameter fitting training and the validation set (20%) used for real-time evaluation of the model's training effect to avoid model overfitting.

[0061] During model training, the learnable weights w are continuously updated and optimized with each training iteration, gradually strengthening the weights of key features. For example, in engineering practice, the connectivity C of structural surfaces has the greatest impact on explosiveness, and its corresponding attention weights are adjusted after training. Typically the highest (generally between 0.15 and 0.2), while the number of drill bit replacements... The impact is relatively small, and the weight is usually below 0.05. Through this dynamic weight allocation, the network can focus on core features and improve the accuracy of explosiveness index calculation.

[0062] S322. By weightedly fusing the network features of the rock mass structural surface, the zoning features of the rock mass wave impedance, and the grouping features of the aperture through the attention-enhanced deep neural network, a comprehensive index of rock mass blastability is obtained.

[0063] First, all features are normalized by using a linear normalization method to map all feature parameters to the 0-1 range.

[0064] Input the characteristic parameters of the current blasting zone into the optimal model. Through forward propagation calculation, output the standardized rock mass blastability comprehensive index and related auxiliary features, including: ① Rock mass blastability comprehensive index I (0≤I≤1): a single quantitative value, retained to 3 decimal places, which intuitively reflects the strength of the blastability of the rock mass in the current blasting zone. At the same time, output the blastability level classification results to facilitate engineers to quickly judge: I<0.3 is extremely difficult to blast (the amount of explosive charge and the borehole parameters need to be increased to ensure the blasting effect), 0.3≤I≤0.7 is a medium blastability level (the requirements can be met by using conventional blasting parameters), and I>0.7 is extremely easy to blast (the amount of explosive charge needs to be reduced to avoid over-blasting and reduce costs); ② Attention weights of each feature ( The input features are 11 in total, each with a weight value (rounded to 3 decimal places). The sum of all weight values ​​is 1. The magnitude of the weight value directly reflects the contribution of the feature to the explosiveness index. For example, the structural surface connectivity C usually has the highest weight (0.15~0.20), followed by the mean wave impedance (0.12~0.18), and the number of drill bit replacements. The weight is the lowest (0.03~0.05), and the weight ranking results are output to clarify the secondary features of the core features that affect explosiveness, providing targeted guidance for subsequent blasting parameter optimization (e.g., when the explosiveness index I is low, prioritize optimizing the parameters corresponding to the core features, such as adjusting the borehole spacing and increasing the degree of structural surface fragmentation).

[0065] S4. Expand the block size distribution samples by generating adversarial networks, and construct the block size correction field of the explosion zone space by combining the simulated block size distribution characteristics and the actual block size distribution characteristics.

[0066] In this embodiment, extracting the simulated block size distribution features includes the following steps: S411. Reconstruct the time-space curve of explosive energy release by combining chemical reaction kinetics.

[0067] First, based on the Arrhenius chemical reaction kinetic model, a rate equation for explosive detonation is established: explosive detonation is essentially a violent exothermic chemical reaction, and the reaction rate is affected by factors such as temperature and the properties of the explosive itself. The Arrhenius model can accurately describe the change in the explosive reactivity over time. (Pre-exponential factor) Determined based on the type of explosive, such as emulsion explosives. The range of values ​​is s -1 Ammonium nitrate explosives The range of values ​​is s -1 The activation energy E reflects the ease or difficulty of the detonation reaction of an explosive; for emulsion explosives, E is typically [value missing]. J / mol, E of ammonium nitrate explosive is J / mol; the gas constant R is a fixed value of 8.314 J / (mol·K); the detonation temperature T (unit: K) is calculated based on the explosive heat of combustion. For emulsion explosives, T is approximately 2800~3200 K, and for ammonium nitrate explosives, T is approximately 2500~2900 K; the reaction order n reflects the relationship between the reaction rate and the concentration of the unreacted portion of the explosive. For most industrial blasting explosives, n is taken as 1.0~1.5, which can be calibrated through laboratory detonation tests.

[0068] Based on the above parameters, the detonation reaction rate equation for explosives is established as follows: Where c represents the reactivity of the explosive, ranging from 0 to 1, where c=0 indicates no reaction and c=1 indicates a complete reaction; dc / dt is the reaction rate of the explosive (unit: s⁻¹), reflecting how fast the reaction proceeds. The larger the dc / dt, the faster the explosive energy is released. By quantifying the effect of temperature on the reaction rate, the entire reaction process of the explosive, from ignition and detonation to the end of the reaction, is accurately captured, avoiding the limitations of fixed-rate models.

[0069] Furthermore, by combining the geometric model of the blast zone, the energy release rate of the explosive is coupled with the spatial position of the borehole to obtain the energy release density E(x,y,z,t) at different times and spatial positions, i.e., the energy spatiotemporal release curve: First, based on the heat of explosion per unit mass Q (unit: J / kg) of the explosive, the reaction rate is converted into the energy release rate, and the derived formula is as follows: Q is determined by actual measurement based on the type of explosive; for emulsion explosives, Q is approximately 4.0 × 10⁻⁶. 6 ~4.5×10 6 J / kg, Q of ammonium nitrate explosive is approximately 3.5 × 10 6 ~4.0×10 6 J / kg and dE / dt represent the energy released per unit mass of explosive per unit time.

[0070] Subsequently, based on the three-dimensional geometric model of the blast zone (including parameters such as borehole location, borehole diameter, borehole depth, and borehole spacing), each borehole is regarded as an independent energy release source. A spatial interpolation method is used to couple the energy release rate of a single borehole to the entire spatial domain of the blast zone: for any spatial coordinate (x, y, z) within the blast zone, the distance from that point to each borehole is calculated. The closer the distance, the greater the energy contribution of that borehole to that point. The energy release density E(x, y, z, t) at that point is obtained by weighted summation.

[0071] Finally, by continuously calculating the energy release density at various spatial points at different times (t from 0 to the time when the explosive reacts completely), an energy time-space release curve was plotted. The curve is divided into three stages: the ignition stage (t=0~t1), where the energy release rate increases slowly and c increases from 0; the detonation stage (t1~t2), where the energy release rate reaches its peak and remains stable, and c rises rapidly to close to 1; and the decay stage (t2~t3), where the energy release rate gradually decreases and c approaches 1 until the reaction ends.

[0072] S412. By using a population balance model-finite element-discrete element cross-scale coupling, the blasting process is simulated from the overall blasting area to rock particles, and the explosive energy release law and block size distribution characteristics are obtained.

[0073] The fracturing of blasted rock mass is a multi-scale evolutionary process from macroscopic (overall stress propagation in the blast zone) to microscopic (collision and slippage of rock particles) and then to block size (particle size classification). This invention adopts a cross-scale coupling method of PBM (population equilibrium model)-FEM (finite element method)-DEM (discrete element method) to progressively simulate the entire fracturing process of rock mass under the action of explosive energy, and realize the precise correlation between energy release and block size distribution.

[0074] Specifically, the macro-scale (FEM) is first based on the three-dimensional geometric model of the blast zone. Tetrahedral elements are used to construct the three-dimensional finite element model of the blast zone. The element size is set according to the scale of the blast zone and the calculation accuracy, usually 0.5~2.0m. The boundary of the blast zone adopts a non-reflective boundary condition (to avoid stress wave reflection at the boundary, which would cause simulation errors). The elements around the blast hole are densified to ensure the accuracy of the stress field simulation near the blast hole.

[0075] Then, input the rock mass mechanics parameters (compressive strength). ,tensile strength The internal friction angle φ and cohesion c were obtained through field core sampling tests and laboratory mechanical tests. Corresponding mechanical parameters were used for rock masses in different wave impedance zones to ensure spatial matching of the parameters. The reconstructed energy-space-time release curve was used as the energy input boundary condition and applied to the borehole walls. An explicit dynamic algorithm (such as the explicit solver in LS-DYNA software) was used to simulate the propagation process of the blasting stress wave: the stress wave originates from the borehole wall and spreads to the surrounding rock mass. When the tensile stress at a certain point in the rock mass exceeds its tensile strength... When the compressive stress exceeds its compressive strength, tensile fracture occurs at that point; When compression fracture occurs, the macroscopic fracture region of the rock mass (i.e., the region where the stress exceeds the rock mass strength limit) is output, clarifying the spatial range and degree of fracture of the fracture region. This helps to define the precise range for subsequent discrete element simulations at the microscale, avoiding the inefficiency caused by blindly expanding the computational domain in microscale simulations.

[0076] Microscale (DEM) simulations model the collision, slippage, and secondary fracturing processes of rock blocks within the macroscopic fracture zone, providing a data source for block-scale PBM simulations. Within the macroscopic fracture zone obtained from FEM calculations, a discrete element particle model is constructed. The particle size matches the spacing of the microscopic fractures in the rock mass, and the particle shape adopts irregular polygons (to better fit the microstructure of the rock mass). The particles are connected through a contact model (such as a linear contact model), and the contact parameters (contact stiffness, friction coefficient) are determined based on micromechanical tests of the rock mass. The friction coefficient typically corresponds to the internal friction angle φ of the rock mass.

[0077] Subsequently, the stress field data at various points within the macroscopic fracture region obtained from the FEM simulation were transformed into the initial force boundary conditions for the DEM particle model. Simultaneously, energy spatiotemporal release curves were coupled to simulate the collision, slippage, and secondary fragmentation processes between particles: the rock mass within the macroscopic fracture region was broken into initial rock blocks (particle aggregates) by stress waves. These rock blocks collided and slipped under the pressure of explosive gases and their own inertia. When the impact force generated during the collision exceeded the rock block's own strength, secondary fragmentation occurred, forming smaller rock particles. During the simulation, the motion trajectory, particle size change, and collision force of each particle were recorded in real time, providing accurate measured data support for the subsequent fragmentation kernel function of the PBM model, achieving an organic connection between the microscopic fragmentation process and block size statistics.

[0078] Population size model (PBM) is used to quantify the grain size distribution of rock fragmentation, transforming the microscopic fragmentation process simulated by DEM into quantifiable fragment size distribution characteristics, thus solving the problem that DEM simulation cannot efficiently statistically analyze the grain size distribution of large-scale rock blocks. First, the rock fragmentation process simulated by DEM is combined with a population balance model (PBM). The PBM model describes the relationship between the probability of rock fragmentation and energy through a fragmentation kernel function, and its core function is to quantify the fragmentation rate and post-fracture grain size distribution of rock blocks of different sizes.

[0079] By combining rock fragmentation data recorded in DEM simulations, the fragmentation kernel function parameters of the PBM model were calibrated, and a quantitative relationship between fragmentation rate and energy input was established. The fragmentation kernel function satisfies: ,in The number of rock blocks of size i reflects the total amount of rock blocks of that size at a certain moment; Let j be the rate at which rock particles of size j break down into i particles. The larger the value, the faster the j-sized rock fragment breaks into i-sized rock fragments. The value of is calibrated based on the fragmentation data simulated by DEM and is positively correlated with the explosive energy input density; the first term on the right side of the formula represents the impact of larger-diameter rock blocks (j>i) breaking into i-diameter rock blocks. The increment, the second term represents the result of i-sized rock fragments breaking into smaller rock fragments (j < i). The difference between the reduction and the change is the rate of change of the number of rock blocks of size i over time. By solving this differential equation, the number distribution of rock blocks of each size at different times can be obtained. Combined with the rock block density conversion, the mass distribution of rock blocks of each size can be obtained. Finally, the transformation from the microscopic crushing process to the macroscopic block size distribution is realized, providing accurate simulation data for subsequent block size prediction and correction.

[0080] Furthermore, the explosive energy time-space release curve includes two core curves: First, the energy-time curve (Et curve), with time t as the abscissa and energy release density E as the ordinate, clearly shows the change in the energy release rate of the explosive from ignition and detonation to the end of the reaction, marking key parameters such as energy peak, peak occurrence time, and reaction duration. The energy peak reflects the intensity of the explosive detonation, the peak occurrence time reflects the speed of energy release, and the reaction duration reflects the time required for the explosive to fully react. Second, the energy-space coordinate curve (E-space coordinate curve), with the spatial coordinates of the blast zone (x, y, z) as the three-dimensional coordinate axis and energy release density E as the numerical axis. It presents the differences in energy distribution at different spatial locations through color cloud maps or three-dimensional surface maps, marking high-energy areas (around the borehole, energy release density ≥ 10). 8 J / m 3 Medium energy region (central part of the explosion zone, energy release density 5×10⁻⁶) 7 ~10 8 J / m 3 ) and low-energy regions (edge ​​of the explosion zone, energy release density < 5 × 10⁻⁶) 7 J / m 3 This provides a basis for subsequent analysis of the correlation between energy distribution and rock mass fracturing.

[0081] Block size distribution characteristics are used to quantify the distribution pattern of rock block grain size, including: ① Proportion of each grain size class: The equivalent grain size of the rock blocks is divided into 4 core levels (which can be adjusted according to engineering needs), namely 0~100mm (fine blocks), 100~300mm (medium blocks), 300~500mm (large blocks), and >500mm (extra-large blocks). The mass proportion of rock blocks in each grain size class is calculated (mass is converted by rock block volume × rock mass density), and a table of the proportion of each grain size class is output to intuitively reflect the uniformity of block size distribution; ② Large block rate: The mass proportion of extra-large rock blocks with a grain size >500mm is statistically analyzed. The proportion of large blocks is a core control indicator in blasting engineering. The lower the value, the better the control over the size of the blasted fragments. This is generally a requirement in engineering projects. The simulation results should highlight this parameter; ③ Block size uniformity coefficient: The uniformity index U is used to quantify the uniformity of block size distribution, satisfying: U=1-σ / D50, where σ is the standard deviation of particle size, reflecting the dispersion of rock blocks of each particle size class. The smaller σ is, the lower the dispersion; D50 is the median particle size, which refers to the particle size corresponding to when the mass of rock blocks reaches 50%. D50 is usually used to measure the average particle size of rock blocks. The value of U ranges from 0 to 1. The closer U is to 1, the more uniform the block size distribution and the more ideal the blasting effect. At the same time, the block size distribution histogram and cumulative distribution curve are output to intuitively present the distribution characteristics of each particle size class, providing visual support for subsequent comparison with actual block size data and construction of correction fields.

[0082] In this embodiment, extracting the actual block size distribution features includes the following steps: S421. Construct a U-Net segmentation model and use the U-Net segmentation model to obtain the equivalent particle size data of each rock block.

[0083] The overall structure of the segmentation model retains the encoder-decoder architecture of U-Net. The encoder is responsible for feature extraction (from shallow texture features to deep semantic features), and the decoder is responsible for feature upsampling and contour restoration. The edge features extracted by the encoder are passed to the decoder through skip connections, which further improves the segmentation accuracy of the rock contour and effectively avoids the interference of background (such as soil, gravel, dust) on the rock segmentation, ensuring that the segmented rock contour is complete and accurate.

[0084] The channel attention mechanism calculates the weights of each feature channel, strengthens the channel weights corresponding to rock features, and weakens the channel weights corresponding to background features, thus highlighting the feature differences between the rock and the background. The spatial attention mechanism focuses on the pixel features of the rock edge region, enhances the response intensity of edge pixels, and solves the problems of blurred rock edges and unclear outlines.

[0085] The effective frame images are input into the trained U-Net segmentation model in batches. The model extracts the rock block features in the image through the encoder, strengthens the edge features and suppresses background interference by combining the CBAM attention mechanism, and then restores the pixel-level contour of the rock block through the decoder upsampling, and outputs the binary segmentation mask of each rock block (1 in the mask represents the rock block area and 0 represents the background area).

[0086] Subsequently, by converting pixels to actual size, the pixel area of ​​the rock block is converted into the actual area. The conversion process is based on the scale in the labeled image. First, the pixel length of the scale in the image is calculated to obtain the pixel-to-actual size conversion coefficient. Then, the pixel area of ​​each rock block is extracted by segmentation mask, and the actual area of ​​the rock block is calculated. Finally, the equivalent grain size of the rock block is derived according to the formula for the area of ​​a circle, and the equivalent grain size data of each rock block is obtained.

[0087] S422. Perform block size gradation statistics based on the equivalent particle size data to obtain actual block size characteristic parameters.

[0088] Based on the actual needs of blasting projects, the equivalent particle size of rock blocks is divided into four core levels: 0~100mm (fine blocks, which can be used directly for subsequent applications without secondary crushing), 100~300mm (medium blocks, which require simple crushing), 300~500mm (large blocks, which require targeted crushing), and >500mm (extra-large blocks, which require intensive crushing). If there are special requirements for the project, the particle size classification range can be flexibly adjusted.

[0089] Subsequently, the equivalent grain size of all the segmented rock blocks was classified and statistically analyzed. First, each rock block was assigned to a corresponding level according to its equivalent grain size, and the number of rock blocks in each level was counted. Then, the mass proportion of rock blocks in each level was calculated using the conversion relationship between rock block volume and rock mass density. Calculate the volume of each rock block. Then, combined with the measured density ρ of the rock mass, the mass m of each rock block is calculated; finally, the total mass of each grain size class is counted, and the mass percentage of rock blocks of that class is calculated (mass percentage of a certain class = total mass of that class / total mass of all rock blocks × 100%), and a block size distribution statistical table is output to intuitively reflect the distribution of rock blocks of different grain sizes.

[0090] Further calculations of the actual large-piece ratio and uniformity coefficient: Actual large-piece ratio It is a core control indicator for blasting engineering, with a focus on the percentage of ultra-large rock blocks with an equivalent particle size >500mm. The lower the value, the better the control over the size of the blasted fragments. This is generally a requirement in engineering projects. If the parameters exceed this range, further optimization of the blasting parameters is required.

[0091] Actual uniformity coefficient Used to quantify the uniformity of rock particle size distribution, satisfying: ,in The standard deviation of the equivalent grain size for all rock blocks reflects the degree of dispersion in rock block grain size. The smaller the value, the smaller the difference in rock grain size. The actual median grain size of the rock block refers to the equivalent grain size when the rock block accounts for 50% of the mass, reflecting the average grain size of the rock block. The value range is 0~1. The closer the value is to 1, the more uniform the rock particle size distribution, the more ideal the blasting effect, and the more effectively the cost and workload of subsequent crushing processes can be reduced.

[0092] In this embodiment, the actual block size characteristic parameters include: ① Actual block size gradation table: presenting the number of rock blocks, total mass of each block size, and mass percentage in tabular form, clearly and intuitively reflecting the actual block size distribution pattern; ② Actual large block rate : 1. Accurately quantify the proportion of extra-large rock blocks as a core error reference indicator for block size prediction correction; 2. Actual uniformity coefficient : Quantify the uniformity of the actual particle size distribution to reflect the control effect of blasting particle size; ④ Actual median particle size It reflects the average grain size of rock blocks and is used to compare the average size difference between simulated block size and actual block size, providing a quantitative basis for error correction.

[0093] The process of expanding the block size distribution samples through a generative adversarial network and constructing a block size correction field for the burst zone space by combining simulated block size distribution features and actual block size distribution features includes the following steps: S431. Expand the block degree distribution samples by generating an adversarial network to obtain a block degree distribution sample set.

[0094] Considering that the block size distribution samples of blasting are limited by engineering scenarios and are prone to problems such as small sample size and uneven sample distribution (e.g., low block size rate and scarcity of high-quality blasting samples with high uniformity), directly using them for error correction will result in insufficient accuracy of the correction field and poor generalization ability. Therefore, it is necessary to use generative adversarial networks (GANs) for data augmentation to expand the high-quality and full-coverage block size distribution sample set, provide sufficient data support for the subsequent construction of the correction field, and ensure that the correction effect is consistent with the actual engineering situation.

[0095] Specifically, a lightweight GAN model is designed based on the characteristics of block distribution (numerical, low-dimensional). The model is divided into two parts: generator G and discriminator D. The two are trained together and compete against each other to achieve the distribution of generated samples consistent with that of real samples.

[0096] The generator G employs a fully connected network structure. The input is random noise z with a dimension of 10 (the random noise follows a normal distribution N(0,1), ensuring the randomness and diversity of the noise and covering different block size distribution scenarios). It is gradually mapped through three fully connected layers (64, 32, and 3 neurons respectively), outputting simulated block size features with the same dimension as the real block size distribution features. Specifically, the output features are [simulated large block ratio]. Simulation uniformity coefficient Simulated median particle size The discriminator D also employs a fully connected network. Its input is either the true block size distribution features (real sample x) or the simulated block size features output by the generator G (generated sample G(z)). Features are extracted through two fully connected layers (32 and 16 neurons respectively), and the discriminant probability (ranging from 0 to 1) is output via a sigmoid activation function. A probability ≥ 0.5 indicates a real sample, while a probability < 0.5 indicates a generated sample, thus achieving accurate discrimination between real and fake samples. All layers of the model use the ReLU activation function (except the discriminator output layer) to alleviate the vanishing gradient problem and accelerate training convergence. Additionally, Dropout layers (with a dropout probability of 0.2) are added between layers to prevent overfitting and improve the diversity of generated samples.

[0097] Using the optimal generator G after training convergence, a large amount of random noise z (the amount of noise is set according to the sample requirements) is input to generate a sufficient number of simulated block size distribution feature samples. The focus is on supplementing small sample scenarios (such as high-quality blasting samples with low large block ratio <5% and high uniformity coefficient >0.8, and low-quality blasting samples with high large block ratio >15% and low uniformity coefficient <0.5), so as to achieve full coverage of sample distribution and solve the problem of imbalance of the original samples.

[0098] After generating samples, they need to be screened and optimized to remove abnormal samples (such as samples with a blockiness rate of <0 or >30%, a uniformity coefficient of <0 or >1, or a median particle size that exceeds the reasonable range for engineering) to ensure the rationality and practicality of the enhanced samples. Then, the screened enhanced samples are merged with the original real samples to construct a complete block size distribution sample set for subsequent error calculation and correction field construction, thereby improving the generalization ability and accuracy of the correction field.

[0099] S432. Based on the block size distribution sample set, calculate the block size error between the simulated block size distribution characteristics and the actual block size distribution characteristics. Using the two-dimensional spatial coordinates of the blasting area as the independent variable, use the Kriging interpolation algorithm to interpolate the calculated block size error to each grid position in the entire blasting area to obtain the error value corresponding to each spatial coordinate point, thus forming the blasting area spatial block size correction field.

[0100] While PBM-FEM-DEM cross-scale coupled simulation can accurately correlate explosive energy with block size distribution, factors such as rock mass heterogeneity and simplification of simulation parameters (e.g., averaging of rock mass mechanical parameters and simplification of explosive energy release) still result in some deviation between simulated and actual block sizes. Directly using simulated block sizes for subsequent parameter optimization can lead to optimization results that deviate from actual engineering requirements. Therefore, it is necessary to accurately calculate the error between simulation and reality, and construct a point-by-point correction field using spatial interpolation to correct the simulated block size position by position, thereby improving the fidelity of block size prediction.

[0101] Specifically, based on the simulated block size distribution characteristics (simulated large block rate) Simulation uniformity coefficient Simulated median particle size ) and actual block size distribution characteristics (actual large block rate) Actual uniformity coefficient Actual median particle size Based on this, the error is calculated one by one according to the corresponding characteristic parameters to ensure the pertinence and accuracy of the error calculation. The error calculation adopts the absolute difference method, and the specific calculation formula is as follows: ① Large block rate error ΔB>0 indicates that the simulated bulky ratio is higher than the actual value, and ΔB<0 indicates that the simulated bulky ratio is lower than the actual value. The larger the absolute value of ΔB, the greater the simulation deviation of the bulky ratio; ② Uniformity coefficient error ΔU>0 indicates that the simulated uniformity is higher than the actual value, and ΔU<0 indicates that the simulated uniformity is lower than the actual value, directly reflecting the simulation deviation of the particle size distribution uniformity; ③ Median particle size error ΔD50>0 indicates that the simulated median grain size is greater than the actual value, while ΔD50<0 indicates that the simulated median grain size is less than the actual value, reflecting the simulation deviation of the average grain size of the rock block.

[0102] After the calculation is completed, the error data is organized, and the spatial coordinates of the blast zone corresponding to each error value are marked to form an error dataset. The spatial distribution pattern of the error is clarified (such as larger error at the edge of the blast zone and smaller error in the middle), which provides a basis for the subsequent construction of the spatial correction field.

[0103] Considering the heterogeneity of the rock mass in the blasting area, the block size error between the simulation and the actual value is not globally constant, but varies with the spatial location of the blasting area (for example, the error is usually larger in areas with dense structural surfaces and areas with abrupt changes in wave impedance). Therefore, it is necessary to construct a spatial correction field to achieve accurate correction of the block size error at each location in the blasting area.

[0104] The correction field is constructed using the two-dimensional spatial coordinates (x,y) of the blast zone as independent variables. The Kriging interpolation algorithm is used to interpolate the calculated three major errors (ΔB, ΔU, ΔD50) to each grid position in the entire blast zone, obtaining the error value corresponding to each spatial coordinate point (x,y), and finally forming the blast zone spatial correction field Δ(x,y), which satisfies: Δ(x,y)=[ΔB(x,y),ΔU(x,y),ΔD50(x,y)], where ΔB(x,y), ΔU(x,y), and ΔD50(x,y) are the blockiness correction error, uniformity coefficient correction error, and median particle size correction error at coordinate (x,y), respectively.

[0105] During Kriging interpolation, the spatial correlation of errors is fully considered (the error values ​​of adjacent grids are strongly correlated). By calculating the spatial variability function of the error, the interpolation parameters (such as range and sill value) are determined to ensure that the interpolated error value conforms to the actual error distribution and avoids interpolation distortion. For grids with no measured error at the edge of the blast zone, boundary extension interpolation is used to ensure that the correction field covers the entire blast zone without any blank areas.

[0106] Furthermore, the simulated block size is corrected point by point to obtain a high-fidelity predicted value: After the correction field is constructed, based on the simulated block size characteristics at each location in the blast zone, the simulated block size at each spatial coordinate point (x,y) is corrected point by point using the correction logic of simulation characteristics - corresponding position error, eliminating the deviation between simulation and reality, and obtaining a high-fidelity predicted value of block size distribution. The specific correction formula is as follows: ① Corrected large block rate By subtracting the large block ratio error at the corresponding location, the predicted large block ratio is made to match the actual large block ratio distribution; ② Corrected uniformity coefficient Correcting the deviation between simulated and actual uniformity improves the accuracy of uniformity prediction; ③ Corrected median particle size This ensures that the predicted median grain size matches the actual average grain size of the rock blocks. During the correction process, the rationality of the corrected predicted values ​​must be verified simultaneously to ensure that the corrected block size predictions conform to actual engineering patterns, ultimately yielding a high-fidelity block size distribution prediction result covering the entire blasting area.

[0107] After completing GAN data augmentation, error calculation, and correction field construction, a standardized, high-fidelity block size distribution prediction correction field is output. The output features specifically include: ① Complete data of the blast zone spatial correction field: containing the three major error values ​​(ΔB(x,y), ΔU(x,y), ΔD50(x,y)) corresponding to each grid coordinate (x,y) of the blast zone, annotated with the interpolation method (Kriging interpolation) and error calculation method for easy subsequent tracking and adjustment; ② High-fidelity block size distribution prediction values: corresponding to each grid coordinate (x,y) of the blast zone, outputting the three corrected core block size features (... ), retain 3 decimal places, and clearly define the meaning and engineering significance of each feature (e.g. The corrected prediction of large block ratio is used to measure the effectiveness of blasting block size control; ( ) The corrected uniformity coefficient reflects the degree of uniformity in block size distribution; ③ A visual cloud map of the block size distribution prediction correction field intuitively presents the spatial distribution pattern of the corrected block size characteristics, which helps engineers quickly grasp the block size prediction situation at different locations in the blasting area and provides targeted guidance for subsequent blasting parameter optimization.

[0108] S5. Based on the executability constraint characteristics of the borehole grid, the comprehensive rock mass blastability index, and the spatial block size correction field of the blasting zone, solve for the Pareto solution set of the optimal blasting parameters.

[0109] The feasible area of ​​the blast zone is the basic premise for the layout of the borehole grid. By combining the terrain boundary of the blast zone and the kinematic characteristics of the explosive vehicle, the workable and non-workable areas are accurately divided to avoid the subsequent design of the borehole grid from exceeding the construction capacity and to ensure the engineering feasibility of the parameter design.

[0110] Specifically, firstly, complete boundary data of the blasting area is collected, including topographic contour coordinates, planar coordinates and dimensions of ground obstacles (such as slopes, existing structures, pipelines, densely vegetated areas, etc.). A two-dimensional planar coordinate system is constructed with the lowest point of the blasting area as the origin (the x-axis is parallel to the long side of the blasting area, and the y-axis is parallel to the short side). A grid discretization method is then used to divide the blasting area plane into uniform two-dimensional grids.

[0111] The grid marking follows the principles of hierarchical labeling and precise differentiation: areas where drilling operations are not possible, such as those within the slope toe line, within the boundary of structures, and within pipeline protection zones, are marked as infeasible grids (marked as 0); areas without obstacles and with a terrain slope ≤ the maximum climbing angle of the explosives truck are marked as infeasible grids. Areas with ground bearing capacity greater than or equal to the drilling rig's operational requirements are marked as feasible grids (marked as 1); areas with terrain slope between and The transition area between these areas is marked as a temporary feasible grid (marked as 0.5), which will be further verified in conjunction with the stability of the explosive vehicle operation. At the same time, the boundary coordinates of the blast zone, obstacle types and sizes are marked on the grid map to form a visualized basic grid map of the blast zone, providing a carrier for subsequent kinematic constraint overlay.

[0112] As the core equipment for borehole construction, the kinematic characteristics (turning radius, climbing angle, and working width) of the explosives truck directly determine the accessibility of the operation. These constraints need to be quantified and superimposed onto the basic grid map, and grids that cannot meet the equipment's operational requirements need to be eliminated to finally obtain the feasible domain for the explosives truck operation.

[0113] Specifically: taking the center of the explosives truck as the origin, based on the minimum turning radius of the explosives truck... (Obtained from the equipment manual, typically 5-8m), draw the arc range covered by the vehicle body when turning; this range is the minimum feasible area for the explosives truck to turn. For areas in the basic grid map that cannot accommodate this arc range, or whose arc range contains infeasible grids, the turning requirements are deemed unmet, and the corresponding grid is marked as an infeasible grid (updated to 0). Simultaneously, consider the explosives truck's maximum climbing angle. Temporary feasible grid (slope) Further screening: If the area has no other obstacles and the ground flatness meets the requirements (undulation difference ≤ 0.3m), it is retained as a feasible grid; if there are excessive undulations or local obstacles, it is updated to an infeasible grid. In addition, a safe distance for explosives truck operations (usually 0.8~1.0m) needs to be reserved, and the edges of feasible grids are shrunk inward by the corresponding distance to avoid collisions between the drilling rig and obstacles during operation. Finally, a complete and accurate feasible domain for explosives truck operations is obtained, and the planar coordinate range of the feasible domain, the grid distribution map, and the obstacle avoidance list are output.

[0114] Furthermore, based on the feasible domain of the explosives truck operation, the borehole grid and the equipment travel path are organically combined through graph theory navigation method. The abstract construction accessibility is transformed into quantitative borehole grid constraint parameters, ensuring that the designed borehole grid is not only within the feasible domain, but also that the travel path between adjacent boreholes meets the kinematic requirements of the explosives truck, while maximizing construction efficiency.

[0115] Specifically, the initial design of the borehole mesh is first processed into nodes, and the center coordinates of each designed borehole are taken as the vertex V of the graph, and the vertex set is... (n is the total number of designed boreholes), each vertex corresponds to the initial coordinates of a borehole; the path between any two adjacent boreholes that the explosives truck can safely travel on (path width ≥ explosives truck operating width W, path slope ≤ (Without obstacles) is transformed into the edge E of the graph, and the edge set E = {e} ij |v i With v j Connectivity is established, and the final borehole mesh graph G=(V,E) is constructed. During construction, connectivity checks are performed: if the straight-line distance between two boreholes exceeds the coverage area of ​​a single operation by the explosives truck (typically 10-15m), or if there are infeasible grid cells in the path, then the corresponding edge e is not constructed. ij If the two blast holes are not connected, temporary operation path nodes need to be added to ensure the connectivity of the entire graph G and avoid isolated blast holes (unreachable blast holes).

[0116] To optimize the explosives truck's construction path and improve operational efficiency, it is necessary to modify each edge e of graph G. ij Calculate the weight ω ij The weights are determined by the principle that shorter travel distances lead to higher construction efficiency, and thus larger weights. Dijkstra's algorithm is then used to find the shortest feasible paths from the detonation point (usually selected as a borehole in the center or edge of the blast zone for easy equipment access) to each borehole vertex: starting from the detonation point, the shortest path length from the starting point to each vertex is calculated sequentially (the path length is the sum of the travel distances of all edges on the path), and the edges and vertices of each shortest path are recorded, outputting a list of shortest paths. The core purpose of this process is to determine the optimal travel route for the explosives truck from the detonation point to each borehole, providing a basis for subsequent borehole construction sequence planning, and verifying the path reachability between adjacent boreholes.

[0117] Based on the connectivity, edge weights, and shortest path results of graph G, and combined with the kinematic parameters of the explosive vehicle, three core quantitative constraint features are extracted as the core basis for subsequent borehole mesh optimization: ① Lower limit of borehole spacing To ensure sufficient operating space for the explosives truck and prevent interference between equipment during the construction of adjacent blast holes, the following settings are provided: The operating width W of the explosives truck (W is obtained from the equipment manual and is usually 2.5~3.5m), meaning the center-to-center distance between any two adjacent blast holes must not be less than [missing value]. Otherwise, it is determined to be an infeasible blast hole; ② Upper limit of row spacing To ensure sufficient space for the explosives truck to turn between adjacent blast hole rows and to prevent equipment jamming or collisions during turns, the following settings are configured: ( The minimum turning radius of the explosives truck, i.e., the spacing between the blast holes, must not exceed [a certain value]. If the row spacing exceeds this range, the number of rows or row spacing of the boreholes needs to be adjusted; ③ Grid feasibility rate F: used to quantify the construction feasibility of the entire borehole grid. The calculation formula is F = number of feasible boreholes / total number of designed boreholes × 100%. The number of feasible boreholes refers to the number of boreholes that are within the feasible domain of the explosive vehicle operation, meet the spacing and row spacing constraints, and have an accessible path. F ≥ 85% is qualified (meets the engineering construction requirements). F < 85% requires re-optimization of the borehole grid design.

[0118] After completing the feasible region modeling and graph theory navigation constraint extraction of the blast zone, the standardized executability constraint features of the borehole mesh are obtained, including: ① Feasible region range of the borehole mesh: outputting the boundary coordinates of the feasible region in planar coordinate form. The system includes: ① Outputting a feasible grid distribution map (marking feasible and infeasible grids and obstacle locations), visually presenting the area where boreholes can be placed; ② Outputting a lower limit for borehole spacing: providing specific quantitative values, clearly indicating the source of the values ​​(explosive truck operating width W), and specifying the constraint requirements (interval between adjacent boreholes ≥ 1 / 2 ≤ ... ); ③ Upper limit of borehole spacing Output specific quantitative values, indicating the source of the values ​​(2 × minimum turning radius of the explosives truck). ), specify the constraint requirements (borehole spacing ≤ ); ④ Mesh Feasibility Rate F: Output the calculation results (retaining 2 decimal places, unit: %), annotate the number of feasible boreholes, the total number of designed boreholes and the calculation formula, and provide the qualification criteria (F≥85% is qualified). If F<85%, add feasibility rate optimization suggestions (such as adjusting borehole positions, removing isolated boreholes); ⑤ Feasible Layout Coordinates of Each Borehole: Verify each designed borehole, filter out feasible boreholes that meet all constraints, output their specific plane coordinates (x,y), annotate the borehole number, and remove infeasible boreholes, explaining the reasons for removal (such as exceeding the feasible domain, not meeting spacing constraints, no reachable path), and finally form a list of feasible borehole coordinates, providing a basic data source for subsequent borehole parameter optimization.

[0119] Furthermore, the step of solving for the optimal Pareto solution set of blasting parameters based on the executability constraint characteristics of the borehole grid, the comprehensive rock mass blastability index, and the spatial block size correction field of the blasting zone includes the following steps: S51. Combining the executability constraint characteristics of the borehole grid, the comprehensive rock mass blastability index, and the spatial block size correction field of the blasting zone, a multi-objective optimization function is constructed.

[0120] First, the six core blasting parameters that have the most significant impact on blasting block size, construction efficiency, and project cost are selected as decision variables, as follows: Hole spacing a: refers to the horizontal distance between the centers of two adjacent blast holes in the same row. It is a key parameter for controlling the distribution of blasting energy and avoiding excessive or insufficient blasting, and directly affects the uniformity of rock fragmentation. Hole row spacing b: refers to the horizontal distance between the centers of two adjacent rows of blast holes. It works in conjunction with the hole spacing to determine the superposition effect of blasting stress waves and affect the block ratio and block size uniformity. Hole diameter d: refers to the actual construction diameter of the borehole, which needs to be matched with the target hole diameter range of the hole diameter group. It directly determines the charge diameter and the charge amount per hole, and also affects the drilling efficiency. Explosive charge diameter The diameter of the explosive charge inside the borehole is usually slightly smaller than the borehole diameter (to allow space for filling with borehole clay). It is positively correlated with the borehole diameter and directly affects the energy transfer efficiency and blasting power of the explosive. Single-hole charge Q: refers to the amount of explosives loaded in a single blast hole. It is a core parameter for controlling the blasting energy and must be matched with the rock mass blastability index I to avoid insufficient energy leading to incomplete blasting or excessive energy causing waste and safety hazards. Detonation time difference Δt: refers to the time interval between detonations of adjacent rows or adjacent boreholes. Setting the detonation time difference reasonably can achieve stress wave superposition, reduce the generation of large blocks, improve block uniformity, and reduce blasting vibration.

[0121] Furthermore, a multi-objective optimization function is constructed, including: maximizing block uniformity. Minimize the large block ratio Minimize explosive consumption q and maximize construction feasibility. Construction feasibility (Value range 0~1) is a key indicator to ensure the feasibility of optimization parameters. Taking into account both the feasibility of borehole mesh and the cost of drill bit replacement, the following must be satisfied: Where F represents the grid feasibility rate (the qualification standard is F≥85%). To minimize the number of drill bit changes. This represents the total number of blast holes in the blast zone. This value quantifies the impact of drill bit replacement on construction efficiency. A higher value indicates fewer drill bit replacements and higher construction efficiency. The optimization direction is to maximize... This ensures that the optimized blasting parameters not only meet the kinematic constraints of the equipment, but also improve construction efficiency and reduce construction losses.

[0122] Furthermore, strict constraints are set to eliminate combinations of decision variables that are not consistent with engineering practice, cannot be implemented, or fail to meet blasting standards, ensuring the practicality of the optimized solution set. Specific constraints are as follows: ①Bore hole spacing and row spacing constraints: The bore hole spacing 'a' must meet the lower limit of the spacing. (Determined by the operating width of the explosives truck) and the upper limit of the row spacing (Determined by the minimum turning radius of the explosives truck), that is To ensure effective superposition of blasting stress waves and avoid blasting dead zones, the row spacing b needs to be matched with the borehole spacing a. It is usually taken as 0.6 to 0.8 times the borehole spacing, i.e., b∈[0.6a,0.8a]. If it exceeds this range, it will lead to uneven distribution of blasting energy and increase the proportion of large blocks.

[0123] ② Aperture Constraint: The aperture d must be within the target aperture range of the aperture group and must not exceed this range. This ensures that the aperture matches the existing drill bit model on site, avoiding problems such as inability to drill, waste of drill bits, or low drilling efficiency. At the same time, it ensures a reasonable match between the aperture and the charge diameter.

[0124] ③ Charge quantity and explosiveness matching constraint: The single-hole charge quantity Q must be precisely matched with the explosiveness index I. The charge quantity is adjusted according to the explosiveness of the rock mass to avoid poor blasting effect caused by energy mismatch: when I < 0.3, the rock mass is of the extremely difficult-to-blast level, and the charge quantity needs to be increased to ensure the fragmentation effect, that is... ( The minimum single-hole charge for extremely difficult-to-blast rock masses is determined by engineering experience and calculations based on rock mechanics parameters. When I > 0.7, the rock mass is classified as extremely easily blastable, and the charge amount needs to be reduced to avoid excessive blasting and energy waste. ( (This refers to the maximum single-hole charge for highly explosive rock masses, determined similarly). When 0.3 ≤ I ≤ 0.7, the rock mass is of moderate explosiveness, and Q can be... It allows for flexible adjustments within a given scope, balancing effectiveness and cost.

[0125] ④ Other auxiliary constraints: charge diameter Must meet (Reserve sufficient space for stemming to ensure blasting safety and energy transfer); the detonation time difference Δt should be within the range of 50~200ms (based on engineering experience, this range allows for effective stress wave superposition while controlling blasting vibration within a safe range); the single-hole charge Q should meet the requirement of Q≥0.5kg (to avoid insufficient charge leading to failure to detonate or insufficient blasting energy), and ( The maximum amount of explosives that a single drilling rig can load in a single operation is determined by the equipment parameters.

[0126] S52. Solve the multi-objective optimization function using a non-dominated sorting genetic algorithm with an elite strategy to obtain the Pareto front solution set, and then use TOPSIS multi-attribute decision-making to select the optimal solution from the solution set.

[0127] An initial population is constructed through random generation. For each decision variable, a value is randomly selected within its constraint range. The combination of values ​​is then verified to ensure all constraints are met. Combinations that do not meet the constraints are eliminated, and new random combinations are added until 100 valid individuals are obtained. Each individual corresponds to a set of feasible blasting parameter combinations. Simultaneously, four objective function values ​​are calculated for each individual. .

[0128] Then, non-dominated sorting is performed to rank each individual in the population and distinguish their relative merits. Non-dominated is defined as follows: if there is no other individual whose objective function values ​​are all no worse than the current individual's, and at least one objective function value is better than the current individual's, then the current individual is a non-dominated individual (i.e., no other individual can surpass it in multi-objective optimization). During the sorting process, all non-dominated individuals in the population are first selected and classified as non-dominated level 1 (optimal level); then, level 1 individuals are removed, and non-dominated individuals are selected again from the remaining individuals and classified as level 2; this process continues until all individuals in the population are classified into their corresponding levels. The smaller the level value, the better the individual's overall optimization performance.

[0129] Next, the parent population (the population in the current iteration) is merged with the offspring population (the new population generated through genetic operations). The merged population is then re-sorted for non-dominated order and the crowding degree is recalculated (crowding degree is used to measure the dispersion of individuals in the population; the higher the crowding degree, the fewer individuals around that individual, and the better the diversity of the solution set). Then, the top 100 individuals with lower rank and higher crowding degree in the merged population are retained as the next generation population to ensure that excellent individuals can be stably inherited, while maintaining the diversity of the population and avoiding the algorithm from getting trapped in local optima.

[0130] Then, by simulating the process of biological heredity and variation, a new offspring population is generated, driving the population towards a better direction of evolution. This mainly includes three steps: selection, crossover, and mutation. The non-dominated sorting, elite preservation, and genetic operations mentioned above are repeated, with a maximum of 50 iterations. When the iteration reaches 50 generations, or when the objective function value of the optimal non-dominated individuals in the population does not change significantly for 10 consecutive generations, the algorithm is considered to have converged, and the iteration stops. The set of non-dominated individuals of level 1 obtained at this time is the Pareto front solution set. This solution set contains multiple non-dominated burst parameter combinations, each of which achieves relative optimality on all four objective functions, but there is no absolutely optimal solution, which needs to be further screened using multi-attribute decision-making methods.

[0131] Since all individuals in the Pareto front solution set are non-dominated solutions, their superiority or inferiority cannot be directly determined by the objective function. Therefore, the TOPSIS (Topology-Based Solution Ranking) multi-attribute decision-making method is adopted, combined with the entropy weight method to determine the weights of each objective function (avoiding the bias of subjective weight setting). The schemes in the Pareto front solution set are comprehensively ranked, and the optimal blasting parameter scheme that best meets the actual engineering needs is selected. The specific implementation steps are as follows: Construct a decision matrix where each row corresponds to a blasting parameter scheme in the Pareto front solution set (assuming the solution set contains M schemes, i.e., M rows), and each column corresponds to one of the four objective functions. That is, to construct an M×4 decision matrix X, with matrix elements... This represents the j-th objective function value of the i-th solution (i=1,2,...,M; j=1,2,3,4, corresponding to 4 objectives respectively). During the construction process, it is necessary to ensure that all objective function values ​​have consistent units and accurate values, all derived from the calculation results during the NSGA-II solution process, while eliminating any abnormal solutions that may exist in the solution set (such as an objective function value that exceeds the reasonable range of engineering).

[0132] Since the four objective functions have different dimensions and ranges, direct comparison will introduce bias. Therefore, the decision matrix needs to be normalized to map all objective function values ​​to the 0-1 interval, eliminating dimensional differences and making the objectives comparable. After normalization, a standardized decision matrix is ​​obtained. The matrix elements are all in the range of 0 to 1.

[0133] Next, the weights are calculated. The entropy weight method is used to determine the weights of the four objective functions. The specific calculation steps are as follows: Calculate the entropy value of the j-th objective function. : Where k = 1 / ln(M), if ,Pick (To avoid the logarithm being meaningless), entropy value The value range is 0~1. The smaller the value, the greater the dispersion of the objective function, the more information it contains, and the greater its impact on the ranking of the schemes.

[0134] Calculate the difference coefficient of the j-th objective function : The larger the difference coefficient, the stronger the discriminative power of the objective function, and the higher the weight should be.

[0135] Calculate the weights of the j-th objective function : The sum of all weights This ensures a reasonable allocation of weights. Based on practical experience in blasting projects, the typically obtained weight range is: Weight 0.35~0.4 (core indicator of explosive effect). 0.25~0.3 (key indicator of blasting effect), q weight 0.15~0.2 (cost indicator). Weighting ranges from 0.1 to 0.15 (construction feasibility index), with typical weighting allocations as follows: This aligns with the project's priority requirements.

[0136] Next, the proximity score is calculated. The proximity score measures how close each solution is to the ideal solution. The higher the proximity score, the closer the solution is to the optimal state. The specific calculation steps are as follows: Determine the positive and negative ideal solutions: The positive ideal solution A+ is the set of optimal values ​​for each objective function, i.e. The ideal solution corresponds to the one that best performs all objective functions; the negative ideal solution A- is the set of worst values ​​for each objective function, i.e. This corresponds to the solution with the worst performance across all objective functions.

[0137] Calculate the Euclidean distance from each solution to the positive and negative ideal solutions: Euclidean distance is used to quantify the spatial distance between a solution and the ideal solution; the closer the distance, the better the solution.

[0138] Next, calculate the proximity. ,satisfy: , The value range is 0~1. The larger the value, the farther the solution is from the negative ideal solution and the closer it is to the positive ideal solution, indicating a better overall performance. When the solution is exactly the same as the ideal solution (theoretically optimal); when When the solution is exactly the same as the negative ideal solution (worst case), the solution is completely consistent with the negative ideal solution.

[0139] Based on the similarity of all options Sort in descending order. The most efficient solution is the one with the best overall blasting parameters; meanwhile, the top 10 solutions with the highest approximation are retained as alternatives to address the needs of different engineering scenarios (e.g., in some scenarios where cost control is prioritized, solutions with smaller q parameters can be selected). The top-ranked solutions are selected. After sorting, it is verified whether the blasting parameters of the optimal and alternative solutions meet all constraints. If any solutions do not meet the constraints, they are removed and the solutions are re-sorted to ensure that all selected solutions are feasible for engineering purposes.

[0140] The final results are divided into two categories, as follows: ① The Pareto front solution set contains detailed information on all schemes in the solution set, specifically: scheme number, 6 core blasting parameters (borehole spacing a, row spacing b, borehole diameter d, charge diameter). The four objective function values ​​are: single-hole charge amount Q, detonation time difference Δt, and total charge amount Q. The non-dominance level and congestion level clearly present the optimization effect and characteristics of each solution, making it easier for engineers to make flexible choices based on actual scenarios.

[0141] ②TOPSIS screening results: Includes the overall optimal blasting parameter scheme (labeled with scheme number, all blasting parameters and objective function values, and proximity). The ranking of alternative solutions, the weight allocation results (entropy, difference coefficient, and weight of the four objective functions), and the positive and negative ideal solutions.

[0142] Please see Figure 2 In an embodiment, to efficiently execute the blasting parameter optimization method for blasting block size control provided by the present invention, the present invention also provides a blasting parameter optimization system for blasting block size control, comprising: an input device 1, an output device 2, a processor 3, and a memory 4, wherein the input device 1, output device 2, processor 3, and memory 4 are interconnected, and the memory 4 stores program instructions for executing the steps of the blasting parameter optimization method for blasting block size control. The blasting parameter optimization system for blasting block size control of the present invention has a compact structure and stable performance, and can stably execute the blasting parameter optimization method for blasting block size control of the present invention, further improving the overall applicability and practical application capability of the present invention.

[0143] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention, and they should all be covered within the scope of the present invention.

Claims

1. A method for optimizing blasting parameters for controlling the size of blasted blocks, characterized in that, Includes the following steps: The network features of rock mass structural surfaces are extracted through digital identification and topological analysis of structural surfaces; The wave impedance zoning characteristics of the rock mass were analyzed by using measurement-while-drilling inversion and wave impedance clustering. Based on the borehole diameter clustering, borehole diameter grouping features are obtained. Combining the rock mass structural surface network features, the rock mass wave impedance zoning features, and the borehole diameter grouping features, the comprehensive rock mass blastability index is calculated using an attention-enhanced deep neural network. By generating adversarial networks to expand block size distribution samples, and combining simulated block size distribution characteristics with actual block size distribution characteristics, a block size correction field for the explosion zone space is constructed. Based on the executability constraint characteristics of the borehole grid, the comprehensive rock mass blastability index, and the spatial block size correction field of the blasting zone, the Pareto solution set of the optimal blasting parameters is obtained.

2. The method for optimizing blasting parameters for controlling blasting block size according to claim 1, characterized in that, The extraction of rock mass structural surface network features through structural surface digital identification and topological analysis includes the following steps: Based on the point cloud dataset of rock mass surface, the geometric parameters of structural surfaces are obtained by segmenting the point cloud using the region growing method. Based on the geometric parameters of the structural surfaces, a topological analysis network of the rock mass structural surfaces is constructed to obtain the network characteristics of the rock mass structural surfaces.

3. The method for optimizing blasting parameters for controlling blasting block size according to claim 1, characterized in that, The method of analyzing the wave impedance zoning characteristics of rock mass using measurement-while-drilling inversion and wave impedance clustering includes the following steps: A wave impedance inversion model was constructed using segmented drilling parameter datasets to obtain wave impedance profile data for all boreholes. The wave impedance profile data is subjected to wave impedance clustering and partitioning to generate a wave impedance partition map of the blast area, thereby obtaining the wave impedance partition characteristics of the rock mass.

4. The method for optimizing blasting parameters for controlling blasting block size according to claim 1, characterized in that, The method of obtaining borehole diameter grouping features based on borehole diameter clustering includes the following steps: Clustering of designed apertures using hierarchical clustering algorithms; Based on the objective function of the number of drill bit replacements, the globally optimal scheduling scheme is obtained through multi-drilling rig collaborative scheduling optimization. By combining the clustering results and the global optimal scheduling scheme, aperture grouping characteristics are obtained.

5. The method for optimizing blasting parameters for controlling blasting block size according to claim 1, characterized in that, The method of combining the rock mass structural surface network characteristics, the rock mass wave impedance zoning characteristics, and the aperture grouping characteristics to calculate the comprehensive rock mass blastability index using an attention-enhanced deep neural network includes the following steps: Construct attention-enhanced deep neural networks; The comprehensive index of rock mass blastability is obtained by weightedly fusing the network features of the rock mass structural surface, the rock mass wave impedance zoning features, and the aperture grouping features through the attention-enhanced deep neural network.

6. The method for optimizing blasting parameters for controlling blasting block size according to claim 1, characterized in that, The process of expanding the block size distribution samples through a generative adversarial network and constructing a block size correction field for the burst zone space by combining simulated block size distribution features and actual block size distribution features includes the following steps: By generating adversarial networks to expand the block degree distribution samples, a block degree distribution sample set is obtained; Based on the block size distribution sample set, the block size error between the simulated block size distribution characteristics and the actual block size distribution characteristics is calculated. Using the two-dimensional spatial coordinates of the blasting area as the independent variable, the Kriging interpolation algorithm is used to interpolate the calculated block size error to each grid position in the entire blasting area, thereby obtaining the error value corresponding to each spatial coordinate point and forming the blasting area spatial block size correction field.

7. The method for optimizing blasting parameters for controlling blasting block size according to claim 6, characterized in that, Extracting the simulated block size distribution features includes the following steps: Reconstructing the time-space curve of explosive energy release by combining chemical reaction kinetics; By using a cross-scale coupling of population balance model, finite element method, and discrete element method, the blasting process is simulated from the overall blasting zone to rock particles, and the energy release law and block size distribution characteristics of explosives are obtained.

8. The method for optimizing blasting parameters for controlling blasting block size according to claim 6, characterized in that, Extracting the actual block size distribution features includes the following steps: A U-Net segmentation model is constructed, and the equivalent grain size data of each rock block is obtained using the U-Net segmentation model; Based on the equivalent particle size data, block size gradation statistics are performed to obtain the actual block size characteristic parameters.

9. The method for optimizing blasting parameters for controlling blasting block size according to claim 1, characterized in that, The process of solving for the optimal Pareto solution set of blasting parameters based on the executability constraint characteristics of the borehole grid, the comprehensive rock mass blastability index, and the spatial block size correction field of the blasting zone includes the following steps: A multi-objective optimization function is constructed by combining the executability constraint characteristics of the borehole grid, the comprehensive index of rock mass blastability, and the spatial block size correction field of the blasting zone. The multi-objective optimization function is solved by a non-dominated sorting genetic algorithm with an elitist strategy to obtain the Pareto front solution set, and then the optimal solution is selected from the solution set by TOPSIS multi-attribute decision-making.

10. A blasting parameter optimization system for controlling blasting block size, characterized in that, The blasting parameter optimization system for controlling blasting block size includes: an input device, an output device, a processor, and a memory, wherein the input device, output device, processor, and memory are interconnected, and the memory stores program instructions for executing the blasting parameter optimization method for controlling blasting block size as described in any one of claims 1-9.

Images