Rock mass damage radius identification method and device based on blasting vibration signal

Abstract

The invention discloses a rock mass damage radius identification method and device based on a blasting vibration signal. The method comprises the following steps: acquiring an initial blasting vibration signal; preprocessing the initial blasting vibration signal to obtain a preprocessed blasting vibration signal; performing parameter optimization processing by utilizing a sparrow search method according to the preprocessed blasting vibration signals to obtain variational mode decomposition parameters; calculating a modal component to be evaluated; component evaluation indexes are calculated; calculating a target modal component; constructing a target blasting vibration signal; constructing a joint multi-domain feature vector; performing principal component analysis on the joint multi-domain feature vector to obtain a joint multi-domain feature matrix; constructing a feature information data set; according to the feature information data set, a sparrow search method is utilized to optimize the decision tree model, and a rock mass damage radius recognition model is obtained; and performing rock mass damage radius identification by using the rock mass damage radius identification model. According to the invention, rock mass damage radius identification is realized, and the efficiency is improved. The method can be widely applied to the technical field of engineering blasting.

Description

Technical Field

[0001] The invention relates to the technical field of engineering blasting, and in particular to a method and a device for identifying rock mass damage radius based on blasting vibration signals. Background Art

[0002] With the development of blasting technology, blasting has been widely used in mining, subway tunnels and other infrastructure construction fields. The blasting vibration signal contains information such as the scale of the blasting, the nature of the rock, and the impact on other surrounding structures. During the blasting process, the stress wave generated by the explosion will cause certain additional damage to the rock mass near the blasting area. In order to accurately evaluate the effect of the blasting and the degree of impact of the blasting on the surrounding environment, it is necessary to determine the radius of rock damage caused by the blasting. The traditional method uses sound wave testing to detect changes in the properties of the rock mass through the wave velocity and amplitude of the sound wave. However, this method requires monitoring the wave velocity changes before and after the rock mass blasting, and the process is cumbersome and inefficient.

[0003] In summary, the technical problems existing in the relevant technologies need to be improved. Summary of the invention

[0004] The embodiment of the present invention provides a method and a device for identifying rock damage radius based on blasting vibration signals, which effectively improves the identification efficiency.

[0005] On the one hand, an embodiment of the present invention provides a method for identifying rock damage radius based on blasting vibration signals, comprising the following steps:

[0006] Acquire initial blasting vibration signal;

[0007] Preprocessing the initial blasting vibration signal to obtain a preprocessed blasting vibration signal;

[0008] According to the pre-processed blasting vibration signal, a parameter optimization process is performed using a sparrow search method to obtain a variational modal decomposition parameter, wherein the variational modal decomposition parameter includes a mode number and a penalty factor;

[0009] Calculating the modal components to be evaluated according to the variational modal decomposition parameters;

[0010] Calculating component evaluation indicators according to the modal component to be evaluated, wherein the component evaluation indicators include multi-scale fuzzy entropy, first variance, first skewness, first kurtosis or correlation coefficient;

[0011] Calculating the target modal component according to the component evaluation index;

[0012] constructing a target blasting vibration signal according to the target modal component;

[0013] constructing a joint multi-domain feature vector according to the target blasting vibration signal;

[0014] Performing principal component analysis on the joint multi-domain feature vector to obtain a joint multi-domain feature matrix;

[0015] Constructing a characteristic information data set, the characteristic information data set including the joint multi-domain characteristic matrix and a rock damage radius, the rock damage radius being determined by comparing a preset change threshold with a wave velocity change obtained by an acoustic wave device;

[0016] According to the characteristic information data set, the decision tree model is optimized by using the sparrow search method to obtain a rock damage radius identification model;

[0017] The rock mass damage radius identification model is used to identify the rock mass damage radius.

[0018] In some embodiments, the preprocessing of the initial blasting vibration signal to obtain a preprocessed blasting vibration signal includes:

[0019] Get the time node of blasting implementation;

[0020] Comparing the blasting implementation time node with the initial blasting vibration signal in time to obtain a blasting vibration signal period;

[0021] The pre-processed blasting vibration signal is extracted from the initial blasting vibration signal according to the blasting vibration signal period.

[0022] In some embodiments, the method of performing parameter optimization processing using a sparrow search optimization method according to the preprocessed blasting vibration signal to obtain variational mode decomposition parameters includes:

[0023] Setting parameters of the sparrow search method, the parameters of the sparrow search method include the value range of the number of modes, the value range of the penalty factor, the population size, the maximum number of iterations, the forager ratio, the guard ratio, the step length control parameter, the early warning movement amplitude coefficient, the disturbance coefficient or the random direction quantity;

[0024] Generating a search population according to the parameters of the sparrow search method, wherein the types of individuals in the search population include foragers, followers and guards;

[0025] Constructing a sparrow search objective function, wherein the optimization target of the sparrow search objective function is to minimize the mean square error between the sum of the modal components to be optimized and the preprocessed blasting vibration signal;

[0026] Constructing a forager position update function according to the search population, the population convergence rate constant, the maximum number of iterations and the direction vector;

[0027] Constructing a follower position update function according to the search population and the random disturbance value;

[0028] Constructing a sentinel position update function according to the search population and the moving speed constant;

[0029] The variational mode decomposition parameters are calculated according to the sparrow search objective function, the forager position update function, the follower position update function and the sentinel position update function.

[0030] In some embodiments, the calculating the modal component to be evaluated according to the variational modal decomposition parameter includes:

[0031] Constructing a variational mode decomposition objective function according to the number of modes, the modal components to be optimized, the center frequency and the time differential operator;

[0032] Constructing conditional constraints according to the modal components to be optimized;

[0033] Constructing a Lagrangian function according to the variational mode decomposition objective function, the conditional constraints, the Lagrangian multiplier and the penalty factor;

[0034] Constructing a modal component update function to be optimized according to the Lagrangian function, the frequency variable and the penalty factor;

[0035] Constructing a center frequency update function according to the modal component to be optimized and the frequency variable;

[0036] Constructing a Lagrange multiplier update function according to the modal components to be optimized and the conditional constraints;

[0037] The modal component to be evaluated is calculated using an alternating direction multiplier method according to the modal component update function to be optimized, the center frequency update function and the Lagrange multiplier update function.

[0038] In some embodiments, the calculating the component evaluation index according to the modal component to be evaluated includes:

[0039] Calculating a coarsening sequence according to the modal components to be evaluated and the scale factors;

[0040] constructing an embedding vector according to a first preset embedding vector dimension, the coarsening sequence and the scale factor;

[0041] According to the embedding vector, calculating the embedding vector distance;

[0042] Calculating fuzzy membership according to the embedding vector distance;

[0043] Calculating a first average value according to the fuzzy membership degree;

[0044] Calculating a second average value according to a second preset embedding vector dimension, where the second preset embedding vector dimension is the first preset embedding vector dimension plus 1;

[0045] Calculating the multi-scale fuzzy entropy according to the first average value and the second average value;

[0046] Calculating the first variance according to the modal component to be evaluated and the modal component average value;

[0047] Calculating the first skewness according to the modal component to be evaluated, the sample mean, the expected value operator and the first variance;

[0048] Calculating the first kurtosis according to the modal component to be evaluated, the sample mean, the expected value operator and the first variance;

[0049] The correlation coefficient is calculated according to the modal component to be evaluated, the preprocessed blasting vibration signal, the covariance calculation function and the standard deviation operator.

[0050] In some embodiments, the calculating the target modal component according to the component evaluation index includes:

[0051] Calculating the evaluation amount according to the component evaluation index;

[0052] Calculate the comprehensive evaluation value according to the evaluation amount and the total number of component evaluation indicators;

[0053] The target modal component is selected from the modal components to be evaluated according to a preset evaluation threshold and the comprehensive evaluation value.

[0054] In some embodiments, constructing a target blasting vibration signal according to the target modal component includes:

[0055] Frequency modulating the target modal component to obtain a unit amplitude analytical signal;

[0056] Calculating a pseudo-distributed spectrum according to the unit amplitude analytical signal, the weighting function, the delay variable, the complex unit and the first frequency;

[0057] Calculating an original effective signal amplitude estimate based on the pseudo-distribution spectrum and the proportionality coefficient;

[0058] According to the scaling factor, the translation factor and the normalization factor, a scaling and translation wavelet function is constructed;

[0059] According to the scaling and translation wavelet function and the conjugate mother wavelet function, the original effective signal amplitude estimation is subjected to wavelet transformation to obtain wavelet transformation coefficients;

[0060] According to the wavelet transform coefficients and the sign function, a soft threshold method is used to perform wavelet analysis on the target modal component to obtain approximate coefficients and detail coefficients;

[0061] Performing wavelet reconstruction on the target modal component according to the approximate coefficient and the detail coefficient to obtain a reconstructed signal;

[0062] Calculate the matrix dimension based on the total number of modal components;

[0063] constructing an initial Hankel matrix according to the reconstructed signal and the matrix dimension;

[0064] Calculating a weighted energy contribution rate according to the initial Hankel matrix;

[0065] Determining the effective order of singular values ​​according to the weighted energy contribution rate;

[0066] According to the effective order of the singular values, the orthogonal matrix and the non-negative diagonal matrix, the initial Hankel matrix is ​​subjected to singular value decomposition to obtain a target Hankel matrix;

[0067] The target blasting vibration signal is constructed according to the target Hankel matrix.

[0068] In some embodiments, constructing a joint multi-domain feature vector according to the target blasting vibration signal includes:

[0069] Calculating a root mean square according to the target blasting vibration signal;

[0070] Calculating a second variance according to the target blasting vibration signal;

[0071] Calculating a second skewness according to the target blasting vibration signal, the sample mean and the second variance;

[0072] Calculating a second kurtosis according to the target blasting vibration signal, the sample mean and the second variance;

[0073] constructing a time domain feature vector according to the root mean square, the second variance, the second skewness and the second kurtosis;

[0074] Calculating spectral skewness according to the target blasting vibration signal, frequency variable, average frequency and frequency standard deviation;

[0075] Calculating spectral kurtosis according to the target blasting vibration signal, the frequency variable, the average frequency and the frequency standard deviation;

[0076] Constructing a frequency domain feature vector according to the spectral skewness and the spectral kurtosis;

[0077] Performing wavelet packet decomposition on the target blasting vibration signal to obtain a plurality of wavelet packet decomposition coefficients;

[0078] Reconstructing the wavelet packet decomposition coefficients to obtain frequency band sub-signals;

[0079] Calculating signal energy according to the frequency band sub-signals;

[0080] Calculating total energy based on multiple signal energies;

[0081] Calculating energy probability according to the signal energy and the total energy;

[0082] Calculating energy entropy based on the energy probability and the total signal energy;

[0083] Constructing a time-frequency domain feature vector according to the plurality of energy entropies;

[0084] The joint multi-domain feature vector is constructed according to the time domain feature vector, the frequency domain feature vector and the time-frequency domain feature vector.

[0085] In some embodiments, performing principal component analysis on the joint multi-domain feature vector to obtain a joint multi-domain feature matrix includes:

[0086] Normalizing the joint multi-domain feature vector to obtain a normalized vector;

[0087] Calculate the covariance matrix based on the sample size and the standardized vector;

[0088] Performing eigenvalue decomposition on the covariance matrix to obtain eigenvalues ​​to be selected and eigenvectors to be selected;

[0089] Calculate the cumulative variance contribution rate according to the eigenvalues ​​to be selected and the total number of eigenvalues;

[0090] Selecting a target eigenvalue from the eigenvalues ​​to be selected according to the cumulative variance contribution rate;

[0091] Selecting a target eigenvector from the eigenvectors to be selected according to the cumulative variance contribution rate and the target eigenvalue;

[0092] The joint multi-domain feature matrix is ​​constructed according to the target feature vector and the joint multi-domain feature vector.

[0093] On the other hand, an embodiment of the present invention provides a rock damage radius identification device based on blasting vibration signals, comprising:

[0094] The first module is used to obtain the initial blasting vibration signal;

[0095] The second module is used to preprocess the initial blasting vibration signal to obtain a preprocessed blasting vibration signal;

[0096] The third module is used to perform parameter optimization processing using the sparrow search method according to the pre-processed blasting vibration signal to obtain variational modal decomposition parameters, wherein the variational modal decomposition parameters include the number of modes and a penalty factor;

[0097] The fourth module is used to calculate the modal component to be evaluated according to the variational modal decomposition parameters;

[0098] A fifth module is used to calculate a component evaluation index according to the modal component to be evaluated, wherein the component evaluation index includes a multi-scale fuzzy entropy, a first variance, a first skewness, a first kurtosis or a correlation coefficient;

[0099] A sixth module is used to calculate the target modal component according to the component evaluation index;

[0100] A seventh module is used to construct a target blasting vibration signal according to the target modal component;

[0101] An eighth module is used to construct a joint multi-domain feature vector according to the target blasting vibration signal;

[0102] A ninth module is used to perform principal component analysis on the joint multi-domain feature vector to obtain a joint multi-domain feature matrix;

[0103] A tenth module is used to construct a feature information data set, wherein the feature information data set includes the joint multi-domain feature matrix and a rock damage radius, and the rock damage radius is obtained by an acoustic wave device;

[0104] The eleventh module is used to optimize the decision tree model using the sparrow search method according to the characteristic information data set to obtain a rock damage radius identification model;

[0105] The twelfth module is used to identify the rock damage radius using the rock damage radius identification model.

[0106] The beneficial effects of the present invention are as follows:

[0107] The embodiment of the present invention first obtains an initial blasting vibration signal, preprocesses the initial blasting vibration signal to obtain a preprocessed blasting vibration signal, then uses a sparrow search method to perform parameter optimization processing based on the preprocessed blasting vibration signal to obtain variational modal decomposition parameters, calculates modal components to be evaluated, calculates component evaluation indicators, calculates target modal components, constructs target blasting vibration signals, and constructs a joint multi-domain feature vector. Then, principal component analysis is performed on the joint multi-domain feature vector to obtain a joint multi-domain feature matrix, and a feature information data set is constructed. Finally, based on the feature information data set, a decision tree model is optimized using a sparrow search method to obtain a rock damage radius identification model, and the rock damage radius identification model is used to identify the rock damage radius, so that the rock damage radius identification can be realized by analyzing the blasting vibration signal, thereby improving efficiency.

[0108] Other features and advantages of the present invention will be described in the following description, and partly become apparent from the description, or understood by practicing the present invention. The purpose and other advantages of the present invention can be realized and obtained by the structures particularly pointed out in the description and the drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0109] In order to more clearly illustrate the technical solutions in the embodiments of the present application, the drawings required for use in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0110] Figure 1 This is a flow chart of a method for identifying rock damage radius based on blasting vibration signals according to an embodiment of the present invention;

[0111] Figure 2 This is a waveform diagram of a vibration signal in the x direction at a certain measuring point in an embodiment of the present invention;

[0112] Figure 3 A schematic diagram of the connection relationship between a blasting vibration tester, an acoustic wave instrument and a blasting area according to an embodiment of the present invention;

[0113] Figure 4 This is a schematic diagram of eight modal components to be evaluated obtained in an embodiment of the present invention;

[0114] Figure 5 A schematic diagram of a comprehensive evaluation value of a modal component to be evaluated according to an embodiment of the present invention;

[0115] Figure 6 It is a schematic diagram of an overall process of analyzing blasting vibration signals and obtaining rock damage radius through model identification according to an embodiment of the present invention;

[0116] Figure 7 The present invention is a schematic structural diagram of a rock damage radius identification device based on blasting vibration signals according to an embodiment of the present invention. DETAILED DESCRIPTION

[0117] In order to make the purpose, technical solutions and advantages of the present application clearer, the present application is further described in detail below in conjunction with the accompanying drawings and examples. It should be understood that the specific embodiments described herein are only used to explain the present application and are not intended to limit the present application. When the following description refers to the drawings, unless otherwise indicated, the same numbers in different drawings represent the same or similar elements. The embodiments described in the following exemplary embodiments do not represent all embodiments consistent with the embodiments of the present application. They are only examples of devices and methods consistent with some aspects of the embodiments of the present application as detailed in the attached claims.

[0118] It is understood that the terms "first", "second", etc. used in this application can be used to describe various concepts in this article, but unless otherwise specified, these concepts are not limited by these terms. These terms are only used to distinguish one concept from another concept. For example, without departing from the scope of the embodiment of the present application, the first information may also be referred to as the second information, and similarly, the second information may also be referred to as the first information. Depending on the context, the words "if" and "if" as used herein can be interpreted as "at the time of" or "when" or "in response to determination".

[0119] The terms "at least one", "multiple", "each", "any", etc. used in this application, at least one includes one, two or more, multiple includes two or more, each refers to each of the corresponding multiple, and any refers to any one of the multiple.

[0120] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as those commonly understood by those skilled in the art to which this application belongs. The terms used herein are only for the purpose of describing the embodiments of this application and are not intended to limit this application.

[0121] Before describing the embodiments of the present application in detail, some nouns and terms involved in the embodiments of the present application are first described. The nouns and terms involved in the embodiments of the present application are subject to the following explanations.

[0122] Blasting vibration signal: refers to the phenomenon that after the explosion of explosives, part of the energy propagates outward in the form of elastic waves, causing the rock mass to vibrate at a certain location. This vibration signal can be recorded and analyzed by monitoring equipment to evaluate its safety impact on the surrounding rock mass and structure.

[0123] In the related technology, with the development of blasting technology, blasting has been widely used in mining, subway tunnels and other infrastructure construction fields. During the blasting process, the stress wave generated by the explosion will cause certain additional damage to the rock mass near the blasting area. In order to accurately evaluate the effect of blasting and the impact of blasting on the surrounding environment, it is necessary to determine the radius of rock damage caused by blasting, and the design of blasting schemes can be provided based on the radius of rock damage. Therefore, it is of great significance to determine the range of rock damage caused by blasting. At present, acoustic wave testing is most widely used in identifying rock damage radius because of its fast and convenient characteristics. The principle is that the propagation characteristics of sound waves in different media are different. When the properties of the rock mass change, the wave velocity and amplitude of the sound wave will respond. However, this method requires monitoring the wave velocity changes before and after the rock blasting, and the process is relatively cumbersome and inefficient.

[0124] In view of this, the embodiment of the present invention utilizes the information contained in the blasting vibration signal, such as the blasting scale, the nature of the rock, the impact on other surrounding structures, etc., to determine the rock damage radius after the blasting by analyzing the information in the blasting vibration signal, thereby making the identification of the rock damage radius simple and efficient.

[0125] The rock damage radius identification method based on blasting vibration signals provided in the embodiment of the present application relates to the technical field of engineering blasting. The rock damage radius identification method based on blasting vibration signals provided in the embodiment of the present application can be applied to a terminal, can also be applied to a server, and can also be software running in a terminal or a server. In some embodiments, the terminal can be a smart phone, a tablet computer, a laptop computer, a desktop computer, a smart speaker, a smart watch, and a car terminal, etc., but is not limited to this; the server side can be configured as an independent physical server, or a server cluster or distributed system composed of multiple physical servers, and can also be configured as a cloud server that provides basic cloud computing services such as cloud services, cloud databases, cloud computing, cloud functions, cloud storage, network services, cloud communications, middleware services, domain name services, security services, CDN, and big data and artificial intelligence platforms. The server can also be a node server in a blockchain network; the software can be an application that implements the rock damage radius identification method based on blasting vibration signals, etc., but is not limited to the above forms.

[0126] The present application can be used in many general or special computer system environments or configurations. For example: personal computers, server computers, handheld or portable devices, tablet devices, multiprocessor systems, microprocessor-based systems, set-top boxes, programmable consumer electronics, network PCs, minicomputers, mainframe computers, distributed computing environments including any of the above systems or devices, etc. The present application can be described in the general context of computer-executable instructions executed by a computer, such as program modules. Generally, program modules include routines, programs, objects, components, data structures, etc. that perform specific tasks or implement specific abstract data types. The present application can also be practiced in distributed computing environments, in which tasks are performed by remote processing devices connected through a communication network. In a distributed computing environment, program modules can be located in local and remote computer storage media including storage devices.

[0127] The following is a detailed explanation of the embodiments of the present application in conjunction with the accompanying drawings:

[0128] Figure 1 is an optional flow chart of a rock damage radius identification method based on blasting vibration signals provided in an embodiment of the present application. Figure 1 The method may include but is not limited to steps S101 to S112.

[0129] Step S101, obtaining an initial blasting vibration signal;

[0130] Step S102, preprocessing the initial blasting vibration signal to obtain a preprocessed blasting vibration signal;

[0131] Step S103, according to the pre-processed blasting vibration signal, a sparrow search method is used to perform parameter optimization processing to obtain variational modal decomposition parameters, where the variational modal decomposition parameters include the number of modes and a penalty factor;

[0132] Step S104, calculating the modal component to be evaluated according to the variational modal decomposition parameters;

[0133] Step S105, calculating component evaluation indicators according to the modal components to be evaluated, the component evaluation indicators including multi-scale fuzzy entropy, first variance, first skewness, first kurtosis or correlation coefficient;

[0134] Step S106, calculating the target modal component according to the component evaluation index;

[0135] Step S107, constructing a target blasting vibration signal according to the target modal component;

[0136] Step S108, constructing a joint multi-domain feature vector according to the target blasting vibration signal;

[0137] Step S109, performing principal component analysis on the joint multi-domain feature vector to obtain a joint multi-domain feature matrix;

[0138] Step S110, constructing a feature information data set, the feature information data set including a joint multi-domain feature matrix and a rock damage radius, the rock damage radius is determined by comparing a preset change threshold with a wave velocity change obtained by an acoustic wave device;

[0139] Step S111, according to the characteristic information data set, the decision tree model is optimized by using the sparrow search method to obtain a rock damage radius identification model;

[0140] Step S112: using the rock mass damage radius identification model to identify the rock mass damage radius.

[0141] Steps S101 to S112 shown in the embodiment of the present application realize the identification of rock damage radius and improve efficiency.

[0142] In step S101 of some embodiments, the initial blasting vibration signal can be obtained by a blasting vibration tester. The initial blasting vibration signal can also be obtained by other methods, but is not limited thereto. For example, the probes can be arranged according to the blasting vibration monitoring scheme, and the blasting vibration tester can monitor the initial blasting vibration signals in three directions of different measuring points to obtain a vibration signal waveform diagram in the x direction of a certain measuring point, as shown in FIG. Figure 2 More specifically, the connection between the blasting vibration tester, the sonic wave instrument and the blasting area is as follows: Figure 3As shown, the computer is connected to the blasting vibration tester and the sonicator, the blasting vibration tester is connected to the blasting vibration measuring point, the sonicator is connected to the sonic hole, there are multiple blasting holes in the blasting area, and the blasting holes can be used as blasting vibration measuring points and sonic holes. It has sonic wave transmitting and receiving equipment, blasting vibration monitoring equipment and data processing and computing equipment. The sonic wave transmitting and receiving equipment is arranged within a certain range around the blasting hole, and is used to measure the sonic wave velocity data between each transmitting-receiving pair before and after the blasting. The blasting vibration monitoring equipment is arranged within a certain range around the blasting area, and is used to monitor the vibration signal generated by the blasting. In the data processing and computing equipment, the sonic wave velocity data measured by the sonic wave transmitting and receiving equipment is calculated by this equipment to calculate the rock damage radius. This equipment can perform variational modal decomposition on the blasting vibration signal obtained by monitoring, and select IMFs with cumulative comprehensive evaluation values ​​exceeding 85% of the total evaluation value through a multi-index comprehensive evaluation method. Then, these IMFs are first subjected to time-frequency peak filtering and then wavelet threshold denoising, and then the processed IMFs are reconstructed, and then singular value decomposition is used for them and then reconstructed. Then, the time domain, frequency domain, and time-frequency domain features of the signal are extracted, and the principal component analysis is used to reduce the dimension of the features. Next, the feature information and its corresponding rock damage radius are stored in the database. Then, some data in the database are selected and divided into training sets and test sets. The GBDT model optimized by the sparrow search algorithm is trained with the training set, and the performance of the model is evaluated using the test set.

[0143] In some embodiments, in step S102, preprocessing the initial blasting vibration signal to obtain a preprocessed blasting vibration signal may include but is not limited to the following steps:

[0144] Get the time node of blasting implementation;

[0145] Compare the blasting implementation time node with the initial blasting vibration signal to obtain the blasting vibration signal period;

[0146] According to the blasting vibration signal period, the pre-processed blasting vibration signal is extracted from the initial blasting vibration signal.

[0147] In some embodiments, the blasting implementation time node (i.e., the time node of the on-site blasting implementation) can be obtained first, and then the blasting implementation time node is compared with the initial blasting vibration signal, and the time node of the recorded initial blasting vibration signal is compared to obtain the blasting vibration signal period, thereby identifying the period of vibration caused by the blasting in the blasting vibration signal. Then, according to the blasting vibration signal period, the pre-processed blasting vibration signal is extracted from the initial blasting vibration signal to eliminate the signal monitored by the blasting vibration tester due to other vibrations.

[0148] In some embodiments, in step S103, according to the preprocessed blasting vibration signal, the sparrow search optimization method is used to perform parameter optimization processing to obtain the variational mode decomposition parameters, which may include but is not limited to the following steps:

[0149] Set the parameters of the sparrow search method, including the range of the number of modes, the range of the penalty factor, the population size, the maximum number of iterations, the forager ratio, the guard ratio, the step size control parameter, the warning movement amplitude coefficient, the disturbance coefficient or the random direction quantity;

[0150] According to the parameters of the sparrow search method, a search population is generated, and the individual types in the search population include foragers, followers and guards;

[0151] Construct a sparrow search objective function, the optimization goal of which is to minimize the mean square error between the sum of the modal components to be optimized and the preprocessed blasting vibration signal;

[0152] Construct the forager position update function according to the search population, population convergence rate constant, maximum number of iterations and direction vector;

[0153] Construct the follower position update function based on the search population and random disturbance value;

[0154] Construct the position update function of the sentinel based on the search population and the moving speed constant;

[0155] According to the sparrow search objective function, the forager position update function, the follower position update function and the guard position update function, the variational mode decomposition parameters are calculated.

[0156] In some embodiments, the sparrow search optimization method can be used to perform parameter optimization processing based on the pre-processed blasting vibration signal to find the optimal parameter combination and obtain the variational modal decomposition parameters, wherein the variational modal decomposition parameters include the number of modes and the penalty factor. The parameters of the sparrow search method can be set first, wherein the parameters of the sparrow search method include the value range of the number of modes K, the value range of the penalty factor α, the population size, the maximum number of iterations, the forager ratio p%, the vigilant ratio a%, the step length control parameter, the warning movement amplitude coefficient, the disturbance coefficient or the random direction quantity. Then, according to the parameters of the sparrow search method, a search population is generated, wherein the individual types in the search population include foragers, followers and vigilants. Then, a sparrow search objective function is constructed, wherein the optimization target of the sparrow search objective function is to minimize the mean square error between the sum of the modal components to be optimized and the pre-processed blasting vibration signal. Exemplarily, a search population containing N sparrow individuals can be generated, and each individual is randomly distributed in the D-dimensional search space. The position of the sparrow individual is represented by X i =(x i1 ,x i1 ,…xid ), and then calculate the fitness value of each individual, that is, construct the sparrow search objective function. The optimization goal of the sparrow search objective function is to minimize the mean square error (MSE) between the sum of the modal components IMF to be optimized after variational mode decomposition (VMD) and the original signal (preprocessed blasting vibration signal). Then, according to the search population, population convergence rate constant, maximum number of iterations and direction vector, the forager position update function is constructed. In the population update strategy, the forager is the individual in the group responsible for exploring new areas, accounting for p% of the search population. The forager position update function is: In the formula, is the position of forager i in generation t+1, is the position of forager i in generation t, α is the population convergence rate constant, which controls the population convergence rate, T max is the maximum number of iterations, Q is a random number that follows a normal distribution and is used to introduce random perturbations, and L is a direction vector that can take values ​​of -1 or 1. Based on the search population and the random perturbation value, a follower position update function is constructed. Followers account for (1-p)% of the search population and mainly rely on foragers for food. Followers update their positions by following foragers. The follower position update function is: In the formula, is the position of follower p in generation t+1, is the position of follower p in generation t, and β is a random number that follows a normal distribution, which is used to introduce random disturbances. According to the search population and the moving speed constant, the position update function of the sentinel is constructed. The sentinel accounts for a% of the search population, and the position update function of the sentinel is: In the formula, is the position of the vigilant o in the t+1 generation, is the position of the alerter o at generation t, S is the moving speed constant, which controls the moving speed of the alerting individual, is a randomly selected sparrow individual with a better position, d oj is the distance between sparrow individuals o and j. Finally, the variational mode decomposition parameters are calculated according to the sparrow search objective function, the forager position update function, the follower position update function and the guard position update function. For example, at the end of each generation of iteration, the fitness values ​​of all individuals (i.e., the sparrow search objective function values) are recalculated, and the current global optimal solution X is updated and recorded. best If the maximum number of iterations is reached, the algorithm terminates, outputs the current optimal solution, and obtains the variational mode decomposition parameters.

[0157] In some embodiments, in step S104, the modal components to be evaluated are calculated according to the variational modal decomposition parameters, which may include but are not limited to the following steps:

[0158] Construct the variational mode decomposition objective function according to the number of modes, the modal components to be optimized, the center frequency and the time differential operator;

[0159] Construct conditional constraints according to the modal components to be optimized;

[0160] Construct a Lagrangian function according to the variational mode decomposition objective function, conditional constraints, Lagrangian multipliers and penalty factors;

[0161] According to the Lagrangian function, frequency variables and penalty factors, the update function of the modal component to be optimized is constructed;

[0162] Construct a center frequency update function according to the modal components and frequency variables to be optimized;

[0163] Construct the Lagrange multiplier update function according to the modal components to be optimized and the conditional constraints;

[0164] According to the modal component update function to be optimized, the center frequency update function and the Lagrange multiplier update function, the modal component to be evaluated is calculated using the alternating direction multiplier method.

[0165] In some embodiments, a variational modal decomposition objective function can be constructed based on the number of modes, the modal components to be optimized, the center frequency, and the time differential operator. Exemplarily, the goal of variational modal decomposition (VMD) is to decompose the preprocessed blasting vibration signal f(t) into K modal components, each of which has a specific center frequency, and these modal components are determined by minimizing the bandwidth. The expression of the variational modal decomposition objective function is: In the formula, is the variational mode decomposition objective function, u k is the kth modal component to be optimized, ω k is the center frequency of the kth modal component to be optimized, K is the number of modes, is the time differential operator, δ(t) is the Dirac delta function, j is the complex unit, t is the time, and * is the convolution operation. Then, according to the modal components to be optimized, conditional constraints are constructed, and the preprocessed blasting vibration signal f(t) can be obtained through all the modal components to be optimized u k (t), the conditional constraint expression is: Then, according to the variational mode decomposition objective function, conditional constraints, Lagrangian multipliers and penalty factors, a Lagrangian function is constructed. For example, in order to solve the optimization problem with constraints, the Lagrangian multiplier λ(t) and the penalty factor α are introduced to construct the Lagrangian function: In the formula, is the Lagrangian function, λ(t) is the Lagrangian multiplier, and α is the penalty factor.

[0166] Then, according to the Lagrangian function, the frequency variable and the penalty factor, a modal component update function to be optimized is constructed. For example, the modal component to be optimized u is updated k (t), keep the center frequency ω k and the Lagrange multiplier λ remains unchanged, minimizing k The Lagrangian function of (t). The expression of the modal component update function to be optimized is: In the formula, is the frequency domain representation of the modal component k to be optimized at the n+1th iteration, is the frequency domain representation of the modal component i to be optimized at the nth iteration, is the Fourier transform of the Lagrange multiplier, is the Fourier transform of the signal, ω is the frequency variable, is the center frequency of the kth modal component to be optimized at the nth iteration. According to the modal component to be optimized and the frequency variable, a center frequency update function is constructed. For example, after the modal component to be optimized is updated, the modal component to be optimized is kept unchanged, and the center frequency of each modal component to be optimized is updated. The expression of the center frequency update function is: In the formula, is the center frequency of the n+1th iteration. According to the modal components to be optimized and the conditional constraints, a Lagrange multiplier update function is constructed. For example, the Lagrange multiplier λ(t) can be updated while ensuring that the conditional constraints are satisfied. The expression of the Lagrange multiplier update function is: In the formula, λ n+1 (t) is the Lagrange multiplier of the n+1th iteration, λ n (t) is the Lagrange multiplier of the nth iteration, and τ is the step size parameter. Finally, according to the modal component update function to be optimized, the center frequency update function and the Lagrange multiplier update function, the alternating direction multiplier method is used to calculate the modal component to be evaluated. Exemplarily, the variational modal decomposition (VMD) solves the optimization problem by the alternating direction multiplier method, and alternately updates the modal components, center frequencies and Lagrange multipliers until the convergence conditions are met. Finally, K modal components to be evaluated and their corresponding center frequencies are output. This embodiment uses the variational modal decomposition optimized by the sparrow search algorithm to process the pre-processed blasting vibration signal obtained by monitoring, and obtains 8 modal components IMF to be evaluated as follows Figure 4 shown.

[0167] In some embodiments, in step S105, the component evaluation index is calculated according to the modal component to be evaluated, which may include but is not limited to the following steps:

[0168] Calculate the coarsening sequence according to the modal components to be evaluated and the scale factor;

[0169] Constructing an embedding vector according to a first preset embedding vector dimension, a coarsening sequence, and a scaling factor;

[0170] According to the embedding vector, calculate the embedding vector distance;

[0171] Calculate the fuzzy membership based on the embedding vector distance;

[0172] Calculate the first average value according to the fuzzy membership degree;

[0173] Calculate a second average value according to a second preset embedding vector dimension, where the second preset embedding vector dimension is the first preset embedding vector dimension plus 1;

[0174] Calculate the multi-scale fuzzy entropy according to the first average value and the second average value;

[0175] Calculate the first variance according to the modal component to be evaluated and the average value of the modal component;

[0176] Calculate the first skewness according to the modal component to be evaluated, the sample mean, the expected value operator and the first variance;

[0177] Calculate the first kurtosis according to the modal component to be evaluated, the sample mean, the expected value operator and the first variance;

[0178] The correlation coefficient is calculated according to the modal component to be evaluated, the preprocessed blasting vibration signal, the covariance calculation function and the standard deviation operator.

[0179] In some embodiments, a component evaluation index can be calculated based on the modal component to be evaluated, wherein the component evaluation index includes multi-scale fuzzy entropy, first variance, first skewness, first kurtosis or correlation coefficient. A coarsening sequence can be calculated based on the modal component to be evaluated and the scale factor. For example, the modal component IMF to be evaluated can be converted into a sequence of multiple different time scales by coarsening. Assume that an IMF component is represented as a time series Where N is the length of the time series. For each scale τ (scale size), the new coarsened sequence is calculated as: In the formula, y (τ) is the coarsening sequence, τ is the scale factor, indicating the degree of coarsening, i is the time index of the time series, j is the index of the time point after coarsening, and y (τ) (j) is the jth data point of the coarsened sequence at the τth scale, and N is the number of data points. Then, an embedding vector is constructed according to the first preset embedding vector dimension, the coarsened sequence and the scale factor. For example, for the swear word sequence y (τ) , construct an m-dimensional embedding vector Y i . Let m be the embedding dimension and r be the similarity tolerance. The embedding vector is calculated as: i =(y(τ) (i),y (τ) (i+1),…,y (τ) (i+m-1)), Where Y i is the embedding vector, and m is the first preset embedding vector dimension. Then, the embedding vector distance is calculated based on the embedding vector, where the calculation formula of the embedding vector distance is: d(Y i ,Y j )=max k＝1,…,m |y (τ) (i+k-1)-y (τ) (j+k-1)|, where d(Y i ,Y j ) is the embedding vector Y i and Y j The embedding vector distance between .

[0180] Then, according to the embedding vector distance, the fuzzy membership is calculated. For example, the fuzzy entropy measures the similarity through the fuzzy membership function. For each pair of embedding vectors Y i and Y j , the expression of fuzzy membership function is: In the formula, μ r (d(Y i ,Y j )) is the fuzzy membership function, r is the similarity tolerance. According to the fuzzy membership, the first average value is calculated. For each i, the average value of the fuzzy membership of the remaining vectors and the fuzzy membership thereof (i.e., the first average value) is calculated, wherein the calculation formula of the first average value is: In the formula, φ m (τ,r) is the first average value. The second average value is calculated according to the second preset embedding vector dimension, where the second preset embedding vector dimension is the first preset embedding vector dimension plus 1. Exemplarily, the embedding dimension is increased from m (the first preset embedding vector dimension) to m+1 (the second preset embedding vector dimension), and the average value calculation process is repeated to obtain φ m+1 (τ,r). According to the first average value and the second average value, the multi-scale fuzzy entropy is calculated. The multi-scale fuzzy entropy is defined as the logarithmic difference of the average values ​​of the fuzzy membership of two adjacent embedded dimensions. The calculation formula of the multi-scale fuzzy entropy is: Where FuzzyEn is the multi-scale fuzzy entropy.

[0181] Then, the first variance is calculated based on the modal component to be evaluated and the average value of the modal component. The calculation formula of the first variance is: In the formula, σ is the first variance, x is the modal component to be evaluated, is the average value of the modal component, and n is the number of data in the modal component to be evaluated. The first skewness is calculated according to the modal component to be evaluated, the sample mean, the expected value operator and the first variance, where the calculation formula of the first skewness is: In the formula, S is, E is the expected value operator, μ is the sample mean, and σ is the first variance. According to the modal component to be evaluated, the sample mean, the expected value operator and the first variance, the first kurtosis is calculated, where the calculation formula of the first kurtosis is: Where K is the first kurtosis. The correlation coefficient is calculated based on the modal component to be evaluated, the preprocessed blasting vibration signal, the covariance calculation function and the standard deviation operator. The calculation formula of the correlation coefficient is: Where CC is the correlation coefficient, y is the preprocessed blasting vibration signal, cov is the covariance calculation function, S is the standard deviation operator, and x is the modal component to be evaluated.

[0182] In some embodiments, in step S106, the target modal component is calculated according to the component evaluation index, which may include but is not limited to the following steps:

[0183] Calculate the evaluation amount according to the component evaluation index;

[0184] Calculate the comprehensive evaluation value based on the total number of evaluation quantities and component evaluation indicators;

[0185] According to the preset evaluation threshold and the comprehensive evaluation value, the target modal component is selected from the modal components to be evaluated.

[0186] In some embodiments, the evaluation quantity can be calculated based on the component evaluation index. For example, the five component evaluation indexes of multi-scale fuzzy entropy, variance, skewness, kurtosis, and correlation coefficient can be set as P1, P2, P3, P4, and P5 respectively. j The rating is Q i1 , Q i2 , Q i3 , Q i4 , Q i5 , the calculation formula of the evaluation quantity is: In the formula, Q ij is the evaluation value of the i-th modal component to be evaluated with respect to the j-th component evaluation index, P ij is the jth component evaluation index of the i-th modal component to be evaluated, and n is the total number of modal components to be evaluated. Then, according to the evaluation quantity and the total number of component evaluation indexes, the comprehensive evaluation value of each modal component IMF to be evaluated is calculated, where the calculation formula of the comprehensive evaluation value is: Where M iis the comprehensive evaluation value, and R is the total number of component evaluation indicators. Finally, according to the preset evaluation threshold and the comprehensive evaluation value, the target modal component is selected from the modal components to be evaluated. For example, the preset evaluation threshold can be set to 85%. If the sum of the comprehensive evaluation values ​​of the first m modal components IMF to be evaluated accounts for more than 85% of the sum of the comprehensive evaluation values ​​of all modal components IMF to be evaluated, the first m modal components IMF to be evaluated are retained as the target modal components. The comprehensive evaluation value of the modal component to be evaluated is shown in the following example: Figure 5 As shown in the figure, the sum of the evaluation indicators of IMF1 to IMF3 accounts for 88.82% of the total evaluation indicators, which has exceeded 85%. Therefore, the first three IMFs are retained as target modal components.

[0187] In some embodiments, in step S107, constructing a target blasting vibration signal according to the target modal component may include but is not limited to the following steps:

[0188] Frequency modulate the target modal component to obtain a unit amplitude analytical signal;

[0189] Calculating a pseudo-distributed spectrum based on the unit amplitude analytic signal, the weighting function, the delay variable, the complex unit and the first frequency;

[0190] Calculate the original effective signal amplitude estimate based on the pseudo-distribution spectrum and the proportionality coefficient;

[0191] According to the scaling factor, the translation factor and the normalization factor, a scaling and translation wavelet function is constructed;

[0192] According to the scaling and translation wavelet function and the conjugate mother wavelet function, the original effective signal amplitude estimation is subjected to wavelet transformation to obtain the wavelet transformation coefficient;

[0193] According to the wavelet transform coefficients and sign functions, the soft threshold method is used to perform wavelet analysis on the target modal components to obtain approximate coefficients and detail coefficients.

[0194] According to the approximate coefficient and the detail coefficient, the target modal component is reconstructed by wavelet to obtain the reconstructed signal;

[0195] Calculate the matrix dimension based on the total number of modal components;

[0196] Construct the initial Hankel matrix according to the reconstructed signal and the matrix dimension;

[0197] According to the initial Hankel matrix, the weighted energy contribution rate is calculated;

[0198] According to the weighted energy contribution rate, the effective order of the singular value is determined;

[0199] According to the effective order of singular values, orthogonal matrix and non-negative diagonal matrix, the initial Hankel matrix is ​​subjected to singular value decomposition to obtain the target Hankel matrix;

[0200] According to the target Hankel matrix, the target blasting vibration signal is constructed.

[0201] In some embodiments, the retained m target modal components may be subjected to time-frequency peak filtering and then to wavelet threshold noise reduction, and then these components may be reconstructed. The reconstructed signal is constructed as a Hankel matrix, and then subjected to singular value decomposition, and then reconstructed to obtain the target blasting vibration signal. The target modal components may be frequency modulated first to obtain a unit amplitude analytical signal. Exemplarily, the essence of time-frequency peak filtering is instantaneous frequency estimation based on Wigner-Ville distribution. The noisy signal to be processed can be expressed as s(t)=x(t)+n(t), where x(t) is a valid signal and n(t) is additive random noise. The noisy signal is encoded to convert it into the form of an analytical signal; the noisy data s(t) is frequency modulated to obtain a unit amplitude analytical signal, where the calculation formula for the unit amplitude analytical signal is: In the formula, z s (t) is the unit amplitude analytical signal, j is the complex unit, μ is the frequency modulation index, t is time, and s(λ) is the instantaneous value of the target modal component. Then, according to the unit amplitude analytical signal, the weighting function, the delay variable, the complex unit and the first frequency, the pseudo Wigner-Ville distribution spectrum of the unit amplitude analytical signal is calculated, where the calculation formula of the pseudo distribution spectrum is: Where W z (t,f) is the pseudo-distributed spectrum of the unit amplitude analytical signal z(t) on the time-frequency plane (t,f), h(τ) is a weighting function, which is usually used to smooth the signal or reduce the influence of cross terms, τ is a delay variable, which represents the interval between two time points, is the unit amplitude analytical signal z(t) in time The value at for The conjugate of , f is the frequency.

[0202] Then, the original effective signal amplitude estimate is calculated based on the pseudo-distribution spectrum and the proportional coefficient. For example, the peak value of the PWVD distribution spectrum of the unit amplitude analytical signal can be obtained according to the maximum likelihood estimation principle as the instantaneous frequency estimate of the analytical signal to obtain the original effective signal amplitude estimate, wherein the calculation formula for the original effective signal amplitude estimate is: In the formula, is the original effective signal amplitude estimation, μ is the proportional coefficient. According to the scaling factor, translation factor and normalization factor, the scaling and translation wavelet function is constructed. In wavelet analysis, various basis functions are obtained by translation and scaling of the base wavelet. The calculation formula of the scaling and translation wavelet function is: In the formula, h a,b is the scaling and translation wavelet function, a is the scaling factor, b is the translation factor, is the normalization factor. According to the scaling and translation wavelet function and the conjugate mother wavelet function, the original effective signal amplitude is estimated by wavelet transform to obtain the wavelet transform coefficient, where the calculation formula of the wavelet transform coefficient is: Where, WT x (a,b) are wavelet transform coefficients, is the conjugate mother wavelet function after scaling and translation.

[0203] Then, based on the wavelet transform coefficients and sign functions, the soft threshold method is used to perform wavelet analysis on the target modal component to obtain approximate coefficients and detail coefficients. For example, wavelet analysis decomposes the signal into two parts: approximate and detail. Therefore, the original signal generates two signals through two complementary filters. The approximation represents the high-scale, low-frequency components of the signal; while the details represent the low-scale, high-frequency parts. The approximate part can be decomposed into the third layer of approximate and detail parts. By repeating this process, the signal can be decomposed into many low-resolution components (i.e., approximate coefficients). For each layer of detail from the first to the Nth layer, a threshold is selected for quantization processing to obtain the estimated wavelet coefficient (i.e., detail coefficient) Using the soft threshold processing method, the expression is: Wherein, D(x) is the detail coefficient after soft threshold processing, which indicates the high-frequency component retained after threshold processing, sign(x) is the sign function used to retain the positive and negative signs of the original coefficients, and λ is the set threshold used to determine whether the coefficient is regarded as noise. According to the approximate coefficients and detail coefficients, the target modal component is reconstructed by wavelet to obtain a reconstructed signal. For example, the signal can be reconstructed by wavelet based on the low-frequency coefficients (approximate coefficients) of the Nth layer of wavelet decomposition and the high-frequency estimates (detail coefficients) from the 1st layer to the Nth layer after quantization processing to obtain a reconstructed signal.

[0204] Then singular value decomposition can be performed to calculate the matrix dimension according to the total number of modal components, and the initial Hankel matrix can be constructed according to the reconstructed signal and the matrix dimension. The calculation formula of the matrix dimension is: n=N+1-m, where N is the total number of modal components, m is the number of matrix rows, and n is the number of matrix columns. The calculation formula for the initial Hankel matrix is: Where A is the initial Hankel matrix, x Nis the data in the reconstructed signal. Then, according to the initial Hankel matrix, the weighted energy contribution rate is calculated, and according to the weighted energy contribution rate, the effective order of the singular value is determined. Exemplarily, the method for determining the effective order of the singular value (i.e., the number of singular values ​​retained k) is the weighted energy contribution rate. The singular values ​​with a PCTE value (weighted energy contribution rate) greater than 0.1% are retained to determine the effective order of the singular value. The calculation formula for the weighted energy contribution rate is: In the formula, PCTE σi is the weighted energy contribution rate, σ i is the i-th singular value, a ij is an element in the matrix. Finally, according to the effective order of singular values, orthogonal matrix and non-negative diagonal matrix, the initial Hankel matrix is ​​decomposed into a singular value to obtain the target Hankel matrix, and the target blasting vibration signal is constructed according to the target Hankel matrix. For example, for any real matrix A∈R m×n , can be decomposed into singular values ​​according to the following formula, and the calculation formula of singular value decomposition is: Where U and V are orthogonal matrices, D is a non-negative diagonal matrix of rank r, and the expression of a non-negative diagonal matrix is: In the formula, S=diag(σ1,σ2,…,σ r ), σ r is the rth singular value, and D is a non-negative diagonal matrix. Retain the first K singular values ​​and make the other singular values ​​0, then a new non-negative diagonal matrix D can be formed. m×n , substituting it into the calculation formula of singular value decomposition and performing inverse transformation, we can get the new matrix It is then reconstructed into the target blasting vibration signal, which has undergone multiple noise reduction processes.

[0205] In some embodiments, in step S108, constructing a joint multi-domain feature vector according to the target blasting vibration signal may include but is not limited to the following steps:

[0206] Calculate the root mean square based on the target blasting vibration signal;

[0207] Calculating a second variance according to the target blasting vibration signal;

[0208] Calculating a second skewness according to the target blasting vibration signal, the sample mean and the second variance;

[0209] Calculating the second kurtosis according to the target blasting vibration signal, the sample mean and the second variance;

[0210] According to the root mean square, the second variance, the second skewness and the second kurtosis, a time domain feature vector is constructed;

[0211] Calculate the spectrum skewness based on the target blasting vibration signal, frequency variable, average frequency and frequency standard deviation;

[0212] Calculate the spectral kurtosis according to the target blasting vibration signal, frequency variable, average frequency and frequency standard deviation;

[0213] According to the spectral skewness and spectral kurtosis, the frequency domain feature vector is constructed;

[0214] Perform wavelet packet decomposition on the target blasting vibration signal to obtain multiple wavelet packet decomposition coefficients;

[0215] Reconstruct the wavelet packet decomposition coefficients to obtain the frequency band sub-signal;

[0216] Calculate signal energy according to the frequency band sub-signals;

[0217] Calculate the total energy based on multiple signal energies;

[0218] Calculate the energy probability based on the signal energy and the total energy;

[0219] Calculate energy entropy based on energy probability and total signal energy;

[0220] According to multiple energy entropies, a time-frequency domain feature vector is constructed;

[0221] A joint multi-domain feature vector is constructed according to the time domain feature vector, the frequency domain feature vector and the time-frequency domain feature vector.

[0222] In some embodiments, the time domain features of the target blasting vibration signal with noise removed may be first extracted. The root mean square is calculated based on the target blasting vibration signal, where the calculation formula of the root mean square is: In the formula, RMS is the root mean square, N is the number of data in the target blasting vibration signal, and x i is the target blasting vibration signal. According to the target blasting vibration signal, the second variance is calculated, wherein the calculation formula of the second variance is: In the formula, x var is the second variance, is the signal average value, and x is the target blasting vibration signal. The second skewness is calculated according to the target blasting vibration signal, the sample mean, and the second variance, where the calculation formula of the second skewness is: In the formula, x skew is the second skewness, μ is the sample mean, and E is the expected calculation function. The second kurtosis is calculated according to the target blasting vibration signal, the sample mean, and the second variance, where the calculation formula of the second kurtosis is: In the formula, x kurtAccording to the root mean square, the second variance, the second skewness and the second kurtosis, the time domain feature vector is constructed, where the expression of the time domain feature vector is: T1 = [t1, t2, ..., t4], where T1 is the time domain feature vector, t i The time domain features of the target blasting vibration signal to remove noise, i = 1, 2, ..., 4. It can be understood that the time domain features may include a root mean square, a second variance, a second skewness, and a second kurtosis.

[0223] Then, the frequency domain features of the target blasting vibration signal with noise removed are extracted. The spectral skewness is calculated based on the target blasting vibration signal, frequency variable, average frequency and frequency standard deviation, where the calculation formula of the spectral skewness is: In the formula, F skew is the spectral skewness, f i is the frequency, f ave is the average frequency, σ is the variance, S(f i ) is the frequency standard deviation, and k is the total number of frequencies. The spectrum kurtosis is calculated based on the target blasting vibration signal, frequency variable, average frequency and frequency standard deviation, where the calculation formula of the spectrum kurtosis is: In the formula, F kurt is the spectral kurtosis. Finally, the frequency domain feature vector is constructed based on the spectral skewness and spectral kurtosis. The expression of the frequency domain feature vector is: T2 = [t5, t6], where T2 is the frequency domain feature vector, t j The frequency domain features of the target blasting vibration signal for removing noise, j = 5, 6. It can be understood that the frequency domain features may include spectral skewness and spectral kurtosis.

[0224] Then, the target blasting vibration signal is decomposed by wavelet packet to obtain multiple wavelet packet decomposition coefficients, and the wavelet packet decomposition coefficients are reconstructed to obtain frequency band sub-signals. For example, the target blasting vibration signal can be decomposed by wavelet packet decomposition, and a three-layer wavelet packet decomposition is performed on it. The wavelet packet decomposition coefficients of 8 frequency bands are selected to reconstruct the wavelet decomposition coefficients, and the signals I in different frequency bands are extracted. 3,0 ~I 3,7 . Select a suitable wavelet basis function, and determine the number of decomposition layers to be 3. In each layer, the original signal or the output signal of the previous layer will be decomposed into two parts: the approximation signal (approximate component, also called the low-frequency part) and the detail signal (detail component, also called the high-frequency part). This process is implemented using a wavelet filter bank, including a low-pass filter and a high-pass filter, which are used to generate the approximation signal and the detail signal. After the signal passes through these two filters, each output is halved. After 3 layers of decomposition, the signal is decomposed into sub-signals of 8 frequency bands. According to the frequency band sub-signals, the signal energy is calculated. For example, after wavelet packet decomposition, 8 frequency band sub-signals I can be obtained. 3,0 ~I 3,7, assuming that the 8 frequency band sub-signals are (i=0, 1, 2…7), the calculation formula for the signal energy of each frequency band sub-signal is: In the formula, E i is the energy of the ith signal, T is the total number of data in the frequency band sub-signal, x i (t) is the ith frequency band sub-signal. The total energy is calculated based on the energy of multiple signals, where the total energy is calculated as follows: In the formula, E t is the total energy, and N is the total number of sub-signals in the frequency band. The energy probability is calculated based on the signal energy and the total energy, where the calculation formula of the energy probability is: In the formula, p i is the i-th energy probability. According to the energy probability and the total signal energy, the energy entropy is calculated, where the energy entropy calculation formula is: Where H is the energy entropy. Based on multiple energy entropies, a time-frequency domain feature vector is constructed. For example, the energy entropy of the eight frequency bands can be calculated, and the corresponding energy entropy is recorded as H 3,j (j=0,1,2,…,7), and use these 8 energy entropies to construct the time-frequency domain feature vector, where the expression of the time-frequency domain feature vector is: T3=[H 3,0 ,H 3,1 ,…,H 3,7 ], where T3 is the feature vector in the time-frequency domain, H 3,j is the energy entropy corresponding to the frequency band (j=0,1,2,…,7).

[0225] Finally, according to the time domain feature vector, frequency domain feature vector and time-frequency domain feature vector, the time domain feature vector, the frequency domain feature vector and the time-frequency domain feature vector can be connected in sequence to construct a joint multi-domain feature vector, where the expression of the joint multi-domain feature vector is: T4=[T1,T2,T3], where T4 is the joint multi-domain feature vector, T1 is the time domain feature vector, T2 is the frequency domain feature vector, and T3 is the time-frequency domain feature vector.

[0226] In some embodiments, in step S109, principal component analysis is performed on the joint multi-domain feature vector to obtain a joint multi-domain feature matrix, which may include but is not limited to the following steps:

[0227] Normalizing the joint multi-domain feature vector to obtain a standardized vector;

[0228] Calculate the covariance matrix based on the sample size and the standardized vector;

[0229] Perform eigenvalue decomposition on the covariance matrix to obtain the eigenvalues ​​and eigenvectors to be selected;

[0230] Calculate the cumulative variance contribution rate based on the eigenvalues ​​to be selected and the total number of eigenvalues;

[0231] According to the cumulative variance contribution rate, the target eigenvalue is selected from the eigenvalues ​​to be selected;

[0232] Selecting a target eigenvector from the eigenvectors to be selected according to the cumulative variance contribution rate and the target eigenvalue;

[0233] A joint multi-domain feature matrix is ​​constructed according to the target feature vector and the joint multi-domain feature vector.

[0234] In some embodiments, principal component analysis can be performed on the joint multi-domain feature vector. When the cumulative variance contribution rate of the first k principal components exceeds 85%, the first k principal components are retained to replace the original joint multi-domain feature vector to obtain a joint multi-domain feature matrix to achieve dimensionality reduction. The joint multi-domain feature vector can be first standardized to obtain a standardized vector, where the calculation formula of the standardized vector is: In the formula, X std is the standardized vector, X is the data in the joint multi-domain feature vector, μ is the mean, and σ is the standard deviation. Then, the covariance matrix is ​​calculated based on the number of samples and the standardized vector. The covariance matrix describes the linear correlation between each pair of features. Assume that the data in the joint multi-domain feature vector is a data set X with n samples and d features, and the covariance matrix C is a d×d square matrix. Each element C in the covariance matrix ij Represents the covariance between the i-th and j-th features. The calculation formula of the covariance matrix is: Where C is the covariance matrix and n is the number of samples. The covariance matrix is ​​subjected to eigenvalue decomposition to obtain the eigenvalues ​​to be selected and the eigenvectors to be selected. For example, the covariance matrix C is subjected to eigenvalue decomposition to obtain the eigenvalues ​​to be selected and the eigenvectors to be selected. The eigenvalues ​​represent the variance of the principal components, and the eigenvectors represent the directions of these principal components. The expression of the eigenvalue decomposition is: Cv i =λ i v i , where λ i is the eigenvalue to be selected, v i is the eigenvector to be selected corresponding to the eigenvalue to be selected. According to the eigenvalue to be selected and the total number of eigenvalues, the cumulative variance contribution rate is calculated, where the calculation formula of the cumulative variance contribution rate is: In the formula, ratio is the cumulative variance contribution rate, and p is the total number of eigenvalues. Finally, according to the cumulative variance contribution rate, the target eigenvalue is selected from the eigenvalues ​​to be selected, and according to the cumulative variance contribution rate and the target eigenvalue, the target eigenvector is selected from the eigenvector to be selected, and the joint multi-domain feature matrix is ​​constructed according to the target eigenvector and the joint multi-domain eigenvector. Exemplarily, the eigenvalues ​​to be selected are sorted from large to small, and the eigenvectors to be selected corresponding to the first k eigenvalues ​​to be selected whose cumulative variance contribution rate exceeds 85% are selected as new basis vectors. These new basis vectors are the principal components of the data in the new space, and they represent the variance that can be explained after the original data is projected to these directions. The original data is projected onto the selected k principal components to obtain the reduced-dimensional data (i.e., the target eigenvector). The matrix V is composed of the selected target eigenvectors. k (where each column is a target feature vector), then the calculation formula for the reduced-dimensional data Z is: Z = XV k , where Z is the joint multi-domain feature matrix, X is the target feature vector, V k is a matrix composed of target feature vectors, and the dimension of the joint multi-domain feature matrix after dimensionality reduction is n×k.

[0235] In some embodiments, in step S110, a characteristic information data set may be constructed, wherein the characteristic information data set includes a joint multi-domain characteristic matrix and a rock damage radius, and the rock damage radius is determined by comparing a preset change threshold with a wave velocity change obtained by an acoustic wave device. Exemplarily, after obtaining the characteristic information (i.e., the joint multi-domain characteristic matrix) of the blasting vibration signals at different measuring points after dimension reduction, a characteristic information data set is constructed. The constructed characteristic information data set includes the characteristic information after dimension reduction (i.e., the joint multi-domain characteristic matrix) and the rock damage radius corresponding thereto. The rock damage radius is determined by measuring the wave velocity change of the rock before and after the blasting by an acoustic wave meter. Each characteristic information of the blasting vibration signal obtained by monitoring after noise reduction, feature extraction, and dimension reduction will have a corresponding rock damage radius. Furthermore, in measuring the rock damage radius, an acoustic wave meter may be used to measure the wave velocity change of the rock before and after the blasting. When the rate of change of the acoustic wave velocity before and after the blasting is greater than 10%, it is considered that the rock has been damaged, and the area that meets this condition is called a rock damage zone; when the rate of change is less than 10%, the rock damage can be ignored. The horizontal damage range of the surrounding rock mass, i.e., the rock mass damage radius, is determined through the monitoring results of the sonic instrument.

[0236] In some embodiments, in step S111, the decision tree model can be optimized by using the sparrow search method according to the feature information data set to obtain a rock damage radius identification model. The feature information after dimension reduction (i.e., the joint multi-domain feature matrix) in the feature information data set is used as the feature variable, and the rock damage radius is used as the target variable. A rock radius damage identification model of the gradient boosting decision tree (GBDT) model optimized by the sparrow search algorithm is established, aiming to identify the rock damage radius corresponding to the feature information through the model. Some data in the feature information data set are selected, and 80% of them are randomly divided into a training set and 20% are divided into a test set. During the training process, the sparrow search algorithm is used to find the number of trees in the gradient boosting decision tree model, the maximum parameter of the tree, and the optimal parameter combination of the learning rate. Then, the gradient boosting decision tree model is trained again with this parameter combined with the training set to obtain the optimal rock radius damage identification model, and then the model performance is tested using the test set. During the training process of the gradient boosting decision tree model, a constant model F0(x) can be initialized first, Then calculate the residual r of the current model i (m) , r i (m) =y i -F m-1 (x i ), where r i (m) represents the residual of the i-th sample in the m-th iteration, y i is the true value of the ith sample, and is the predicted value of the model for the ith sample in the m-1th iteration. The residual is used as the new target value to train the next decision tree, minimize the loss of the residual, and correct the current model. In the formula, h m (x) is the model function of the mth tree; L(r i (m) ,h(x i )) is the loss function, which calculates the difference between the residual of each sample and the current model, h(x i ) is the predicted value of the i-th sample through the current decision tree model. Multiple iterations are performed, and each iteration updates the predicted value of the model by multiplying the predicted value of the new tree by the learning rate η. m (x) = F m-1 (x)+ηh m (x), where F m (x) is the prediction model of the mth round, η is the learning rate, h m (x) is the model function of the mth tree. After multiple rounds of iterations, the final model will get a gradually convergent result. The expression is as follows: In the formula, F M (x) is the model function after M rounds of iterations.

[0237] For example, the data in the feature information data set were screened, 80 sets of data were used to train the GBDT decision tree model, and 20 sets of data were used to test the model performance. The four performance indicators of the model, prediction accuracy (AR), mean square error (MSE), root mean square error (RMSE), and determination coefficient (R2), were 95.75%, 0.056, 0.237, and 0.938, respectively. The calculation formulas for the four performance indicators are as follows: In the formula, i is the i-th test sample, y i is the actual value, is the predicted value, is the mean of the predicted values ​​and n is the number of all test samples.

[0238] In some embodiments, in step S112, a rock damage radius identification model can be used to identify the rock damage radius. For example, a group of feature data is randomly selected from the feature information data set and input into the rock damage radius identification model. The rock blasting damage radius obtained by model identification is 6.89m, and the actual radius measured by the sonicator test is 6.46m. The model identification result has high accuracy and high identification efficiency.

[0239] In some embodiments, the overall process of analyzing the blasting vibration signal and obtaining the rock mass damage radius through model identification is as follows: Figure 6 As shown in the figure, the blasting vibration signals in three directions can be obtained first, the particle swarm algorithm can be used to calculate the optimal decomposition parameters, and variational mode decomposition can be performed to obtain multiple modal components to be evaluated. Then, multi-scale fuzzy entropy, correlation coefficient, variance, skewness and kurtosis are used for evaluation, and the k modal components ranked first are selected. Then, time-frequency peak filtering, wavelet threshold denoising and wavelet reconstruction are performed to obtain the denoised target blasting vibration signal. Then, according to the time domain characteristics, frequency domain characteristics and time-frequency domain characteristics corresponding to the target blasting vibration signal, a joint feature vector is constructed, and principal component analysis is performed to obtain the feature information after dimensionality reduction. The feature information after dimensionality reduction is corresponding to the rock damage radius obtained by the acoustic wave test to obtain the feature information data set, which is then divided into a training set and a test set. The sparrow search algorithm is used to optimize the gradient boosting decision tree model, and the rock loss radius recognition model is trained to obtain the rock loss radius recognition model. Finally, the rock loss radius recognition model is used for recognition to obtain the rock loss radius recognition result.

[0240] The beneficial effects of implementing the embodiments of the present invention include: the embodiments of the present invention first obtain the initial blasting vibration signal, pre-process the initial blasting vibration signal to obtain the pre-processed blasting vibration signal, then use the sparrow search method to perform parameter optimization processing based on the pre-processed blasting vibration signal, obtain the variational modal decomposition parameters, calculate the modal components to be evaluated, calculate the component evaluation index, calculate the target modal components, construct the target blasting vibration signal, and construct the joint multi-domain feature vector, then perform principal component analysis on the joint multi-domain feature vector to obtain the joint multi-domain feature matrix, and construct the feature information data set, and finally use the sparrow search method to optimize the decision tree model based on the feature information data set to obtain the rock damage radius identification model, and use the rock damage radius identification model to identify the rock damage radius, so that the rock damage radius can be identified by analyzing the blasting vibration signal, thereby improving efficiency. At the same time, this embodiment does not need to perform acoustic wave testing on the rock before and after blasting, simplifies the process of rock damage radius identification, reduces the consumption of manpower and material resources, and makes the determination of rock damage radius fast and convenient. This embodiment makes full use of the information contained in the blasting vibration signal, establishes the relationship between the blasting vibration signal and the rock damage radius, and can provide effective guidance for the design of blasting schemes.

[0241] like Figure 7 As shown, the embodiment of the present invention also provides a rock damage radius identification device based on blasting vibration signals, comprising:

[0242] The first module 801 is used to obtain an initial blasting vibration signal;

[0243] The second module 802 is used to pre-process the initial blasting vibration signal to obtain a pre-processed blasting vibration signal;

[0244] The third module 803 is used to perform parameter optimization processing using the sparrow search method according to the pre-processed blasting vibration signal to obtain variational modal decomposition parameters, where the variational modal decomposition parameters include the number of modes and a penalty factor;

[0245] The fourth module 804 is used to calculate the modal component to be evaluated according to the variational modal decomposition parameters;

[0246] The fifth module 805 is used to calculate the component evaluation index according to the modal component to be evaluated, and the component evaluation index includes multi-scale fuzzy entropy, first variance, first skewness, first kurtosis or correlation coefficient;

[0247] The sixth module 806 is used to calculate the target modal component according to the component evaluation index;

[0248] The seventh module 807 is used to construct a target blasting vibration signal according to the target modal component;

[0249] The eighth module 808 is used to construct a joint multi-domain feature vector according to the target blasting vibration signal;

[0250] The ninth module 809 is used to perform principal component analysis on the joint multi-domain feature vector to obtain a joint multi-domain feature matrix;

[0251] A tenth module 810 is used to construct a feature information data set, the feature information data set includes a joint multi-domain feature matrix and a rock damage radius, and the rock damage radius is obtained by an acoustic wave device;

[0252] The eleventh module 811 is used to optimize the decision tree model using the sparrow search method according to the characteristic information data set to obtain a rock damage radius identification model;

[0253] The twelfth module 812 is used to identify the rock damage radius using the rock damage radius identification model.

[0254] The contents of the above method embodiments are all applicable to the present device embodiments. The functions specifically implemented by the present device embodiments are the same as those of the above method embodiments, and the beneficial effects achieved are also the same as those achieved by the above method embodiments.

[0255] The preferred embodiments of the present invention are described above with reference to the accompanying drawings, but the scope of the rights of the present invention is not limited thereto. Any modification, equivalent substitution and improvement made by a person skilled in the art without departing from the scope and essence of the present invention should be within the scope of the rights of the present invention.

Claims

1. A rock damage radius identification method based on blasting vibration signals, characterized in that: The following steps are involved: Acquire initial blasting vibration signal; Preprocessing the initial blasting vibration signal to obtain a preprocessed blasting vibration signal; According to the pre-processed blasting vibration signal, a parameter optimization process is performed using a sparrow search method to obtain a variational modal decomposition parameter, wherein the variational modal decomposition parameter includes a mode number and a penalty factor; Calculating the modal components to be evaluated according to the variational modal decomposition parameters; Calculating component evaluation indicators according to the modal component to be evaluated, wherein the component evaluation indicators include multi-scale fuzzy entropy, first variance, first skewness, first kurtosis or correlation coefficient; Calculating the target modal component according to the component evaluation index; constructing a target blasting vibration signal according to the target modal component; constructing a joint multi-domain feature vector according to the target blasting vibration signal; Performing principal component analysis on the joint multi-domain feature vector to obtain a joint multi-domain feature matrix; Constructing a characteristic information data set, the characteristic information data set including the joint multi-domain characteristic matrix and a rock damage radius, the rock damage radius being determined by comparing a preset change threshold with a wave velocity change obtained by an acoustic wave device; According to the characteristic information data set, the decision tree model is optimized by using the sparrow search method to obtain a rock damage radius identification model; The rock mass damage radius identification model is used to identify the rock mass damage radius.

2. The method according to claim 1, characterized in that The preprocessing of the initial blasting vibration signal to obtain a preprocessed blasting vibration signal includes: Get the time node of blasting implementation; Comparing the blasting implementation time node with the initial blasting vibration signal in time to obtain a blasting vibration signal period; The pre-processed blasting vibration signal is extracted from the initial blasting vibration signal according to the blasting vibration signal period.

3. The method according to claim 1, characterized in that The method of performing parameter optimization processing by using the sparrow search optimization method according to the preprocessed blasting vibration signal to obtain variational mode decomposition parameters includes: Setting parameters of the sparrow search method, the parameters of the sparrow search method include the value range of the number of modes, the value range of the penalty factor, the population size, the maximum number of iterations, the forager ratio, the guard ratio, the step length control parameter, the early warning movement amplitude coefficient, the disturbance coefficient or the random direction quantity; Generating a search population according to the parameters of the sparrow search method, wherein the types of individuals in the search population include foragers, followers and guards; Constructing a sparrow search objective function, wherein the optimization target of the sparrow search objective function is to minimize the mean square error between the sum of the modal components to be optimized and the preprocessed blasting vibration signal; Constructing a forager position update function according to the search population, the population convergence rate constant, the maximum number of iterations and the direction vector; Constructing a follower position update function according to the search population and the random disturbance value; Constructing a sentinel position update function according to the search population and the moving speed constant; The variational mode decomposition parameters are calculated according to the sparrow search objective function, the forager position update function, the follower position update function and the sentinel position update function.

4. The method according to claim 1, characterized in that The step of calculating the modal component to be evaluated according to the variational modal decomposition parameter comprises: Constructing a variational mode decomposition objective function according to the number of modes, the modal components to be optimized, the center frequency and the time differential operator; Constructing conditional constraints according to the modal components to be optimized; Constructing a Lagrangian function according to the variational mode decomposition objective function, the conditional constraints, the Lagrangian multiplier and the penalty factor; Constructing a modal component update function to be optimized according to the Lagrangian function, the frequency variable and the penalty factor; Constructing a center frequency update function according to the modal component to be optimized and the frequency variable; Constructing a Lagrange multiplier update function according to the modal components to be optimized and the conditional constraints; The modal component to be evaluated is calculated using an alternating direction multiplier method according to the modal component update function to be optimized, the center frequency update function and the Lagrange multiplier update function.

5. The method according to claim 1, characterized in that The step of calculating the component evaluation index according to the modal component to be evaluated includes: Calculating a coarsening sequence according to the modal components to be evaluated and the scale factors; constructing an embedding vector according to a first preset embedding vector dimension, the coarsening sequence and the scale factor; According to the embedding vector, calculating the embedding vector distance; Calculating fuzzy membership according to the embedding vector distance; Calculating a first average value according to the fuzzy membership degree; Calculating a second average value according to a second preset embedding vector dimension, where the second preset embedding vector dimension is the first preset embedding vector dimension plus 1; Calculating the multi-scale fuzzy entropy according to the first average value and the second average value; Calculating the first variance according to the modal component to be evaluated and the modal component average value; Calculating the first skewness according to the modal component to be evaluated, the sample mean, the expected value operator and the first variance; Calculating the first kurtosis according to the modal component to be evaluated, the sample mean, the expected value operator and the first variance; The correlation coefficient is calculated according to the modal component to be evaluated, the preprocessed blasting vibration signal, the covariance calculation function and the standard deviation operator.

6. The method according to claim 1, characterized in that The step of calculating the target modal component according to the component evaluation index comprises: Calculating the evaluation amount according to the component evaluation index; Calculate the comprehensive evaluation value according to the evaluation amount and the total number of component evaluation indicators; The target modal component is selected from the modal components to be evaluated according to a preset evaluation threshold and the comprehensive evaluation value.

7. The method according to claim 1, characterized in that The step of constructing a target blasting vibration signal according to the target modal component comprises: Frequency modulating the target modal component to obtain a unit amplitude analytical signal; Calculating a pseudo-distributed spectrum according to the unit amplitude analytical signal, the weighting function, the delay variable, the complex unit and the first frequency; Calculating an original effective signal amplitude estimate based on the pseudo-distribution spectrum and the proportionality coefficient; According to the scaling factor, the translation factor and the normalization factor, a scaling and translation wavelet function is constructed; According to the scaling and translation wavelet function and the conjugate mother wavelet function, the original effective signal amplitude estimation is subjected to wavelet transformation to obtain wavelet transformation coefficients; According to the wavelet transform coefficients and the sign function, a soft threshold method is used to perform wavelet analysis on the target modal component to obtain approximate coefficients and detail coefficients; Performing wavelet reconstruction on the target modal component according to the approximate coefficient and the detail coefficient to obtain a reconstructed signal; Calculate the matrix dimension based on the total number of modal components; constructing an initial Hankel matrix according to the reconstructed signal and the matrix dimension; Calculating a weighted energy contribution rate according to the initial Hankel matrix; Determining the effective order of singular values ​​according to the weighted energy contribution rate; According to the effective order of the singular values, the orthogonal matrix and the non-negative diagonal matrix, the initial Hankel matrix is ​​subjected to singular value decomposition to obtain a target Hankel matrix; The target blasting vibration signal is constructed according to the target Hankel matrix.

8. The method according to claim 1, characterized in that The step of constructing a joint multi-domain feature vector according to the target blasting vibration signal includes: Calculating a root mean square according to the target blasting vibration signal; Calculating a second variance according to the target blasting vibration signal; Calculating a second skewness according to the target blasting vibration signal, the sample mean and the second variance; Calculating a second kurtosis according to the target blasting vibration signal, the sample mean and the second variance; constructing a time domain feature vector according to the root mean square, the second variance, the second skewness and the second kurtosis; Calculating spectral skewness according to the target blasting vibration signal, frequency variable, average frequency and frequency standard deviation; Calculating spectral kurtosis according to the target blasting vibration signal, the frequency variable, the average frequency and the frequency standard deviation; Constructing a frequency domain feature vector according to the spectral skewness and the spectral kurtosis; Performing wavelet packet decomposition on the target blasting vibration signal to obtain a plurality of wavelet packet decomposition coefficients; Reconstructing the wavelet packet decomposition coefficients to obtain frequency band sub-signals; Calculating signal energy according to the frequency band sub-signals; Calculating total energy based on multiple signal energies; Calculating energy probability according to the signal energy and the total energy; Calculating energy entropy based on the energy probability and the total signal energy; Constructing a time-frequency domain feature vector according to the plurality of energy entropies; The joint multi-domain feature vector is constructed according to the time domain feature vector, the frequency domain feature vector and the time-frequency domain feature vector.

9. The method according to claim 1, characterized in that: The principal component analysis is performed on the joint multi-domain feature vector to obtain a joint multi-domain feature matrix, including: Normalizing the joint multi-domain feature vector to obtain a normalized vector; Calculate the covariance matrix based on the sample size and the standardized vector; Performing eigenvalue decomposition on the covariance matrix to obtain eigenvalues ​​to be selected and eigenvectors to be selected; Calculate the cumulative variance contribution rate according to the eigenvalues ​​to be selected and the total number of eigenvalues; Selecting a target eigenvalue from the eigenvalues ​​to be selected according to the cumulative variance contribution rate; Selecting a target eigenvector from the eigenvectors to be selected according to the cumulative variance contribution rate and the target eigenvalue; The joint multi-domain feature matrix is ​​constructed according to the target feature vector and the joint multi-domain feature vector.

10. A rock damage radius identification device based on blasting vibration signals, characterized in that: include: The first module is used to obtain the initial blasting vibration signal; The second module is used to preprocess the initial blasting vibration signal to obtain a preprocessed blasting vibration signal; The third module is used to perform parameter optimization processing using the sparrow search method according to the pre-processed blasting vibration signal to obtain variational modal decomposition parameters, wherein the variational modal decomposition parameters include the number of modes and a penalty factor; The fourth module is used to calculate the modal component to be evaluated according to the variational modal decomposition parameters; A fifth module is used to calculate a component evaluation index according to the modal component to be evaluated, wherein the component evaluation index includes a multi-scale fuzzy entropy, a first variance, a first skewness, a first kurtosis or a correlation coefficient; A sixth module is used to calculate the target modal component according to the component evaluation index; A seventh module is used to construct a target blasting vibration signal according to the target modal component; An eighth module is used to construct a joint multi-domain feature vector according to the target blasting vibration signal; A ninth module is used to perform principal component analysis on the joint multi-domain feature vector to obtain a joint multi-domain feature matrix; A tenth module is used to construct a feature information data set, wherein the feature information data set includes the joint multi-domain feature matrix and a rock damage radius, and the rock damage radius is obtained by an acoustic wave device; The eleventh module is used to optimize the decision tree model using the sparrow search method according to the characteristic information data set to obtain a rock damage radius identification model; The twelfth module is used to identify the rock damage radius using the rock damage radius identification model.

Images