Tunnel blasting construction safety risk intelligent early warning method, system, device and medium

By using a distributed sensor network and an intelligent early warning system, the problem of real-time monitoring of the multi-physics coupling process in tunnel blasting construction has been solved, enabling accurate early warning of safety risks in harsh environments and improving the safety and efficiency of tunnel blasting construction.

CN122328211APending Publication Date: 2026-07-03CHINA RAILWAY NO 2 ENG GROUP CO LTD
View PDF 0 Cites 0 Cited by

Patent Information

Authority / Receiving Office
CN · China
Patent Type
Applications(China)
Current Assignee / Owner
CHINA RAILWAY NO 2 ENG GROUP CO LTD
Filing Date
2026-04-01
Publication Date
2026-07-03

Smart Images

  • Figure CN122328211A_ABST
    Figure CN122328211A_ABST
Patent Text Reader

Abstract

This application relates to an intelligent early warning method, system, equipment, and medium for safety risks in tunnel blasting construction. The method includes: based on a trigger signal, hard-triggering all intelligent acquisition nodes to simultaneously begin data acquisition, and recording the local hardware timestamp of each intelligent acquisition node as the zero moment of blasting; based on the acceleration time-series data of the entire blasting process, combined with a crystal oscillator acceleration sensitivity model, calculating the dynamic clock drift error, and performing time-axis correction on the acquired data to obtain a dynamically corrected time axis and local time-based physical field data; registering each local time-based physical field data onto the dynamically corrected time axis to obtain a multi-physics dataset, and extracting evolutionary and coupling features to obtain a blasting cycle feature vector; inputting the blasting cycle feature vector into a long short-term memory neural network early warning model to obtain the safety risk level and handling recommendations. This method enables synchronous data acquisition in tunnels without GPS signals and in environments with strong electromagnetic interference during blasting.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] This invention belongs to the field of blasting risk early warning technology, and in particular relates to an intelligent early warning method, system, equipment and medium for safety risks in tunnel blasting construction. Background Technology

[0002] With the development of tunnel blasting excavation technology, safety risk assessment technology for tunnel blasting construction has emerged. Based on the static geological survey data carried out before blasting, and combined with the professional experience of on-site technical personnel, subjective judgments are made to generate risk assessment conclusions.

[0003] In traditional technologies, construction monitoring uses conventional methods such as ground-penetrating radar for advanced detection, post-blasting rock convergence measurement, and vibration velocity meter single-point recording. These methods focus only on geological surveys before blasting and data collection of lag effects after blasting to carry out risk management.

[0004] However, the aforementioned methods, with their assessment models contradicting the dynamic nature of blasting operations, cannot cover the entire process of multi-physics field coupling evolution during blasting, making it difficult to accurately predict chain risks such as gas outbursts, surrounding rock instability, and structural damage. Conventional monitoring methods cannot capture the transient physical field evolution characteristics of blasting at the second to millisecond level, directly leading to a serious lag in safety warnings. Hidden dangers such as gas exceeding limits, surrounding rock deformation, and structural damage are often only discovered after they have become apparent, significantly increasing the probability of accidents. Problems such as satellite timing failure, clock distortion caused by strong impacts and strong electromagnetic interference, and asynchronous amplification of data from multiple types of sensors become prominent in extreme environments during tunnel blasting. Existing time synchronization technologies cannot adapt to harsh construction environments, hindering the implementation of real-time quantitative monitoring and accurate early warning of multi-physics field coupling processes. Summary of the Invention

[0005] Therefore, it is necessary to provide an intelligent early warning method, system, equipment, and medium for tunnel blasting construction safety risks that can overcome the lack of GPS signals in tunnels and resist strong electromagnetic interference at the moment of blasting, in order to address the above-mentioned technical problems.

[0006] Firstly, this application provides an intelligent early warning method for safety risks in tunnel blasting construction, including:

[0007] The frequency drift coefficient and phase offset of the clock of each slave node relative to the master node in the distributed sensor network of the tunnel to be blasted area are determined by static clock drift calibration, and the clock synchronization baseline is obtained; the distributed sensor network includes multiple intelligent acquisition nodes;

[0008] Based on the non-contact electromagnetic induction trigger on the detonation network, in response to the trigger signal generated at the moment of detonation, all intelligent acquisition nodes are hard-triggered to start data acquisition simultaneously, obtain the acquired data, and record the local hardware timestamp of each intelligent acquisition node as the zero moment of detonation; after the field programmable gate array corresponding to each intelligent acquisition node detects the rising edge of the interrupt pin caused by the trigger signal, the hard-triggered node freezes all analog-to-digital conversion channels of the corresponding intelligent acquisition node and latches the local hardware counter value as the local hardware timestamp corresponding to the detonation trigger.

[0009] Based on the acceleration time series data of the entire blasting process recorded by the accelerometers built into each intelligent acquisition node, combined with the pre-trained crystal oscillator acceleration sensitive model, the dynamic clock drift error caused by impact vibration is calculated, and the time axis of the acquired data is corrected by the dynamic clock drift error to obtain the dynamic correction time axis of each intelligent acquisition node and the local time base physical field data with the blasting zero time corresponding to each intelligent acquisition node as the origin.

[0010] Based on the clock synchronization baseline, the local time base physical field data are uniformly registered to the dynamic correction time axis of the master node to obtain a multi-physics dataset.

[0011] The evolution characteristics of the gas emission physical field, vibration propagation physical field, and stress redistribution physical field are extracted from the multiphysics dataset, and the coupling characteristics of the physical field pairs are combined to obtain the blasting cycle feature vector.

[0012] The feature vector of the blasting cycle is input into a pre-trained long short-term memory neural network early warning model to obtain the safety risk level of the current blasting cycle and the corresponding handling suggestions.

[0013] In one embodiment, the clock frequency drift coefficient and phase offset of each slave node relative to the master node in the distributed sensor network of the tunnel to be blasted area are determined by static clock drift calibration to obtain a clock synchronization baseline, including:

[0014] Based on the periodic broadcast synchronization messages of the master node, the frequency drift coefficient and phase offset of each slave node's clock relative to the master node are estimated online using a constant gain recursive least squares algorithm for the timestamp pairs of multiple consecutive synchronization messages recorded by each slave node. The recursive formula for the constant gain recursive least squares algorithm is as follows: , , ,in, For the first The clock parameter vector obtained by the recursive estimation. This is an estimate of the frequency drift coefficient. This is the estimated phase offset value. For the first The timestamp recorded by the slave node in the secondary synchronization message. , For the first The master node's clock in the secondary synchronization message. Here is the gain matrix. Let be the error covariance matrix. Forgetting factor;

[0015] Based on the frequency drift coefficient and phase offset, a linear correction model for the local time of the slave node relative to the clock of the master node is calculated; the expression for the linear correction model is: ;

[0016] Based on the linear correction model of all slave nodes, a hierarchical clock synchronization tree with the master node as the root node is constructed to obtain the clock synchronization baseline.

[0017] In one embodiment, based on the acceleration time-series data of the entire blasting process recorded by the accelerometers built into each intelligent acquisition node, and combined with a pre-trained crystal oscillator acceleration-sensitive model, the dynamic clock drift error caused by impact vibration is calculated. The acquired data is then corrected on the time axis using the dynamic clock drift error, resulting in the dynamically corrected time axis of each intelligent acquisition node and the local time-based physical field data of each intelligent acquisition node with the blasting zero moment corresponding to the intelligent acquisition node as the origin, including:

[0018] Based on the crystal oscillator acceleration sensitivity model, the relative frequency deviation of the crystal oscillator is calculated using acceleration time-series data from the entire blasting process; the expression for the crystal oscillator acceleration sensitivity model is: ,in, for The relative frequency deviation of the crystal oscillator at any given time. This is the acceleration sensitivity coefficient matrix. For memory effect kernel function, for Acceleration time-series data for the entire blasting process at any given moment;

[0019] The dynamic clock drift error caused by impact vibration is calculated based on the relative frequency deviation of the crystal oscillator; the formula for calculating the dynamic clock drift error is as follows: ;

[0020] The time axis of the acquired data is nonlinearly corrected using dynamic clock drift error to obtain a dynamically corrected time axis and local time-based physical field data with the zero time of the explosion as the origin; the expression for nonlinear correction is: ,in, These are the time axis coordinates of each sampling point on the time axis. To dynamically correct the time axis coordinates, This refers to dynamic clock drift error.

[0021] In one embodiment, based on a clock synchronization baseline, the local time-based physical field data are uniformly registered onto the dynamic correction time axis of the master node to obtain a multi-physics dataset, including:

[0022] Using the dynamic correction time axis of the master node as a reference, the local time-based physical field data of each slave node are regarded as an asynchronous observation sequence;

[0023] For each physical quantity data in the asynchronous observation sequence, optimal state estimation is performed using a Kalman filter to obtain the filtered values ​​of each physical quantity;

[0024] Based on the uniform sampling time point sequence on the dynamic correction time axis of the master node, the filtered values ​​of all physical quantities corresponding to each slave node are calculated by Lagrange interpolation to obtain the physical quantity estimates of the corresponding uniform sampling time point sequence.

[0025] By aligning all physical quantity estimates from all intelligent acquisition nodes with the dynamic correction time axis of the master node, a multiphysics dataset is obtained.

[0026] In one embodiment, the evolution features of the gas emission physical field, vibration propagation physical field, and stress redistribution physical field are extracted from the multiphysics dataset, and combined with the coupling features of the physical field pairs to obtain the blasting cycle feature vector, including:

[0027] Based on the time series data of gas concentration after blasting from the multiphysics dataset, the initial rise slope, peak concentration, time to reach the peak and decay time of gas emission are calculated to obtain the evolution characteristics of the gas emission physical field.

[0028] Based on the blasting vibration signal from the multiphysics dataset, the arrival times of longitudinal and transverse waves are identified, and the peak velocity, dominant frequency distribution, energy duration, and response spectrum intensity are calculated to obtain the evolution characteristics of the vibration propagation physical field.

[0029] Based on the stress-strain data of the surrounding rock from the multiphysics dataset, the adjustment magnitude, adjustment rate and number of abrupt change points of stress redistribution are calculated to obtain the evolution characteristics of the stress redistribution physical field.

[0030] Based on the time delay between the vibration energy accumulation curve and the starting point of gas emission growth, the gas-vibration coupling characteristics are obtained;

[0031] Based on the correlation coefficient between the stress unloading rate and the maximum peak value of the vibration wave, the stress-vibration coupling characteristics are obtained;

[0032] Based on the evolution characteristics of the vibration propagation physical field, the evolution characteristics of the gas emission physical field, the stress redistribution physical field, the gas-vibration coupling characteristics, and the stress-vibration coupling characteristics, a blasting cycle feature vector is constructed.

[0033] Secondly, this application also provides an intelligent early warning system for safety risks in tunnel blasting construction, including:

[0034] The pre-blasting static calibration module is used to determine the frequency drift coefficient and phase offset of the clock of each slave node relative to the master node in the distributed sensor network of the tunnel to be blasted area through static clock drift calibration, and to obtain the clock synchronization baseline; the distributed sensor network includes multiple intelligent acquisition nodes.

[0035] The blasting synchronization module is used to respond to the trigger signal generated at the moment of blasting initiation by the non-contact electromagnetic induction trigger on the initiation network. It hard-triggers all intelligent acquisition nodes to start data acquisition simultaneously, obtain the acquired data, and record the local hardware timestamp of each intelligent acquisition node as the zero moment of blasting. After the field programmable gate array of each intelligent acquisition node detects the rising edge of the interrupt pin caused by the trigger signal, it freezes all analog-to-digital conversion channels of the corresponding intelligent acquisition node and latches the local hardware counter value as the local hardware timestamp corresponding to the blasting trigger.

[0036] The post-blast drift compensation module is used to calculate the dynamic clock drift error caused by impact vibration based on the acceleration time series data of the entire blasting process recorded by the accelerometers built into each intelligent acquisition node, combined with the pre-trained crystal oscillator acceleration sensitive model. The module then corrects the time axis of the acquired data through the dynamic clock drift error, thereby obtaining the dynamic correction time axis of each intelligent acquisition node and the local time base physical field data of each intelligent acquisition node with the blasting zero time corresponding to the intelligent acquisition node as the origin.

[0037] The configuration module is used to uniformly register the local time-base physical field data to the dynamic correction time axis of the master node based on the clock synchronization baseline, so as to obtain a multi-physics dataset.

[0038] The feature module is used to extract the evolution features of the gas emission physical field, vibration propagation physical field, and stress redistribution physical field from the multi-physics dataset, and combine them with the coupling features of the physical field pairs to obtain the blasting cycle feature vector.

[0039] The early warning module is used to input the feature vector of the blasting cycle into a pre-trained long short-term memory neural network early warning model to obtain the safety risk level of the current blasting cycle and the corresponding handling suggestions.

[0040] Thirdly, this application also provides a computer device, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the steps of any of the above-described intelligent early warning methods for safety risks in tunnel blasting construction.

[0041] Fourthly, this application also provides a computer-readable storage medium having a computer program stored thereon, wherein the computer program, when executed by a processor, implements the steps of any of the above-described intelligent early warning methods for safety risks in tunnel blasting construction.

[0042] The aforementioned intelligent early warning method, system, equipment, and medium for tunnel blasting construction safety risks establish a baseline for the time relationship between nodes by constructing a distributed sensor network and performing static clock drift calibration. Using a hardware interrupt signal generated by a non-contact electromagnetic induction trigger on the detonation network, all intelligent acquisition nodes' field-programmable gate arrays are hard-triggered to simultaneously latch their local timestamps as a unified blasting zero moment and begin data acquisition. Combining the acceleration time-series data of the entire blasting process recorded by accelerometers with a pre-trained crystal oscillator acceleration-sensitive model, the dynamic clock drift error caused by impact vibration is calculated and compensated, thereby obtaining the dynamic correction time axis and local time-based physical field data of each node with the blasting zero moment as the origin. Based on the clock synchronization baseline obtained by static calibration, the local time-based physical field data of all nodes are uniformly registered to the dynamic correction time axis of the master node, forming a strictly synchronized multi-physics dataset. The evolution characteristics of gas outburst, vibration propagation, and stress redistribution, as well as the coupling characteristics between them, are extracted from this dataset to form a blasting cycle feature vector. This vector is input into a pre-trained long short-term memory neural network early warning model, which outputs the safety risk level of the current blasting cycle and corresponding handling suggestions. By combining hardware triggering and dynamic compensation, the challenge of millisecond-level synchronous acquisition of multiple sensors in tunnels without GPS signals and under conditions of strong electromagnetic interference during blasting was solved. This enabled precise capture and quantitative analysis of the coupled evolution of gas outburst, vibration propagation, and stress redistribution during blasting. As a result, tunnel blasting safety risk assessment was transformed from a lagging model relying on static geological surveys and manual experience to a dynamic intelligent early warning model based on real-time synchronous data. This significantly improved the accuracy and timeliness of risk identification and provided actionable early warning results with clear disposal suggestions for on-site construction. Attached Figure Description

[0043] To more clearly illustrate the technical solutions in the embodiments or related technologies of this application, the accompanying drawings used in the description of the embodiments or related technologies will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0044] Figure 1 This is a flowchart illustrating the intelligent early warning method for safety risks in tunnel blasting construction according to the present invention.

[0045] Figure 2 This is a structural diagram of the intelligent early warning system for safety risks in tunnel blasting construction according to the present invention. Detailed Implementation

[0046] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this application.

[0047] In one embodiment, such as Figure 1 As shown, an intelligent early warning method for safety risks in tunnel blasting construction is provided. This embodiment illustrates the method's application to a terminal, but it is understood that the method can also be applied to a server, or to a system including both a terminal and a server, and implemented through interaction between the terminal and the server. In this embodiment, the method includes the following steps:

[0048] S101. The frequency drift coefficient and phase offset of the clock of each slave node relative to the master node in the distributed sensor network of the tunnel to be blasted area are determined by static clock drift calibration, and the clock synchronization baseline is obtained; the distributed sensor network includes multiple intelligent acquisition nodes.

[0049] Indicatively, the distributed sensor network is an architecture with strong resistance to electromagnetic interference. Each intelligent acquisition node is deployed in the stable rock mass section behind the tunnel face to be blasted, in the pre-set blast hole area around the tunnel face, and at the adjacent structural protection target. Each intelligent acquisition node is a miniature measurement and control unit that integrates multiple types of sensor arrays, a high-stability local clock module, and an FPGA embedded processing core. The multiple types of sensor arrays include gas monitoring sensors, triaxial vibration monitoring detectors, and surrounding rock stress and strain sensors. The high-stability local clock module uses a temperature-controlled crystal oscillator as the clock source and is independently electromagnetically shielded. All nodes are connected through a wired fiber optic network to avoid the influence of electromagnetic interference.

[0050] Specifically, during the calibration process, the node in the stable rock mass area with the least blasting disturbance is selected as the master clock node. Initial time synchronization of the master node is completed at the tunnel entrance or in a signal-enabled area using a portable high-precision time synchronization device. Furthermore, the network-wide constant-gain recursive least squares time synchronization algorithm is initiated. For example, the recursive formula corresponding to the constant-gain recursive least squares time synchronization algorithm is: , , ,in Let be the vector of clock parameters to be estimated. To observe the time from the node, For the observation matrix, Here is the gain matrix. Let be the error covariance matrix. Forgetting factor, This is the frequency drift coefficient. Phase offset. Through iterative algorithm calculation, the static frequency drift coefficient and phase offset of each slave node relative to the master node are obtained to form a clock synchronization baseline.

[0051] S102. The non-contact electromagnetic induction trigger based on the detonation network responds to the trigger signal generated at the moment of detonation, hard-triggers all intelligent acquisition nodes to start data acquisition simultaneously, obtains the acquired data, and records the local hardware timestamp of each intelligent acquisition node as the zero moment of detonation; after the field programmable gate array corresponding to each intelligent acquisition node detects the rising edge of the interrupt pin caused by the trigger signal, it freezes all analog-to-digital conversion channels of the corresponding intelligent acquisition node and latches the local hardware counter value as the local hardware timestamp corresponding to the detonation trigger.

[0052] Furthermore, a non-contact electromagnetic induction trigger is installed on the main line of the blasting initiation network. When the detonator discharges and generates an initiation current flowing through the busbar, it instantly generates a strong pulsed magnetic field. The trigger's induction coil is excited by the magnetic field, generating a steep rising edge trigger signal. This steep rising edge trigger signal is directly connected to the FPGA (Field Programmable Gate Array) hardware interrupt pin of each intelligent acquisition node, and the FPGA is programmed to set this hardware interrupt as the highest response priority. When the FPGA detects a rising edge signal that reaches the threshold, it freezes all analog-to-digital conversion acquisition channels of this node within a nanosecond time period, and simultaneously latches the real-time value of the local hardware counter, recording this value as the local hardware timestamp for each node, thus marking the zero moment of blasting. The hard triggering process is executed at the hardware level, without the need for software protocol stack forwarding processing, completely avoiding the interference of software processing delay, network transmission delay, and electromagnetic interference on the transmission of synchronization messages, enabling all nodes to start data acquisition at the same physical moment and ensuring the initial synchronization of the data acquired by each node.

[0053] S103. Based on the acceleration time series data of the entire blasting process recorded by the accelerometers built into each intelligent acquisition node, and combined with the pre-trained crystal oscillator acceleration sensitive model, calculate the dynamic clock drift error caused by impact vibration, and correct the time axis of the acquired data through the dynamic clock drift error to obtain the dynamic correction time axis of each intelligent acquisition node and the local time base physical field data of each intelligent acquisition node with the blasting zero time corresponding to the intelligent acquisition node as the origin.

[0054] Specifically, each intelligent acquisition node incorporates a high-sensitivity triaxial micromechanical accelerometer, rigidly mounted on the same circuit board as a temperature-controlled crystal oscillator. It acquires triaxial acceleration time-series data throughout the entire blasting process at a sampling rate far exceeding the dominant frequency of seismic waves. This triaxial acceleration time-series data directly reflects the disturbance state of the nodes caused by the blasting impact vibration. The pre-trained crystal oscillator acceleration sensitivity model was obtained through laboratory vibration table testing and calibration. ,in This refers to the relative frequency deviation of the crystal oscillator. The acceleration sensitivity coefficient matrix characterizes the instantaneous frequency drift caused by transient impact. It is the instantaneous acceleration vector. The kernel function represents the memory effect, characterizing the hysteresis recovery properties of the crystal oscillator after vibration. Based on the crystal oscillator acceleration sensitivity model and collected acceleration time-series data, the dynamic frequency deviation of the crystal oscillator at each moment during the blast impact period is calculated through integral operations. The cumulative phase time deviation is obtained, and based on this, the original acquired data is nonlinearly resampled and corrected to eliminate the time axis distortion caused by mechanical impact. Finally, the dynamic corrected time axis of each node is generated, as well as the local time-based gas, vibration, stress and strain physical field data with the zero time of its own explosion as the origin.

[0055] S104. Based on the clock synchronization baseline, the local time base physical field data are uniformly registered to the dynamic correction time axis of the master node to obtain a multi-physics dataset.

[0056] Optionally, after dynamic clock drift error correction, the local data streams of each intelligent acquisition node have achieved accurate and uniform acquisition. However, due to individual crystal oscillator differences and compensation residuals, there is still a slight sub-millisecond deviation in the time axis of each node relative to the master clock. Time registration across the entire network can be completed by relying on the clock synchronization baseline. The registration process employs an asynchronous data fusion preprocessing method combining Kalman filtering and Lagrange interpolation. First, the Kalman filter state equation is constructed. With observation equation ,in For the actual physical quantity to be estimated, Here is the state transition matrix. For the observation matrix, For process noise, To mitigate noise, Kalman filtering is used to smooth and delay the data from each slave node, eliminating measurement errors and inherent sensor response biases. Furthermore, Lagrange interpolation is employed with the uniform sampling time of the master clock node as the target time. Complete data interpolation, where The estimated physical quantity of the target time of the master clock. For the physical quantity of the neighboring filtered data points, These represent the target time and the data point sampling time, respectively. All physical field data after filtering and interpolation from the slave nodes are uniformly calculated to the dynamic correction time axis of the master node, integrating multi-dimensional data such as gas concentration, vibration velocity, and surrounding rock stress and strain to form a standardized multi-physics dataset with equally spaced sampling and strict temporal synchronization.

[0057] S105. Extract the evolution features of the gas emission physical field, vibration propagation physical field, and stress redistribution physical field from the multi-physics dataset, and combine them with the coupling features of the physical field pairs to obtain the blasting cycle feature vector.

[0058] For example, for a standardized multiphysics dataset, the evolution features of individual physics fields are first extracted class by class. Specifically, for the vibration propagation physics field, the arrival times of P-waves and S-waves are identified, and peak vibration velocity, dominant frequency distribution, energy duration, and response spectrum intensity features are extracted. For the gas emission physics field, the initial concentration rise slope, peak concentration attainment time, peak concentration value, and concentration decay time constant features are extracted. For the stress redistribution physics field, stress adjustment amplitude, adjustment rate, and stress abrupt change state features are extracted. Furthermore, the focus is on extracting multiphysics coupling features. By quantitatively analyzing the intrinsic correlation characteristics between physics fields, including the time delay between vibration energy accumulation and gas emission growth, the correlation between the surrounding rock stress unloading rate and the peak value of the vibration wave, and the temporal correlation between stress redistribution and rock mass fracture, this approach can overcome the limitations of judging by single physical quantity thresholds and accurately reflect the risk nature of the transient multiphysics coupling evolution during blasting. The extracted single-physics field evolution features and multi-physics field coupling features are regularized and fused, redundant features are eliminated, and feature dimensions are unified to construct a blasting cycle feature vector that represents the full-dimensional safety state of a single blasting cycle.

[0059] S106. Input the blasting cycle feature vector into the pre-trained long short-term memory neural network early warning model to obtain the safety risk level of the current blasting cycle and the corresponding handling suggestions.

[0060] Optionally, a long short-term memory neural network early warning model is used. This model selects synchronous multiphysics monitoring data from historical blasting shifts and corresponding risk assessment results as training samples. The prediction error is calculated using the cross-entropy loss function, and the model parameters are iteratively optimized using an adaptive moment estimation algorithm until the model converges. In the real-time early warning phase, the feature vector sequence of the current blasting cycle is input into the trained model. The model, through internal neuron operations, outputs the probability of each risk level: no risk, low risk, medium risk, and high risk. The level with the highest probability is selected as the final safety risk level for the current blasting cycle. Simultaneously, corresponding engineering response suggestions are output: low risk corresponds to normal construction and routine monitoring and control; medium risk corresponds to optimizing blasting parameters and increasing monitoring frequency; and high risk corresponds to emergency response measures such as personnel evacuation, emergency ventilation, and suspension of construction. This transforms blasting construction risk analysis from post-event analysis to precise real-time early warning, comprehensively improving the efficiency of safety management in tunnel blasting construction.

[0061] In the aforementioned intelligent early warning method for safety risks in tunnel blasting construction, static clock drift calibration is performed on the distributed sensor network in the tunnel to be blasted area to determine the clock parameters of each slave node relative to the master node and obtain the clock synchronization baseline. Relying on the non-contact electromagnetic induction triggers on the initiation network to respond to the blasting initiation signal, all intelligent acquisition nodes are hard-triggered to synchronously collect data and record the local hardware timestamp as the zero moment of blasting. Combining the acceleration time-series data of the entire blasting process collected by the built-in accelerometers of each intelligent acquisition node with the pre-trained crystal oscillator acceleration sensitivity model, the dynamic clock drift error is calculated and the time axis correction of the collected data is completed, obtaining the dynamic correction time axis of each node and the local time-based physical field data. Based on the clock synchronization baseline, various local time-based physical field data are integrated... A standardized multiphysics dataset is formed by dynamically correcting the timeline from the master node. Multiphysics evolution features and physics pair coupling features are extracted from this dataset, and a blasting cycle feature vector is constructed. This feature vector is then input into a pre-trained long short-term memory neural network early warning model, which outputs the safety risk level of the current blasting cycle and corresponding handling suggestions. This solves the problem of simultaneous multi-sensor data acquisition in tunnels without GPS, under strong electromagnetic interference, and with high impact environments. It abandons the traditional static, experience-based assessment model, achieving dynamic quantitative analysis and intelligent early warning of tunnel blasting construction safety risks. This overcomes the limitations of single physical quantity threshold judgment, significantly improving the accuracy, timeliness, and comprehensiveness of risk warnings, and providing reliable support for the safety management of tunnel blasting construction.

[0062] In one embodiment, the clock frequency drift coefficient and phase offset of each slave node relative to the master node in the distributed sensor network of the tunnel to be blasted area are determined by static clock drift calibration to obtain a clock synchronization baseline, including:

[0063] S11. Based on the periodic broadcast synchronization messages of the master node, for the timestamp pairs of multiple consecutive synchronization messages recorded by each slave node, the constant gain recursive least squares algorithm is used to estimate the frequency drift coefficient and phase offset of the clock of each slave node relative to the master node online; the recursive formula of the constant gain recursive least squares algorithm is as follows: , , ,in, For the first The clock parameter vector obtained by the recursive estimation. This is an estimate of the frequency drift coefficient. This is the estimated phase offset value. For the first The timestamp recorded by the slave node in the secondary synchronization message. , For the first The master node's clock in the secondary synchronization message. Here is the gain matrix. Let be the error covariance matrix. It is a forgetting factor.

[0064] As an illustration, the master node is configured to broadcast synchronization messages to the entire network at fixed intervals during the static calibration phase. Each slave node, upon receiving the message, accurately records its local timestamp at the corresponding reception time, forming multiple master-slave node timestamp pairs. Specifically, Indicates the first The clock parameter vector obtained by the recursive estimation is composed of the frequency drift coefficient estimate and the phase offset estimate. For the first The frequency drift coefficient estimate obtained in the next iteration represents the degree of frequency deviation between the slave node clock and the master node. For the first The phase offset estimate obtained in the next iteration represents the initial time deviation of the slave node clock relative to the master node; For the first During the transmission of the secondary synchronization message, the local receiving timestamp recorded by the node; Let be the observation vector, given by the first... The clock value of the master node in the secondary synchronization message is composed of a combination of the clock value of the master node and the constant 1; For the first The master node clock value carried in the secondary synchronization message; For the first The gain matrix of each iteration controls the magnitude of the iteration correction. For the first The error covariance matrix of the next iteration reflects the error level of the parameter estimation; This is a forgetting factor used to mitigate the influence of historical data on the current estimation results, improving the algorithm's ability to track slow clock drift. Through continuous iterative calculations, it gradually approximates the true frequency drift coefficient and phase offset of the slave node's clock relative to the master node, achieving accurate online estimation of clock parameters.

[0065] S12. Based on the frequency drift coefficient and phase offset, calculate the linear correction model of the slave node's local time relative to the master node's clock; the expression for the linear correction model is: .

[0066] Furthermore, based on the estimated stable frequency drift coefficient and phase offset, a linear mapping relationship between the master and slave node clocks is constructed. This linear correction model can intuitively reflect the quantization correlation between the slave node's local time and the master node's reference time. This represents the time value of the local clock on the slave node, which is the original time data to be corrected. The final estimated value of the frequency drift coefficient after iterative convergence is the constant correction coefficient. This indicates the time value of the master node's reference clock, which serves as a unified time reference for the entire network. The final estimated value of the phase shift after iterative convergence is the constant correction intercept.

[0067] S13. Based on the linear correction model of all slave nodes, construct a hierarchical clock synchronization tree with the master node as the root node to obtain the clock synchronization baseline.

[0068] Specifically, based on the deployment location and transmission topology of the tunnel distributed sensor network, each slave node is hierarchically divided. Slave nodes located at the core transmission nodes, where signal transmission is stable, are classified as mid-level nodes, while those deployed at the edge and directly collecting physical field data are classified as end-level nodes. The master node serves as the top-level root node of the entire clock synchronization system, and a hierarchical clock synchronization tree is built level by level. During the construction process, the linear correction model corresponding to each slave node is embedded into the node control logic of the corresponding level, clarifying the clock transmission and correction relationships between the master node and mid-level nodes, and between mid-level nodes and end-level nodes, forming a hierarchical clock correction system covering all intelligent acquisition nodes in the entire network. The hierarchical clock synchronization tree is the final clock synchronization baseline, clarifying the clock correction rules and hierarchical relationships of all nodes in the network. This ensures the stability of clock synchronization in a static environment and adapts to the clock error compensation requirements during subsequent dynamic blasting processes, providing a unified benchmark framework for the time registration of multi-source sensor data and avoiding synchronization failures caused by clock hierarchy confusion.

[0069] In one embodiment, based on the acceleration time-series data of the entire blasting process recorded by the accelerometers built into each intelligent acquisition node, and combined with a pre-trained crystal oscillator acceleration-sensitive model, the dynamic clock drift error caused by impact vibration is calculated. The acquired data is then corrected on the time axis using the dynamic clock drift error, resulting in the dynamically corrected time axis of each intelligent acquisition node and the local time-based physical field data of each intelligent acquisition node with the blasting zero moment corresponding to the intelligent acquisition node as the origin, including:

[0070] S21. Based on the crystal oscillator acceleration sensitivity model, the relative frequency deviation of the crystal oscillator is calculated using acceleration time-series data from the entire blasting process; the expression for the crystal oscillator acceleration sensitivity model is: ,in, for The relative frequency deviation of the crystal oscillator at any given time. This is the acceleration sensitivity coefficient matrix. For memory effect kernel function, for Acceleration time-series data for the entire blasting process at any given moment.

[0071] For example, for The relative frequency deviation of the crystal oscillator at a given time represents the degree to which the actual oscillation frequency of the crystal oscillator deviates from its nominal frequency at that time. The acceleration sensitivity coefficient matrix is ​​determined by the crystal oscillator's factory characteristics and laboratory vibration table calibration, and is used to quantify the instantaneous disturbance effect of transient impact acceleration on the crystal oscillator frequency. This is the memory effect kernel function, reflecting the hysteresis recovery characteristics of the crystal oscillator's frequency deviation after being subjected to vibrational excitation. This refers to the time lag. The acceleration time series data for the entire blasting process at time t is given. , , These correspond to the instantaneous acceleration values ​​collected by the accelerometer in the x, y, and z axes, respectively, fully covering the vibration effects of blasting impact on the intelligent acquisition node in the three spatial dimensions.

[0072] S22. The dynamic clock drift error caused by impact vibration is calculated based on the relative frequency deviation of the crystal oscillator; the formula for calculating the dynamic clock drift error is as follows: .

[0073] Specifically, for The dynamic clock drift error caused by impact vibration is characterized from the zero moment of blasting. arrive At any given moment, the total time offset formed by the cumulative frequency deviation of the crystal oscillator due to impact vibration; for The relative frequency deviation of the crystal oscillator at any given moment is the core integrand in the integration operation; the lower limit of integration, 0, corresponds to the zero moment of the explosion, and the upper limit... For the target moment when the drift error is to be calculated, the integration process fully covers the entire process from the occurrence of the blast impact to... The cumulative process of frequency deviation throughout the entire time period.

[0074] S23. The time axis of the acquired data is nonlinearly corrected using dynamic clock drift error to obtain a dynamically corrected time axis and local time-based physical field data with the zero moment of the explosion as the origin; the expression for nonlinear correction is: ,in, These are the time axis coordinates of each sampling point on the time axis. To dynamically correct the time axis coordinates, This refers to dynamic clock drift error.

[0075] The time axis of the raw acquired data exhibits nonlinear distortion due to the dynamic drift of the crystal oscillator. The standard for nonlinear correction is... Inverse compensation is performed on the original time coordinates to achieve nonlinear correction of the time axis. The time axis coordinates of each sampling point in the original data acquisition time axis are the sampling times recorded by the local clock of the intelligent acquisition node, without eliminating clock drift caused by impact vibration. The time axis coordinates for dynamic correction are the time values ​​that reflect the actual evolution of the physical field after correction. For the corresponding original time point Dynamic clock drift error.

[0076] In one embodiment, based on a clock synchronization baseline, the local time-based physical field data are uniformly registered onto the dynamic correction time axis of the master node to obtain a multi-physics dataset, including:

[0077] S31. Using the dynamic correction time axis of the master node as a reference, the local time-based physical field data of each slave node are regarded as an asynchronous observation sequence.

[0078] Indicatively, the master node's dynamic calibration timeline is a unified time base formed after static clock drift calibration and dynamic clock drift error correction, reflecting the true evolution of the physical field. The sampling times on this timeline have a unified time scale across the entire network. Although each slave node completes dynamic calibration of its own local timeline, due to factors such as sensor hardware response characteristics and differences in data transmission delays, the sampling times of its local time-based physical field data cannot be perfectly aligned with the sampling times of the master node's dynamic calibration timeline. Therefore, the physical field data of slave nodes exhibiting temporal asynchrony are defined as asynchronous observation sequences.

[0079] S32. For each physical quantity data in the asynchronous observation sequence, perform optimal state estimation using a Kalman filter to obtain the filtered value of each physical quantity.

[0080] Specifically, Kalman filtering, based on the system's state equation and observation equation, and combining the statistical characteristics of observation noise and process noise, can optimally estimate physical quantity data in asynchronous observation sequences, eliminating the interference of random noise on the data. The core equations of the Kalman filter include the state equation and the observation equation. The state equation is as follows: The observation equation is ,in, for The true state value of the physical quantity at any given time is the core parameter to be estimated; The state transition matrix represents the physical quantity transitioning from... Time's up The evolutionary pattern of time; It is process noise, which follows a zero-mean Gaussian distribution and reflects random perturbations in the evolution of the physical field; for Physical quantity observations in a time-asynchronous observation sequence; For the observation matrix, establish a linear mapping relationship between the true state of physical quantities and the observed values; To observe noise, a zero-mean Gaussian distribution is used to reflect measurement errors during sensor acquisition. Through a Kalman filter prediction-update iterative process, optimal state estimation is performed on the physical quantity data of the asynchronous observation sequence of each slave node, ultimately obtaining the filtered values ​​of each physical quantity, effectively reducing data bias caused by noise interference.

[0081] S33. Based on the uniform sampling time point sequence on the dynamic correction time axis of the master node, the filtered values ​​of all physical quantities corresponding to each slave node are calculated by Lagrange interpolation to obtain the physical quantity estimates of the corresponding uniform sampling time point sequence.

[0082] Furthermore, the dynamic correction time axis of the master node is divided into a uniform sampling time point sequence according to a preset sampling frequency, serving as the target time node for data registration across the entire network. Lagrange interpolation can solve for the estimated physical quantity value at any target time point based on discrete filtered physical quantity values, i.e. ,in, Uniform sampling time points of the target on the master node dynamic correction time axis The corresponding estimated physical quantity; The number of neighboring filtered values ​​used in the interpolation calculation; The j-th neighboring filtered value is the value obtained from the filtered values ​​of the node physical quantities at the target time point. The closest discrete value; This represents the original sampling time corresponding to the j-th neighboring filtered value; This represents the original sampling time corresponding to all neighboring filtered values ​​except for the j-th value. Through Lagrange interpolation, the discrete physical quantity filtered values ​​of each slave node are mapped to the uniform sampling time points of the master node's dynamic correction time axis, achieving time-series matching between slave node data and the master node's time axis.

[0083] S34. Align all physical quantity estimates of all smart acquisition nodes with the dynamic correction time axis of the master node to obtain a multiphysics dataset.

[0084] After completing the Lagrange interpolation of the filtered physical quantity values ​​of each slave node, the physical quantity data such as gas emission, vibration propagation, and stress redistribution from all intelligent acquisition nodes have been mapped to the same set of uniform sampling time points on the dynamic correction time axis of the master node. Based on the sampling time order of this time axis, the estimated values ​​of various physical quantities from all nodes are integrated and aligned to form a structured dataset containing time dimension, node spatial dimension, and physical field type dimension—that is, a multiphysics dataset.

[0085] In one embodiment, the evolution features of the gas emission physical field, vibration propagation physical field, and stress redistribution physical field are extracted from the multiphysics dataset, and combined with the coupling features of the physical field pairs to obtain the blasting cycle feature vector, including:

[0086] S41. Based on the time series data of gas concentration after blasting from the multi-physics dataset, calculate the initial rise slope, peak concentration, time to reach the peak and decay time of gas emission, and obtain the evolution characteristics of the gas emission physical field.

[0087] Indicatively, relying on the centralized and normalized gas concentration time series data of the multiphysics field dataset, quantitative calculations are carried out on the full-cycle variation law of gas emission after blasting. The initial rise slope is obtained by differential calculation of the initial rise segment of gas concentration, which characterizes the gas emission rate in the early stage of rock mass fracture penetration after blasting, directly reflecting the initial intensity of gas migration; the peak concentration is the maximum value in the gas concentration time series data, which reflects the upper limit of the total gas emission induced by blasting; the time to reach the peak is the time from the zero moment of blasting to the gas concentration rising to the peak, which characterizes the response speed of gas emission and fracture penetration efficiency; the decay time is the time it takes for the gas concentration to fall from the peak to the stable baseline value, which reflects the dissipation rate and ventilation effect after gas emission.

[0088] S42. Based on the blasting vibration signal from the multi-physics dataset, identify the arrival times of longitudinal and transverse waves, and calculate the peak velocity, dominant frequency distribution, energy duration, and response spectrum intensity to obtain the evolution characteristics of the vibration propagation physical field.

[0089] Optionally, by detecting abrupt amplitude changes and analyzing waveform characteristics of vibration signals, the initial arrival times of longitudinal and transverse waves can be accurately identified, clarifying the propagation sequence and laws of blasting stress waves. Furthermore, multi-dimensional parameter extraction is carried out based on the time-domain and frequency-domain characteristics of the vibration signals. The peak velocity of vibration is the extreme value of the vibration time-domain signal, characterizing the core impact intensity of the blasting vibration; the dominant frequency distribution is obtained through spectrum transformation calculation, reflecting the energy concentration frequency band of the blasting vibration, which can determine the impact characteristics of the vibration on the surrounding rock and adjacent structures; the energy duration is the duration for which the vibration signal maintains an effective amplitude, reflecting the period of action of the blasting vibration; the response spectrum intensity is calculated in combination with dynamic response characteristics, characterizing the dynamic damage potential of blasting vibration on tunnel support structures and surface buildings, comprehensively depicting the evolution characteristics of the physical field of vibration propagation.

[0090] S43. Based on the stress-strain data of the surrounding rock from the multi-physics dataset, calculate the adjustment magnitude, adjustment rate and number of abrupt change points of stress redistribution to obtain the evolution characteristics of the stress redistribution physical field.

[0091] For example, regarding the dynamic process of stress redistribution in the surrounding rock after blasting, quantitative analysis is conducted based on the normalized stress-strain time series data of the surrounding rock. The adjustment range of stress redistribution is the absolute value of the difference between the stress in the surrounding rock before and after blasting, reflecting the fluctuation range of the stress in the surrounding rock under blasting disturbance; the adjustment rate is the ratio of the stress adjustment range to the corresponding adjustment time, characterizing the reconstruction speed of the stress field in the surrounding rock and reflecting the stress adjustment capability of the surrounding rock; the number of abrupt change points is obtained by the slope abrupt change detection algorithm of the stress-strain curve to determine whether there is an abnormal sudden change in the stress field. The more abrupt change points there are, the more unstable the stress state of the surrounding rock is, and the higher the probability of inducing disasters such as rock bursts and collapses.

[0092] S44. Based on the time delay between the vibration energy accumulation curve and the starting point of gas emission growth, the gas-vibration coupling characteristics are obtained.

[0093] Furthermore, energy integration is first performed on the blasting vibration signal to generate a vibration energy accumulation curve, clarifying the accumulation law of vibration energy over time and locking in the key time point when the vibration energy reaches the critical threshold. By detecting the inflection point of the gas emission time series data, the initial time point when the gas emission begins to increase is located. The time difference between the above two key time points is calculated. This time delay is the gas-vibration coupling characteristic, which quantifies the lag time of blasting vibration inducing rock mass fracture and opening gas migration channels. It intuitively reflects the intrinsic relationship between vibration and gas emission. The shorter the time delay, the stronger the coupling between rock mass fracture and gas emission, and the higher the risk of gas outburst and exceeding limits.

[0094] S45. Based on the correlation coefficient between the stress unloading rate and the maximum peak value of the vibration wave, the stress-vibration coupling characteristics are obtained.

[0095] Optionally, the Pearson correlation coefficient method can be used to perform a quantitative correlation analysis on the two sets of data, i.e. ,in, This is the correlation coefficient between the stress unloading rate and the maximum peak value of the vibration wave. For the first Sample values ​​of group stress unloading rate This represents the sample mean of the stress unloading rate. The maximum peak value of the i-th vibration wave sample is given. The mean value of the sample of the maximum peak value of the vibration wave. This represents the total number of valid samples. The larger the absolute value of the correlation coefficient, the stronger the coupling relationship between stress unloading and vibration propagation. Under a high positive correlation coefficient, rapid stress unloading in the surrounding rock is accompanied by high-intensity vibration, which is a typical precursor to rockburst and surrounding rock instability.

[0096] S46. Based on the evolution characteristics of the vibration propagation physical field, the evolution characteristics of the gas emission physical field, the stress redistribution physical field, the gas-vibration coupling characteristics, and the stress-vibration coupling characteristics, the blasting cycle feature vector is constructed.

[0097] Specifically, the extracted single physical field evolution features and multi-field coupling features are first regularized, redundant features are removed and the feature dimensions are unified to eliminate the interference of dimension differences on subsequent model calculations. Furthermore, according to the preset feature sequence order, the various effective features are arranged in an orderly manner to finally construct a standardized blasting cycle feature vector.

[0098] It should be understood that although the steps in the flowcharts of the embodiments described above are shown sequentially according to the arrows, these steps are not necessarily executed in the order indicated by the arrows. Unless explicitly stated herein, there is no strict order restriction on the execution of these steps, and they can be executed in other orders. Moreover, at least some steps in the flowcharts of the embodiments described above may include multiple steps or multiple stages. These steps or stages are not necessarily completed at the same time, but can be executed at different times. The execution order of these steps or stages is not necessarily sequential, but can be performed alternately or in turn with other steps or at least some of the steps or stages of other steps.

[0099] Based on the same inventive concept, this application also provides a tunnel blasting construction safety risk intelligent early warning system for implementing the above-mentioned intelligent early warning method for tunnel blasting construction safety risks. The solution provided by this system is similar to the solution described in the above method. Therefore, the specific limitations of one or more tunnel blasting construction safety risk intelligent early warning system embodiments provided below can be found in the limitations of the tunnel blasting construction safety risk intelligent early warning method described above, and will not be repeated here.

[0100] In one exemplary embodiment, such as Figure 2 As shown, an intelligent early warning system for safety risks in tunnel blasting construction is provided, including:

[0101] The pre-blasting static calibration module 201 is used to determine the frequency drift coefficient and phase offset of the clock of each slave node relative to the master node in the distributed sensor network of the tunnel to be blasted area through static clock drift calibration, and to obtain the clock synchronization baseline; the distributed sensor network includes multiple intelligent acquisition nodes.

[0102] The blasting synchronization module 202 is used to respond to the trigger signal generated at the moment of blasting initiation by the non-contact electromagnetic induction trigger on the initiation network, hard triggering all intelligent acquisition nodes to start data acquisition simultaneously, obtain the acquired data, and record the local hardware timestamp of each intelligent acquisition node as the zero moment of blasting; after the field programmable gate array corresponding to each intelligent acquisition node detects the rising edge of the interrupt pin caused by the trigger signal, it freezes all analog-to-digital conversion channels of the corresponding intelligent acquisition node and latches the local hardware counter value as the local hardware timestamp corresponding to the blasting trigger;

[0103] The post-blast drift compensation module 203 is used to calculate the dynamic clock drift error caused by impact vibration based on the acceleration time series data of the entire blasting process recorded by the accelerometers built into each intelligent acquisition node, combined with the pre-trained crystal oscillator acceleration sensitive model, and to correct the time axis of the acquired data through the dynamic clock drift error, so as to obtain the dynamic correction time axis of each intelligent acquisition node and the local time base physical field data with the blasting zero time corresponding to the intelligent acquisition node as the origin;

[0104] Configuration module 204 is used to uniformly register each local time base physical field data to the dynamic correction time axis of the master node based on the clock synchronization baseline, so as to obtain a multi-physics dataset.

[0105] Feature module 205 is used to extract the evolution features of gas emission physical field, vibration propagation physical field and stress redistribution physical field from the multi-physics dataset, and combine them with the coupling features of the physical field pairs to obtain the blasting cycle feature vector.

[0106] The early warning module 206 is used to input the blasting cycle feature vector into a pre-trained long short-term memory neural network early warning model to obtain the safety risk level of the current blasting cycle and the corresponding handling suggestions.

[0107] In one embodiment, the pre-blasting static calibration module 201 is further configured to:

[0108] Based on the periodic broadcast synchronization messages of the master node, the frequency drift coefficient and phase offset of each slave node's clock relative to the master node are estimated online using a constant gain recursive least squares algorithm for the timestamp pairs of multiple consecutive synchronization messages recorded by each slave node. The recursive formula for the constant gain recursive least squares algorithm is as follows: , , ,in, For the first The clock parameter vector obtained by the recursive estimation. This is an estimate of the frequency drift coefficient. This is the estimated phase offset value. For the first The timestamp recorded by the slave node in the secondary synchronization message. , For the first The master node's clock in the secondary synchronization message. Here is the gain matrix. Let be the error covariance matrix. Forgetting factor;

[0109] Based on the frequency drift coefficient and phase offset, a linear correction model for the local time of the slave node relative to the clock of the master node is calculated; the expression for the linear correction model is: ;

[0110] Based on the linear correction model of all slave nodes, a hierarchical clock synchronization tree with the master node as the root node is constructed to obtain the clock synchronization baseline.

[0111] In one embodiment, the post-blast drift compensation module 203 is further configured to:

[0112] Based on the crystal oscillator acceleration sensitivity model, the relative frequency deviation of the crystal oscillator is calculated using acceleration time-series data from the entire blasting process; the expression for the crystal oscillator acceleration sensitivity model is: ,in, for The relative frequency deviation of the crystal oscillator at any given time. This is the acceleration sensitivity coefficient matrix. For memory effect kernel function, for Acceleration time-series data for the entire blasting process at any given moment;

[0113] The dynamic clock drift error caused by impact vibration is calculated based on the relative frequency deviation of the crystal oscillator; the formula for calculating the dynamic clock drift error is as follows: ;

[0114] The time axis of the acquired data is nonlinearly corrected using dynamic clock drift error to obtain a dynamically corrected time axis and local time-based physical field data with the zero time of the explosion as the origin; the expression for nonlinear correction is: ,in, These are the time axis coordinates of each sampling point on the time axis. To dynamically correct the time axis coordinates, This refers to dynamic clock drift error.

[0115] In one embodiment, the configuration module 204 is further configured to:

[0116] Using the dynamic correction time axis of the master node as a reference, the local time-based physical field data of each slave node are regarded as an asynchronous observation sequence;

[0117] For each physical quantity data in the asynchronous observation sequence, optimal state estimation is performed using a Kalman filter to obtain the filtered values ​​of each physical quantity;

[0118] Based on the uniform sampling time point sequence on the dynamic correction time axis of the master node, the filtered values ​​of all physical quantities corresponding to each slave node are calculated by Lagrange interpolation to obtain the physical quantity estimates of the corresponding uniform sampling time point sequence.

[0119] By aligning all physical quantity estimates from all intelligent acquisition nodes with the dynamic correction time axis of the master node, a multiphysics dataset is obtained.

[0120] In one embodiment, the feature module 205 is further configured to:

[0121] Based on the time series data of gas concentration after blasting from the multiphysics dataset, the initial rise slope, peak concentration, time to reach the peak and decay time of gas emission are calculated to obtain the evolution characteristics of the gas emission physical field.

[0122] Based on the blasting vibration signal from the multiphysics dataset, the arrival times of longitudinal and transverse waves are identified, and the peak velocity, dominant frequency distribution, energy duration, and response spectrum intensity are calculated to obtain the evolution characteristics of the vibration propagation physical field.

[0123] Based on the stress-strain data of the surrounding rock from the multiphysics dataset, the adjustment magnitude, adjustment rate and number of abrupt change points of stress redistribution are calculated to obtain the evolution characteristics of the stress redistribution physical field.

[0124] Based on the time delay between the vibration energy accumulation curve and the starting point of gas emission growth, the gas-vibration coupling characteristics are obtained;

[0125] Based on the correlation coefficient between the stress unloading rate and the maximum peak value of the vibration wave, the stress-vibration coupling characteristics are obtained;

[0126] Based on the evolution characteristics of the vibration propagation physical field, the evolution characteristics of the gas emission physical field, the stress redistribution physical field, the gas-vibration coupling characteristics, and the stress-vibration coupling characteristics, a blasting cycle feature vector is constructed.

[0127] In one embodiment, a computer device is provided, including a memory and a processor, the memory storing a computer program, the processor executing the computer program to implement the steps in the above method embodiments.

[0128] In one embodiment, a computer-readable storage medium is provided having a computer program stored thereon, which, when executed by a processor, implements the steps in the above method embodiments.

[0129] For the device embodiments, since they basically correspond to the method embodiments, the relevant parts can be referred to in the description of the method embodiments. The device embodiments described above are merely illustrative. The components described as separate parts may or may not be physically separate, and the components shown as units may or may not be physical units, that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this disclosure according to actual needs. Those skilled in the art can understand and implement this without creative effort.

[0130] The above-described embodiments are merely illustrative of several implementation methods of the embodiments of this application, and their descriptions are relatively specific and detailed. However, they should not be construed as limiting the scope of the patent application. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the embodiments of this application, and these modifications and improvements all fall within the protection scope of the embodiments of this application.

Claims

1. A method for intelligent early warning of safety risks in tunnel blasting construction, characterized in that, The method includes: The frequency drift coefficient and phase offset of the clock of each slave node relative to the master node in the distributed sensor network of the tunnel to be blasted area are determined by static clock drift calibration, and the clock synchronization baseline is obtained; the distributed sensor network includes multiple intelligent acquisition nodes. Based on the non-contact electromagnetic induction trigger on the detonation network, in response to the trigger signal generated at the moment of detonation, all the aforementioned intelligent acquisition nodes are hard-triggered to start data acquisition simultaneously, obtain the acquired data, and record the local hardware timestamp of each of the aforementioned intelligent acquisition nodes as the zero moment of detonation; after the field programmable gate array corresponding to each of the aforementioned intelligent acquisition nodes detects the rising edge of the interrupt pin caused by the trigger signal, it freezes all analog-to-digital conversion channels of the corresponding intelligent acquisition node and latches the local hardware counter value as the local hardware timestamp corresponding to the detonation trigger. Based on the acceleration time series data of the entire blasting process recorded by the accelerometers built into each of the intelligent acquisition nodes, combined with the pre-trained crystal oscillator acceleration sensitive model, the dynamic clock drift error caused by impact vibration is calculated, and the acquisition data is corrected on the time axis through the dynamic clock drift error, so as to obtain the dynamic correction time axis of each intelligent acquisition node and the local time base physical field data of each intelligent acquisition node with the blasting zero time corresponding to the intelligent acquisition node as the origin. Based on the clock synchronization baseline, each of the local time-based physical field data is uniformly registered to the dynamic correction time axis of the master node to obtain a multi-physics dataset. The evolution characteristics of the gas emission physical field, vibration propagation physical field, and stress redistribution physical field are extracted from the multiphysics dataset, and the coupling characteristics of the physical field pairs are combined to obtain the blasting cycle feature vector. The blasting cycle feature vector is input into a pre-trained long short-term memory neural network early warning model to obtain the safety risk level of the current blasting cycle and the corresponding handling suggestions.

2. The method according to claim 1, characterized in that, The process of determining the clock frequency drift coefficient and phase offset of each slave node relative to the master node in the distributed sensor network of the tunnel to be blasted area through static clock drift calibration, and obtaining the clock synchronization baseline, includes: Based on the periodic broadcast synchronization messages of the master node, the frequency drift coefficient and phase offset of the clock of each slave node relative to the master node are estimated online using a constant gain recursive least squares algorithm for the timestamp pairs of consecutive synchronization messages recorded by each slave node; the recursive formula of the constant gain recursive least squares algorithm is as follows: , , ,in, For the first The clock parameter vector obtained by the recursive estimation. This is an estimate of the frequency drift coefficient. This is the estimated phase offset value. For the first The timestamp recorded by the slave node in the secondary synchronization message. , For the first The master node's clock in the secondary synchronization message. Here is the gain matrix. Let be the error covariance matrix. Forgetting factor; Based on the frequency drift coefficient and the phase offset, a linear correction model for the local time of the slave node relative to the clock of the master node is calculated; the expression for the linear correction model is: ; Based on the linear correction model of all the slave nodes, a hierarchical clock synchronization tree with the master node as the root node is constructed to obtain the clock synchronization baseline.

3. The method according to claim 1, characterized in that, The acceleration time-series data of the entire blasting process recorded by the accelerometers built into each of the intelligent acquisition nodes, combined with a pre-trained crystal oscillator acceleration-sensitive model, is used to calculate the dynamic clock drift error caused by impact vibration. The acquired data is then corrected on the time axis using this dynamic clock drift error, resulting in the dynamically corrected time axis of each intelligent acquisition node and local time-based physical field data with the blasting zero moment corresponding to each intelligent acquisition node as the origin. This includes: Based on the crystal oscillator acceleration sensitivity model, the relative frequency deviation of the crystal oscillator is calculated using the acceleration time-series data of the entire blasting process; the expression of the crystal oscillator acceleration sensitivity model is: ,in, for The relative frequency deviation of the crystal oscillator at any given moment. This is the acceleration sensitivity coefficient matrix. For memory effect kernel function, for Acceleration time-series data for the entire blasting process at any given moment; The dynamic clock drift error caused by impact vibration is calculated based on the relative frequency deviation of the crystal oscillator; the formula for calculating the dynamic clock drift error is as follows: ; The time axis of the acquired data is nonlinearly corrected using the dynamic clock drift error to obtain the dynamically corrected time axis and the local time-based physical field data with the zero moment of the explosion as the origin; the expression for the nonlinear correction is: ,in, These are the time axis coordinates of each sampling point in the time axis. The time axis coordinates of the dynamically corrected time axis. This refers to dynamic clock drift error.

4. The method according to claim 1, characterized in that, Based on the clock synchronization baseline, each of the local time-based physical field data is uniformly registered to the dynamic correction time axis of the master node to obtain a multi-physics dataset, including: Using the dynamic correction time axis of the master node as a reference, the local time-based physical field data of each slave node are regarded as an asynchronous observation sequence; For each physical quantity data in the asynchronous observation sequence, optimal state estimation is performed using a Kalman filter to obtain the filtered value of each physical quantity; Based on the uniform sampling time point sequence on the dynamic correction time axis of the master node, the filtered values ​​of all physical quantities corresponding to each slave node are calculated by Lagrange interpolation to obtain the physical quantity estimate value corresponding to the uniform sampling time point sequence. The multiphysics dataset is obtained by aligning all the physical quantity estimates of all the intelligent acquisition nodes according to the dynamic correction time axis of the master node.

5. The method according to claim 1, characterized in that, The evolution features of the gas emission physical field, vibration propagation physical field, and stress redistribution physical field are extracted from the multiphysics dataset, and combined with the coupling features of the physical field pairs to obtain the blasting cycle feature vector, including: Based on the time series data of gas concentration after blasting from the multiphysics dataset, the initial rise slope, peak concentration, time to reach peak concentration, and decay time of gas emission are calculated to obtain the evolution characteristics of the gas emission physical field. Based on the blasting vibration signal from the multiphysics dataset, the arrival times of longitudinal and transverse waves are identified, and the peak velocity, dominant frequency distribution, energy duration, and response spectrum intensity are calculated to obtain the evolution characteristics of the vibration propagation physical field. Based on the surrounding rock stress-strain data of the multiphysics dataset, the adjustment magnitude, adjustment rate and number of abrupt change points of stress redistribution are calculated to obtain the evolution characteristics of the stress redistribution physical field. Based on the time delay between the vibration energy accumulation curve and the starting point of gas emission growth, the gas-vibration coupling characteristics are obtained; Based on the correlation coefficient between the stress unloading rate and the maximum peak value of the vibration wave, the stress-vibration coupling characteristics are obtained; Based on the evolution characteristics of the vibration propagation physical field, the evolution characteristics of the gas emission physical field, the stress redistribution physical field, the gas-vibration coupling characteristics, and the stress-vibration coupling characteristics, the blasting cycle feature vector is constructed.

6. An intelligent early warning system for safety risks in tunnel blasting construction, characterized in that, The system includes: The pre-blasting static calibration module is used to determine the frequency drift coefficient and phase offset of the clock of each slave node relative to the master node in the distributed sensor network of the tunnel to be blasted area through static clock drift calibration, and to obtain the clock synchronization baseline; the distributed sensor network includes multiple intelligent acquisition nodes; The blasting synchronization module is used to hard-trigger all the aforementioned intelligent acquisition nodes to simultaneously start data acquisition based on the trigger signal generated at the moment of blasting initiation by the non-contact electromagnetic induction trigger on the initiation network, obtain the acquired data, and record the local hardware timestamp of each of the aforementioned intelligent acquisition nodes as the zero moment of blasting; after the field programmable gate array corresponding to each of the aforementioned intelligent acquisition nodes detects the rising edge of the interrupt pin caused by the trigger signal, it freezes all analog-to-digital conversion channels of the corresponding intelligent acquisition node and latches the local hardware counter value as the local hardware timestamp corresponding to the blasting trigger; The post-blast drift compensation module is used to calculate the dynamic clock drift error caused by impact vibration based on the acceleration time series data of the entire blasting process recorded by the accelerometers built into each of the intelligent acquisition nodes, combined with a pre-trained crystal oscillator acceleration sensitive model, and to perform time axis correction on the acquired data through the dynamic clock drift error, so as to obtain the dynamic correction time axis of each intelligent acquisition node and the local time base physical field data of each intelligent acquisition node with the blasting zero time corresponding to the intelligent acquisition node as the origin; The configuration module is used to uniformly register each of the local time-based physical field data to the dynamic correction time axis of the master node based on the clock synchronization baseline, so as to obtain a multi-physics dataset. The feature module is used to extract the evolution features of the gas emission physical field, vibration propagation physical field, and stress redistribution physical field from the multi-physics dataset, and combine them with the coupling features of the physical field pairs to obtain the blasting cycle feature vector. The early warning module is used to input the feature vector of the blasting cycle into a pre-trained long short-term memory neural network early warning model to obtain the safety risk level of the current blasting cycle and the corresponding handling suggestions.

7. A computer device comprising a memory and a processor, wherein the memory stores a computer program, characterized in that, When the processor executes the computer program, it implements the method of any one of claims 1 to 5.

8. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the method of any one of claims 1 to 5.