Blasting zoning explosion control method for long and large tunnel in environment sensitive area

Abstract

The invention provides a zoning explosion control method for blasting of a long and large tunnel in an environment-sensitive area, and relates to the technical field of blasting of the long and large tunnel in the environment-sensitive area, and the method comprises the following steps: collecting earth surface data, tunnel data and auxiliary data, preprocessing the data, performing rock mass intelligent zoning based on the preprocessed data, and aiming at each zone, performing rock mass intelligent zoning; and carrying out blasting parameter optimization to obtain blasting parameters, and carrying out blasting subarea blasting control based on the blasting parameters of each area. By means of the method, zoned explosion control of the explosion area can be achieved, the explosion control accuracy is improved, and accidents are prevented.

Description

Technical Field

[0001] The present invention relates to the technical field of long tunnel blasting in environmentally sensitive areas, and particularly to a method for controlled blasting in zones for long tunnel blasting in environmentally sensitive areas. Background Art

[0002] In traditional long tunnel blasting projects, the identification of geological conditions highly relies on manual experience and has systematic defects:

[0003] 1. Low efficiency of manual survey: Geological engineers collect data on site through tools such as compasses and hammer blows. It takes 2 - 3 hours to complete the basic judgment for a single heading face. According to statistics, the error rate of manually judging the hardness of rock masses is as high as 25% - 40%. Especially in areas with concealed structures, the proportion of cases where misjudgment leads to a deviation in explosive charge of more than 30% reaches 17%.

[0004] 2. Lack of quantitative analysis means: The existing technology lacks the ability to digitally model fracture networks and joint densities. For example, in a certain railway tunnel project, due to the failure to identify a concealed fracture group, the amount of collapse after blasting reached 520m 3 , resulting in a direct economic loss of more than 8 million yuan.

[0005] 3. Serious fragmentation of data: Surface topography, tunnel internal structure, and historical monitoring data are scattered in paper records and isolated systems, lacking multi-source data fusion analysis. In a certain cross-river tunnel project, due to the failure to correlate surface water systems and tunnel burial depth data, blasting caused a sudden drop in the groundwater level, resulting in a surface settlement accident.

[0006] Therefore, it is very necessary to design a method for controlled blasting in zones for long tunnel blasting in environmentally sensitive areas. Summary of the Invention

[0007] In order to overcome the deficiencies of the prior art, the purpose of the present invention is to provide a method for controlled blasting in zones for long tunnel blasting in environmentally sensitive areas.

[0008] To achieve the above purpose, the present invention provides the following solutions:

[0009] The present invention provides a method for controlled blasting in zones for long tunnel blasting in environmentally sensitive areas, including:

[0010] Step 1: Collect surface data, tunnel data, and auxiliary data, and perform preprocessing on them;

[0011] Step 2: Perform intelligent zoning of rock masses based on the preprocessed data;

[0012] Step 3: Optimize blasting parameters for each zone to obtain blasting parameters;

[0013] Step 4: Perform controlled blasting in zones based on the blasting parameters of each zone.

[0014] Preferably, in step 1, surface data, tunnel data and auxiliary data are collected, specifically:

[0015] Collect surface data based on UAV lidar and multispectral cameras. Among them, the surface data includes terrain elevation model, surface lithology spectral characteristics, surface fracture distribution density and GPS coordinates of surface sensitive targets;

[0016] Collect tunnel data based on 3D laser scanners, ground-penetrating radars and digital borehole cameras. Among them, the tunnel data includes face rock mass hardness, joint density, rock layer dip / strike, groundwater permeability coefficient and 3D fracture network model;

[0017] Obtain auxiliary data based on the historical database. Among them, the auxiliary data includes mechanical parameters of geological borehole cores, historical blasting vibration monitoring records, regional seismic wave propagation characteristics and explosive performance parameter library.

[0018] Preferably, in step 1, the surface data, tunnel data and auxiliary data are preprocessed, specifically:

[0019] Obtain the surface curvature C based on the terrain elevation model through slope / curvature calculation and terrain roughness index calculation surf ;

[0020] Obtain the lithology code based on the surface lithology spectral characteristics through multispectral band fusion and support vector machine lithology classification;

[0021] Obtain the surface fracture density F based on the surface fracture distribution density through UAV image fracture extraction and density kernel estimation calculation surf ;

[0022] Calculate the distance D from the sensitive area based on the GPS coordinates of the surface sensitive targets through spatial topological relationship calculation and buffer analysis;

[0023] Obtain the standard deviation σ of the fracture spacing based on the 3D fracture network model through point cloud skeleton extraction and fracture spacing statistical analysis.

[0024] Preferably, in step 2, intelligent rock mass zoning is carried out based on the preprocessed data, specifically:

[0025] Step 201: Standardize the face rock mass hardness f, groundwater permeability coefficient K and the distance D from the sensitive area;

[0026] Step 202: Calculate the rock mass integrity index RII, environmental sensitivity coefficient ES and fracture development index FDI;

[0027] Step 203: Based on the Rock Mass Integrity Index (RII), Environmental Sensitivity Coefficient (ES), Fracture Development Index (FDI), the longitudinal wave velocity Vp in the regional seismic wave propagation characteristics, and the historical PPV mean value in the historical blasting vibration monitoring records, automatically divide the intelligent rock mass zones through a spectral clustering algorithm improved based on the phase velocity.

[0028] Preferably, in step 201, standardize the hardness f of the tunnel face rock mass, the groundwater permeability coefficient K, and the distance D from the sensitive area, specifically as follows:

[0029] Standardize the hardness f of the tunnel face rock mass based on the Min-Max normalization method;

[0030] Standardize the groundwater permeability coefficient K based on the logarithmic transformation normalization method;

[0031] Standardize the distance D from the sensitive area based on the piecewise function mapping.

[0032] Preferably, in step 202, calculate the Rock Mass Integrity Index (RII), Environmental Sensitivity Coefficient (ES), and Fracture Development Index (FDI), specifically as follows:

[0033] The Rock Mass Integrity Index (RII) is calculated through the hardness f of the tunnel face rock mass and the joint density J d and is:

[0034] RII = 0.6tanh(0.1f) + 0.4exp(-0.02J d );

[0035] The Environmental Sensitivity Coefficient (ES) is calculated through the distance D from the sensitive area and the surface curvature C surf and is:

[0036] ES = 0.7×(1 / D) + 0.3×C surf ;

[0037] The Fracture Development Index (FDI) is calculated through the standard deviation σ of the fracture spacing and the rock layer dip angle θ and is:

[0038] FDI = σ / (1 + θ / 90).

[0039] Preferably, in step 203, based on the Rock Mass Integrity Index (RII), Environmental Sensitivity Coefficient (ES), Fracture Development Index (FDI), the longitudinal wave velocity Vp in the regional seismic wave propagation characteristics, and the historical PPV mean value in the historical blasting vibration monitoring records, automatically divide the intelligent rock mass zones through a spectral clustering algorithm improved based on the phase velocity, specifically as follows:

[0040] The input vector is determined as the rock mass integrity index RII, the environmental sensitivity coefficient ES, the fracture development index FDI, the longitudinal wave velocity Vp in the regional seismic wave propagation characteristics, and the historical PPV mean value in the historical blasting vibration monitoring records, that is:

[0041] V = [RII, ES, FDI, V p , historical PPV];

[0042] Based on the input vector V, automatic division of intelligent rock mass zones is carried out through an accelerated decomposition spectral clustering algorithm based on cosine similarity.

[0043] Preferably, in step 3, for each zone, blasting parameter optimization is carried out to obtain blasting parameters, specifically:

[0044] Construct a data set;

[0045] Construct a blasting parameter optimization model based on the Attention-Stacked LSTM network structure;

[0046] Train the blasting parameter optimization model based on the data set;

[0047] Input the input parameters of each zone into the trained blasting parameter optimization model to obtain optimized blasting parameters.

[0048] According to the specific embodiments provided by the present invention, the following technical effects are disclosed by the present invention:

[0049] The present invention provides a method for blasting zone control blasting of long tunnels in environmentally sensitive areas. The method includes: collecting surface data, tunnel data and auxiliary data, and preprocessing them. Based on the preprocessed data, intelligent rock mass zoning is carried out. For each zone, blasting parameter optimization is carried out to obtain blasting parameters. Based on the blasting parameters of each zone, blasting zone control blasting is carried out. The present invention can achieve blasting zone control blasting in the blasting area, improve the accuracy of blasting control, and prevent accidents. Description of the Drawings

[0050] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for use in the embodiments. Obviously, the drawings in the following description are only some embodiments of the present invention. For those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative efforts.

[0051] Figure 1 It is a flowchart of the method provided by the embodiment of the present invention;

[0052] Figure 2Schematic diagram of the Encoder-Decoder model structure with the Attention mechanism added;

[0053] Figure 3 Schematic diagram of the Stacked LSTM model structure;

[0054] Figure 4 Schematic diagram of the Attention-Stacked LSTM network structure. Specific implementation manners

[0055] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0056] The purpose of the present invention is to provide a method for controlled blasting in sections of long tunnels in environmentally sensitive areas, which can achieve controlled blasting in sections of the blasting area, improve the accuracy of controlled blasting, and prevent accidents.

[0057] To make the above objects, features, and advantages of the present invention more obvious and understandable, the present invention will be further described in detail below in conjunction with the accompanying drawings and specific implementation manners.

[0058] Figure 1 For the method flow chart provided in the embodiment of the present invention, as Figure 1 shown, the present invention provides a method for controlled blasting in sections of long tunnels in environmentally sensitive areas, including:

[0059] Step 1: Collect surface data, tunnel data, and auxiliary data, and preprocess them;

[0060] Step 2: Perform intelligent zoning of rock masses based on the preprocessed data;

[0061] Step 3: Optimize blasting parameters for each zone to obtain blasting parameters;

[0062] Step 4: Perform controlled blasting in sections based on the blasting parameters of each zone.

[0063] In Step 1, collecting surface data, tunnel data, and auxiliary data specifically includes:

[0064] Collect surface data based on unmanned aerial vehicle lidar and multispectral cameras. Among them, the surface data includes terrain elevation models, surface lithological spectral characteristics, surface fracture distribution densities, and GPS coordinates of surface sensitive targets;

[0065] Tunnel data is collected based on 3D laser scanners, ground-penetrating radars, and digital borehole cameras. Among them, the tunnel data includes the hardness of the face rock mass, joint density, rock layer dip / strike, groundwater permeability coefficient, and 3D fracture network model;

[0066] Auxiliary data is obtained based on the historical database. Among them, the auxiliary data includes the mechanical parameters of geological borehole cores, historical blasting vibration monitoring records, regional seismic wave propagation characteristics, and explosive performance parameter libraries;

[0067] It should be noted that the collected data can be processed according to specific requirements to improve its accuracy. Some embodiments are given in the present invention. For example:

[0068] For the point cloud data in the surface data and tunnel data, the processing objectives are spatial registration and feature enhancement. The specific process includes:

[0069] 1. Noise filtering is performed based on the improved RANSAC algorithm (the dynamic distance threshold δ can be set to 0.05 m);

[0070] 2. Coordinate unification is performed based on the multi-scale ICP algorithm combined with the FPFH feature descriptor;

[0071] For the image data in the data, the Brown-Conrady model can be used to correct the lens distortion;

[0072] For the data-type data in the data, the following operations can be taken to process it:

[0073] 1. Missing value processing is performed based on KNN imputation;

[0074] 2. Outlier detection is performed based on the isolation forest algorithm;

[0075] For the time-series monitoring data in the data, the following operations can be taken to process it:

[0076] 1. Denoising processing is performed based on the improved wavelet threshold denoising. Among them, its adaptive threshold formula is:

[0077] λ = σ × sqrt(2 × log(N)) / log(j + 1);

[0078] In the formula, σ is the noise standard deviation, and j is the wavelet decomposition level.

[0079] In step 1, the surface data, tunnel data, and auxiliary data are preprocessed. Specifically:

[0080] The surface curvature C is obtained based on the terrain elevation model through slope / curvature calculation and terrain roughness index calculation surf ;

[0081] Based on the spectral characteristics of surface lithology, through multispectral band fusion and lithology classification by support vector machine, a lithology code is obtained;

[0082] Based on the distribution density of surface fissures, the surface fissure density F is calculated through the extraction of fissures from UAV images and density kernel estimation; surf ;

[0083] Based on the GPS coordinates of surface sensitive targets, the distance D from the sensitive area is calculated through spatial topological relationship calculation and buffer analysis;

[0084] Based on the three-dimensional fissure network model, the standard deviation σ of the fissure spacing is obtained through point cloud skeleton extraction and fissure spacing statistical analysis.

[0085] In step 2, based on the preprocessed data, intelligent zoning of rock masses is carried out, specifically:

[0086] Step 201: Standardize the hardness f of the rock mass at the tunnel face, the groundwater permeability coefficient K, and the distance D from the sensitive area;

[0087] Step 202: Calculate the rock mass integrity index RII, the environmental sensitivity coefficient ES, and the fissure development index FDI;

[0088] Step 203: Based on the rock mass integrity index RII, the environmental sensitivity coefficient ES, the fissure development index FDI, the longitudinal wave velocity V in the regional seismic wave propagation characteristics p and the historical PPV mean value in the historical blasting vibration monitoring records, automatic division of the intelligent zoning of rock masses is carried out through a spectral clustering algorithm improved based on the phase velocity.

[0089] In step 201, standardizing the hardness f of the rock mass at the tunnel face, the groundwater permeability coefficient K, and the distance D from the sensitive area, specifically:

[0090] Standardize the hardness f of the rock mass at the tunnel face based on the Min-Max normalization method;

[0091] Standardize the groundwater permeability coefficient K based on the logarithmic transformation normalization method;

[0092] Standardize the distance D from the sensitive area based on piecewise function mapping.

[0093] In step 202, calculate the rock mass integrity index RII, the environmental sensitivity coefficient ES, and the fissure development index FDI, specifically:

[0094] The rock mass integrity index RII is calculated through the hardness f of the rock mass at the tunnel face and the joint density J d and is:

[0095] RII = 0.6tanh(0.1f) + 0.4exp(-0.02J d );

[0096] The environmental sensitivity coefficient ES is calculated through the distance D from the sensitive area and the surface curvature C surf and is:

[0097] ES = 0.7×(1 / D) + 0.3×C surf ;

[0098] The fracture development index FDI is calculated through the standard deviation σ of the fracture spacing and the rock stratum dip angle θ and is:

[0099] FDI = σ / (1 + θ / 90).

[0100] In step 203, based on the rock mass integrity index RII, the environmental sensitivity coefficient ES, the fracture development index FDI, the longitudinal wave velocity Vp in the regional seismic wave propagation characteristics, and the historical PPV mean value in the historical blasting vibration monitoring records, an automatic division of the rock mass intelligent zone is carried out through a spectral clustering algorithm improved based on the phase velocity, specifically:

[0101] Determine that the input vector is the rock mass integrity index RII, the environmental sensitivity coefficient ES, the fracture development index FDI, the longitudinal wave velocity V p in the regional seismic wave propagation characteristics, and the historical PPV mean value in the historical blasting vibration monitoring records, that is:

[0102] V = [RII, ES, FDI, V p , historical PPV];

[0103] Based on the input vector V, an automatic division of the rock mass intelligent zone is carried out through an accelerated decomposition spectral clustering algorithm based on the cosine similarity;

[0104] Next, a detailed description of the accelerated decomposition spectral clustering algorithm based on the cosine similarity is given:

[0105] Traditional spectral clustering algorithms have defects in aspects such as the construction of the similarity matrix, the operation complexity of the power method algorithm, and the instability of subspace clustering. Therefore, the present invention proposes an accelerated decomposition spectral clustering algorithm (CLSC) based on the cosine similarity. Compared with traditional spectral clustering algorithms, the present invention uses the cosine distance to replace the original Euclidean distance when constructing the similarity matrix, and constructs the similarity matrix only relying on the self-information of the data points, avoiding the problems brought by manual parameter adjustment. At the same time, the present invention constructs the feature vector of the data points jointly by the cosine similarity and the coordinate position. The present invention introduces the Lanczos iterative algorithm, which can accelerate the solution of the eigen-space of large sparse matrices and obtain the clustering results faster and more accurately.

[0106] First, introduce the improvement idea:

[0107] The accelerated decomposition spectral clustering algorithm based on cosine similarity mainly proposes improvements in the similarity matrix and the calculation of the eigenvectors and eigenvalues of the Laplacian matrix. Most traditional spectral clustering algorithms use the Euclidean distance to measure the similarity between samples. The scale parameter σ needs to be set manually. An inappropriate selection affects the construction of the similarity matrix and, in severe cases, affects the accuracy of the entire algorithm. Using the cosine distance to measure similarity avoids the parameter setting step. When performing spectral decomposition, the Lanczos iterative algorithm is used to approximately estimate the similarity matrix and its leading eigenvectors, reducing its time complexity and effectively improving the running speed. Moreover, the reduction in space complexity enables it to handle large-scale data better and more stably.

[0108] Next, introduce the algorithm model:

[0109] 1. Construction of the cosine similarity matrix:

[0110] The sample set is decomposed into data points. The similarity matrix is a matrix composed of the similarities between each data point in the sample, and its form is as follows:

[0111]

[0112] Among them, each element w ij represents the similarity between the i-th sample and the j-th sample. The similarity between samples can be calculated by the Euclidean distance formula W ij = d ij Let the sample be represented by a set X = (x1, x2,... x n ) with n data points. Then the Euclidean distance between two data points is:

[0113]

[0114] Expressed in terms of cosine similarity:

[0115]

[0116] 2. Iterative calculation of the Lanczos algorithm

[0117] The Lanczos algorithm is an iterative algorithm for approximately solving the eigenvalue problem of symmetric matrices. This algorithm is based on a Krylov subspace, which consists of an initial vector and several matrix product vectors. In each iteration step, the algorithm calculates a new product vector and adds it to the Krylov subspace. Then, the Gram - Schmidt orthogonalization process is used to ensure the linear independence between the basis vectors. In each iteration step, the algorithm calculates a symmetric tridiagonal matrix, which is generated by the basis vectors and the matrix product vectors within the Krylov subspace. By calculating the eigenvalues and eigenvectors of this tridiagonal matrix, the eigenvalues and eigenvectors of the original matrix can be obtained.

[0118] There are various algorithms for the eigenvalue problem of real symmetric matrices. Most solutions first diagonalize matrix A. Due to the symmetry of A, in the first step, only an orthogonal matrix Q needs to be found such that:

[0119] T = Q T AQ;

[0120] To obtain the eigenpairs of matrix A, it can first be transformed into a tridiagonal matrix T, and then the eigenvalues and corresponding eigenvectors of T are solved. Compared with directly solving the eigenvalue problem of A, this method is simpler. The commonly used implementation is to use the Householder transformation to transform matrix A into a tridiagonal matrix T. In the process of transforming matrix A into a tridiagonal matrix, common algorithms include QR iteration, Rayleigh quotient iteration, divide - and - conquer method, bisection method, inverse iteration algorithm, and Jacobi method, etc. However, when dealing with large - scale sparse matrix A, these algorithms will generate a large number of intermediate matrices, resulting in the loss of sparsity, occupying a large amount of memory, and bringing difficulties to calculating eigenvalues. Although for some special types of sparse matrix A, using Jacobi transformation for tridiagonalization can reduce the loss of sparsity to a certain extent, through continuous research, many improved algorithms have been proposed, which have improved in terms of convergence and accuracy. However, in most cases, it is difficult for any step - by - step transformation method to maintain the sparsity of the intermediate matrices during the similarity transformation process.

[0121] To make full use of sparsity and minimize the memory occupation as much as possible, a method of directly calculating matrix T and elements is needed to achieve tridiagonalization. Based on this idea, Lanczos first proposed a method in 1950 to gradually transform a general matrix into a tridiagonal matrix, namely the famous Lanczos iteration.

[0122] 3. The Accelerated Approximation Strategy of the Lanczos Algorithm

[0123] Specifically: Let the similarity matrix W ∈ R n×n, decomposing the matrices Q and T in the above formula gives:

[0124] Q = [q1, q2, …, q n ;

[0125]

[0126] Since the matrix Q is an orthogonal matrix, so Q -1 = Q T , we get:

[0127] WQ = QT;

[0128] Rewriting it in vector form gives:

[0129] Wq i = β i-1 q i-1 + a i q i + β i q i+1 , i = 1, 2, …, n;

[0130] Let β0q0 = β n q n = 0, according to the orthogonality of the vectors q i , it is deduced that:

[0131]

[0132] On the other side, let any given q I ∈ R n and α i q I = 1, from the above three formulas, the matrices Q and T are deduced. In summary, performing iteration gives the well-known Lanczos iteration, which is the main body of the Lanczos algorithm. The vector q i in it is called the Lanczos vector. From the vector formula, the matrix formula is deduced as:

[0133]

[0134] At this time, the matrix T j is called the j-order Lanczos matrix, and Q j is the Lanczos vector matrix. If the β i ≠ 0 (i = 1, 2, …, n - 1) always holds in the Lanczos iteration process, then the generated matrices Q n and T n will satisfy the equation:

[0135]

[0136] After the matrix \(w\) is tridiagonalized, when \(\beta\) i0 = 0 in a certain step of the Lanczos iteration, the iteration process stops. This helps to solve the eigenvalue problem. Based on this, it can be deduced that before the \(i_0\)-th step, an invariant subspace of the matrix \(w\) appears, and thus the first \(i_0 - 1\) eigenvalues of \(W\) can be obtained. It should be noted that the matrix \(w\) remains unchanged throughout the Lanczos iteration, and only involves the product of \(w\) and vectors. Therefore, the Lanczos iteration can make full use of the sparsity of \(w\) and is applicable to the tridiagonalization of sparse matrices. Assume that the orthogonal matrix \(Q=[Q\) k ,Q\) u , where \(Q\) k is an \(n\times k\) matrix, and \(Q\) u is an \(n\times(n - k)\) matrix. In fact, after \(k\) steps of Lanczos iteration, \(Q\) k has been calculated, while \(Q\) u has not been determined.

[0137]

[0138] Among them, \(T\) k is a \(k\times k\) matrix, and \(T\) u is an \((n - k)\times(n - k)\) matrix. The Rayleigh-Ritz method process uses eigenvalues \(\lambda_1\leq\lambda_2\leq.......\leq\lambda\) k , which are also approximate eigenvalues of \(W\). These approximate values are called Ritz values. Let \(T\) k = \(VAV\) T be the eigenvalue decomposition of \(T\) k . The corresponding eigenvectors are the columns \(y\) k of \(Q\) i \(V\), \(y\) k \(v\) i , \(y\) i approximate the eigenvectors of \(W\) and are called Ritz vectors; The Lanczos acceleration decomposition algorithm combines the Lanczos iteration and the Rayleigh-Ritz method to solve the eigenvalues of the original matrix. The algorithm first constructs an orthogonal matrix \(Q\) k = \([q_1,q_2,......,q\) k of Lanczos vectors by the Lanczos iteration, and then obtains the eigenvalues of , that is, the Ritz values, which are approximate eigenvalues of the original matrix.

[0139] Next, the specific steps of the accelerated decomposition spectral clustering algorithm (CLSC) based on cosine similarity are introduced:

[0140] (1) Read the sample data and convert it into an array, and further reshape the data into a matrix with samples as rows;

[0141] (2) Calculate the similarity matrix W between data samples using the cosine similarity formula;

[0142] (3) Based on the similarity matrix W, normalize it; use the similarity matrix W to calculate the degree matrix D; use the similarity matrix W and the degree matrix D to perform regular normalization to calculate the Laplacian matrix L;

[0143] (4) Introduce the Lanczos acceleration decomposition algorithm during the spectral decomposition of L to obtain the eigenvectors approximating the original matrix;

[0144] (5) Extract the columns y of Q k V as the new feature space, take the eigenvectors u1, u2, …, u corresponding to the top k largest eigenvalues i , construct the matrix U, U = (u1, u2, …, u k ), ∈ R k ; n×k ;

[0145] (6) Normalize the row vectors of the matrix U to obtain the matrix Z ∈ R N×k :

[0146] (7) Take each row of the matrix z as a point in the space and perform K - means++ clustering to obtain the clustering result;

[0147] Finally, obtain the partition result based on the clustering result.

[0148] The present invention also sets up a traceability mechanism. Through metadata management, a data lineage tracking system is established. Any partition parameter can be traced back to: the original data file (timestamp + device number), the version of the processing algorithm (Git commit ID), the operator information, etc.

[0149] In step 3, for each partition, optimize the blasting parameters to obtain the blasting parameters, specifically:

[0150] Step 301: Construct a data set;

[0151] Step 302: Construct a blasting parameter optimization model based on the Attention - Stacked LSTM network structure;

[0152] Step 303: Train the blasting parameter optimization model based on the data set;

[0153] Step 304: Input the input parameters of each partition into the trained blasting parameter optimization model to obtain the optimized blasting parameters.

[0154] In step 301, the data set includes the partition-average RII, the partition-maximum joint density J dmax , the minimum distance D from the sensitive area min , the buried depth H, the average value of historical blasting PPV, and the rock stratum dip angle θ. The specific composition is shown in Table 1;

[0155] Table 1 Composition table of the data set

[0156]

[0157] In step 302, a blasting parameter optimization model is constructed based on the Attention-Stacked LSTM network structure. Specifically:

[0158] As shown in Table 1, the input parameters of the present invention are the partition-average RII, the partition-maximum joint density J dmax , the minimum distance D from the sensitive area min , the buried depth H, the average value of historical blasting PPV, and the rock stratum dip angle θ, all of which can be obtained from the collected data;

[0159] The output parameters of the present invention are blasting parameters, including: single-hole charge Q (kg), hole spacing S (m), and millisecond delay time Δt (ms);

[0160] Next, the Attention-Stacked LSTM network structure will be introduced in detail:

[0161] Each module in the network uses LSTM as the working unit. Stacking two LSTM units enhances data processing and, by introducing the attention mechanism, largely solves the problem of long-distance information loss;

[0162] First, the attention mechanism will be introduced:

[0163] The attention mechanism Attention in the neural network aims to, when the computing resources are limited, allocate weights to all computing information through weights, and then allocate computing resources according to the weight distribution, while solving the problem of information overload. It is a resource allocation scheme.

[0164] 1. Attention distribution

[0165] The input information volume X can be regarded as an information memory for storing information. Given a query vector q, query the information in this form and select a certain information in the input information volume X. During the query and selection process, it is necessary to know the index position of the selected information. A "soft" mechanism is used for selection, and a part of the information is extracted from all the input information as the input of the attention mechanism, and the extraction is carried out proportionally according to the relevance.

[0166]

[0167] Among them, an attention variable is defined represents the index position of the selected information. When i = 1, it means the first data is input. Given X and q, α i represents the probability of the i-th input information. The probability formed by α i is called the attention distribution. s(x i , q) is the attention scoring function. This model includes an additive model, a dot product model, a scaled dot product model, and a bilinear model. The formulas are as follows:

[0168] s(x i , q) = v T tanh(Wx i + Uq);

[0169]

[0170] 2. Weighted average of attention

[0171] The attention distribution α i represents the degree of relevance between the i-th information in the input information vector X and the query q when the query q is given. A "soft" information selection mechanism is adopted to give the result obtained by the query. This process is to summarize the input information in a weighted average manner to obtain the Attention value.

[0172]

[0173] 3. Key-value pair attention pattern

[0174] In the attention data input mechanism, a pair of key-value pairs can be used to represent a pair of keys and values for information input, and then N input messages can represent information input. The "key" is used to calculate the attention distribution, and the "value" is used to calculate the information summary. Then, it can be regarded as a soft processing mechanism in the attention process: the input message x stores the content in the memory, and the element consists of a key and a value. Currently, there is a key = query, and the goal is to retrieve the corresponding Value value in the memory, including obtaining the corresponding memory value, that is, the Attention value. In soft addressing, it is not necessary to strictly satisfy the condition of Key = Query to retrieve the stored information, but to calculate the similarity between the query Query and the address Key of the item in the memory to determine how much content should be retrieved from the corresponding Value item. The value corresponding to each address Key will be extracted and summed, which is equivalent to calculating the weight of each value according to the similarity between the query and the keyword, and then performing a weighted sum of these values. The weighted sum obtains the final value of the value, that is, the value of attention;

[0175] The above calculations can be summarized into three processes:

[0176] (1) Calculate the similarity between Query and Key. It can be calculated using the additive model, dot product model, or cosine similarity listed above to obtain the attention score s i ;

[0177] s i = F(Q, k i );

[0178] (2) Numerically transform the attention score using the softmax function. On the one hand, it can be normalized to obtain a probability distribution with the sum of all weight coefficients equal to 1, and on the other hand, the characteristics of the softmax function can be used to highlight the weights of important elements;

[0179]

[0180] (3) Perform a weighted sum of Value according to the weight coefficients:

[0181]

[0182] According to the above formula, we get:

[0183]

[0184] According to the formula change, it is the mathematical change principle of the soft attention mechanism.

[0185] 4. Attention Mechanism and Encoder-Decoder Framework

[0186] As a general idea, the attention mechanism itself does not depend on the framework. At present, it is used in combination with the Encoder-Decoder framework in most cases. Figure 2 shown.

[0187] Among them, X t is the input sequence; Y t is the output sequence, C i is the weighted sum of the hidden quantity h in the encoder. Weight coefficient a ij and the encoder state h at each moment j And the decoder state h' at the previous moment i-l The input sequence is calculated by the encoder in the Encoder and then input into the Decoder, and the output is finally obtained after weight calculation.

[0188] 5. Self-Attention Mechanism

[0189] The traditional Attention mechanism occurs between the elements of the target and all the elements of the source. That is, in the Encoder-Decoder model, the calculation of the Attention weight requires not only the hidden state in the Encoder but also the hidden state in the Decoder. In the actual training process, the input received by the neural network is a lot of vectors of different sizes, and there is a certain relationship between different vectors. However, in the actual model training, the relationship between these input vectors cannot be fully utilized, which also leads to poor training model effect.

[0190] In practical applications, the sequences obtained are often of unequal length. In order to make the sequence length consistent when inputting the network, convolutional or recurrent neural network encoding is often used. However, both methods use local information in the sequence to encode. In order to obtain long-distance dependencies in the sequence, other encoding methods need to be sought. Using the attention mechanism to "dynamically" generate different connection weights, the influence between different events in the sequence can be directly obtained, and long-distance dependencies between events can be easily established. This is the self-attention mechanism model, which has been successfully applied in many fields.

[0191] The self-attention mechanism usually uses the query-key-value model, assuming that the input sequence has been obtained. The output sequence is The calculation process of the self-attention mechanism model is as follows:

[0192] (1) For each input vector x i, first map x through three different linear transformations i to three different linear spaces, and the corresponding results are: the query vector the key vector the value vector For the input vector sequence X, the matrices Q = [q1; q2;..; q N composed of query vectors, the matrix K = [k1; k2;..; k N composed of key vectors, and the matrix V = [v1; v2;..; v N composed of value vectors are obtained through three different linear transformations. The formulas for Q, K, and V are as follows:

[0193]

[0194] where are the parameter matrices obtained during different linear transformations respectively.

[0195] (2) When the query vector is determined as q n ∈Q, the output vector h n can be obtained through a series of calculations:

[0196]

[0197] where n, j ∈ [1, 2,..., N] are the positions of the input sequence and the output sequence, and a nj represents the weight that the j-th input receives attention from the n-th input.

[0198] 6. Multi-Head Self-Attention Mechanism

[0199] Since the self-attention mechanism was proposed, it has been studied by many scholars. The self-attention mechanism was proposed to replace convolutional and recurrent neural networks. Generally, when using the self-attention mechanism, an input is mapped to different projection spaces through different linear transformations to extract more available information and make the model achieve better results. This is the multi-head self-attention. Given an input vector sequence encoded with the self-attention mechanism, when the attention scoring function s(q, x) selects the dot product model, it is expressed as:

[0200]

[0201] where the self-attention encoding can be regarded as establishing an interaction relationship between different vectors in the input vector X. To extract more available interaction information, the multi-head attention mechanism can be used to map the input vector X to different projection spaces and then capture different interaction information. Suppose the self-attention encoding is applied in M different projection spaces respectively, then there is

[0202] MultiHead(H) = (head1;...; head M )W o ;

[0203] head m = satt(Q m , K m , V m );

[0204]

[0205] where, m ∈ {1,…, M}.

[0206] 7. Stacked Long Short-Term Memory Networks

[0207] The success of deep neural networks is usually attributed to the hierarchical structure caused by multiple layers. In each layer, a certain part of the problem is solved, and the results and the remaining problems are passed to the next layer for further processing. Deep learning is based on the assumption that a deep, hierarchical model can represent certain functions better than some large shallow models. Additional hidden layers can be added to a multi-layer perceptron neural network to deepen the network. The additional hidden layers are understood to re-learn the representation of the combined results from the previous layers and create new representations at a high level of abstraction.

[0208] Stacked LSTM is now a stable technology for challenging sequence prediction problems. The stacked LSTM architecture can be defined as an LSTM model composed of multiple LSTM layers. The upper LSTM layer provides a sequence output instead of a single-value output to the lower LSTM layer. Specifically, there is one output for each input time step, rather than one output time step for all input time steps. The stacked LSTM hidden layer makes the model deeper and can be more accurately described as a deep learning technology. The structure of the stacked LSTM model is as Figure 3 shown.

[0209] 8. Introduction to the construction of the Attention-StackedLSTM network structure:

[0210] Take an LSTM network as the encoder to process the input sequence and encode the entire sequence into a fixed-length vector, that is, the input vector; then use a custom Attention layer as the decoder, which focuses on learning key features and performs a fully connected process on the hidden layer features based on the key features to improve the model's learning ability, thus forming an input unit, and adding another unit to form the Attention-StackedLSTM network structure. Its structure diagram is as Figure 4 shown.

[0211] Among them, X t is the input sequence; Y t is the output sequence; h j is the hidden state, and the weight of each hidden state h j is a i j, indicating the influence of this feature on the output result. Due to the introduction of the Attention mechanism, the upper and lower data c j needs to comprehensively calculate the weight information and the hidden information, such as i as

[0212]

[0213] In step 4, based on the blasting parameters of each area, blast zone controlled blasting is carried out, specifically:

[0214] According to the obtained single-hole charge amount Q, hole spacing S, and millisecond delay time Δt, intelligent hole layout and precise charging are carried out for each zone, and a timing control initiation is carried out by adopting an electronic detonator networking system according to specific requirements to achieve blast zone controlled blasting.

[0215] In this specification, each embodiment is described in a progressive manner. What each embodiment focuses on explaining is the difference from other embodiments. The same or similar parts among the embodiments can be referred to each other.

[0216] In this article, specific examples are used to elaborate on the principle and implementation manner of the present invention. The descriptions of the above embodiments are only used to help understand the method and its core idea of the present invention; at the same time, for those of ordinary skill in the art, according to the idea of the present invention, there will be changes in the specific implementation manner and application scope. In summary, the content of this specification should not be construed as a limitation to the present invention.

Claims

1. A method for controlled blasting in blasting zones of long tunnels in environmentally sensitive areas, characterized in that, Including: Step 1: Collect surface data, tunnel data and auxiliary data, and preprocess them; Step 2: Conduct intelligent rock mass zoning based on the preprocessed data; Step 3: Optimize blasting parameters for each zone to obtain blasting parameters; Step 4: Conduct blasting zone controlled blasting based on the blasting parameters of each zone.

2. The method according to claim 1, characterized in that In Step 1, the surface data, tunnel data and auxiliary data are collected, specifically: Collect surface data based on unmanned aerial vehicle lidar and multispectral cameras. Among them, the surface data includes terrain elevation model, surface lithology spectral characteristics, surface fracture distribution density and GPS coordinates of surface sensitive targets; Collect tunnel data based on 3D laser scanners, ground-penetrating radars and digital borehole cameras. Among them, the tunnel data includes face rock mass hardness, joint density, rock layer dip / strike, groundwater permeability coefficient and 3D fracture network model; Obtain auxiliary data based on the historical database. Among them, the auxiliary data includes mechanical parameters of geological borehole cores, historical blasting vibration monitoring records, regional seismic wave propagation characteristics and explosive performance parameter library.

3. The method according to claim 2, wherein In Step 1, the surface data, tunnel data and auxiliary data are preprocessed, specifically: The surface curvature C is obtained by slope / curvature calculation and terrain roughness index calculation based on the terrain elevation model surf ; Based on the surface lithology spectral characteristics, through multispectral band fusion and support vector machine lithology classification, obtain lithology codes; Calculating the surface fissure density F through fissure extraction from UAV images and density kernel estimation based on the distribution density of surface fissures surf ; Based on the GPS coordinates of surface sensitive targets, calculate the distance D from the sensitive area through spatial topological relationship calculation and buffer analysis; Based on the 3D fracture network model, obtain the standard deviation σ of fracture spacing through point cloud skeleton extraction and fracture spacing statistical analysis.

4. The method according to claim 3, characterized in that, In Step 2, conduct intelligent rock mass zoning based on the preprocessed data, specifically: Step 201: Standardize the face rock mass hardness f, groundwater permeability coefficient K and the distance D from the sensitive area; Step 202: Calculate the rock mass integrity index RII, environmental sensitivity coefficient ES and fracture development index FDI; Step 203: Based on the rock mass integrity index RII, environmental sensitivity coefficient ES, fracture development index FDI, the longitudinal wave velocity Vp in the regional seismic wave propagation characteristics and the historical PPV mean value in the historical blasting vibration monitoring records, conduct automatic division of intelligent rock mass zones through the spectral clustering algorithm improved based on phase velocity.

5. The method according to claim 4, characterized in that, In Step 201, the face rock mass hardness f, groundwater permeability coefficient K and the distance D from the sensitive area are standardized, specifically: Standardize the face rock mass hardness f based on the Min-Max normalization method; Standardize the groundwater permeability coefficient K based on the logarithmic transformation normalization method; Standardize the distance D from the sensitive area based on piecewise function mapping.

6. The method according to claim 5, characterized in that, In Step 202, calculate the rock mass integrity index RII, environmental sensitivity coefficient ES and fracture development index FDI, specifically: The rock mass integrity index RII is calculated through the hardness f of the rock mass at the tunnel face and the joint density J d and is given by: RII = 0.6tanh(0.1f) + 0.4exp(-0.02J d ); The environmental sensitivity coefficient ES is calculated through the distance D from the sensitive area and the surface curvature C surf and is obtained as follows: ES = 0.7×(1 / D) + 0.3×C surf ; The fracture development index FDI is calculated through the standard deviation σ of fracture spacing and the rock layer dip θ, and is: FDI = σ / (1 + θ / 90).

7. The method according to claim 6, characterized in that, In step 203, based on the rock mass integrity index RII, the environmental sensitivity coefficient ES, the fracture development index FDI, the longitudinal wave velocity Vp in the regional seismic wave propagation characteristics, and the historical PPV mean value in the historical blasting vibration monitoring records, an intelligent partition of the rock mass is automatically divided by the spectral clustering algorithm improved based on the phase velocity, specifically: Determine that the input vector is the rock mass integrity index RII, the environmental sensitivity coefficient ES, the fracture development index FDI, the longitudinal wave velocity Vp in the regional seismic wave propagation characteristics, and the historical PPV mean value in the historical blasting vibration monitoring records, that is: V = [RII, ES, FDI, V p , historical PPV]; Based on the input vector V, an intelligent partition of the rock mass is automatically divided by the accelerated decomposition spectral clustering algorithm based on the cosine similarity.

8. The method according to claim 7, wherein In step 3, for each partition, optimize the blasting parameters to obtain the blasting parameters, specifically: Construct a data set; Construct a blasting parameter optimization model based on the Attention-Stacked LSTM network structure; Train the blasting parameter optimization model based on the data set; Input the input parameters of each partition into the trained blasting parameter optimization model to obtain the optimized blasting parameters.

Images