Wave velocity non-contact testing and inversion characterization method for rock mass loading failure catastrophe early warning

Abstract

The application discloses a wave velocity non-contact testing and inversion characterization method for rock mass loading destruction disaster early warning, which uses water as a coupling medium, and realizes non-contact sound wave collection of the rock mass through a ring-shaped ultrasonic array, so that the problems of easy decoupling and signal interruption of the contact sensor are solved; STA / LTA combined with AIC criterion is adopted to pick up the wave form first arrival, a three-section water-rock refraction model is established, and Newton-Raphson iteration is used to correct the travel time error; the Tikhonov regularization SIRT algorithm is introduced to invert the two-dimensional wave velocity field of the rock mass; a ConvLSTM model is built to predict the wave velocity time evolution, a three-level early warning mechanism is established based on the wave velocity attenuation rate, spatial entropy and other indexes; a scalar damage variable is constructed from the wave velocity, a second-order damage tensor is obtained through multi-directional ray ellipsoid fitting, three stages of damage evolution are tracked, a disaster probability graph is generated through random forward propagation, and the shear, tension and other damage modes are automatically distinguished; the application can continuously monitor the damage evolution of the deep rock mass, identify the disaster high-risk area in advance, and is suitable for the instability prevention and control of the deep mineral engineering rock mass.

Description

Technical Field

[0001] This invention belongs to the field of rock mechanics and engineering safety monitoring technology. Specifically, it relates to a method for early warning of rock mass disasters and identification of dangerous areas by using ultrasonic waves for non-contact testing during the loading and failure process of rock mass, combined with wave velocity inversion and spatiotemporal index analysis. Background Technology

[0002] With the increasing depletion of shallow mineral resources, the development of deep resources has become a strategic development direction in my country's energy and mineral sectors. Deep rock masses exist in complex environments characterized by high stress, high temperature, and high permeability. This extreme environment significantly alters the mechanical response behavior of rocks, making them more susceptible to dynamic disasters such as rock bursts and large deformations during excavation. Therefore, studying the mechanical behavior of rocks in in-situ deep environments and establishing methods for monitoring and early warning of rock mass instability disasters are prerequisites for ensuring the safe and efficient development of deep resources.

[0003] Currently, wave velocity testing methods for rock mass disaster early warning mainly rely on traditional contact monitoring methods. Traditional wave velocity testing often uses sensors to directly contact the rock surface. During the rock mass loading and failure process, surface deformation, crack propagation, and debris shedding can easily lead to the failure of the coupling interface, causing signal interruption and making it impossible to obtain continuous wave velocity data throughout the entire failure process. Moreover, existing technologies mostly use a single transmission and reception path to obtain average wave velocity, lacking multi-angle, multi-channel acoustic travel time data acquisition and inversion methods. This makes it impossible to reconstruct the spatial distribution of the wave velocity field inside the rock mass, and it is also difficult to identify disaster precursors and locate dangerous areas through multi-index fusion.

[0004] To address the aforementioned problems, this invention proposes a non-contact wave velocity testing and inversion characterization method for early warning of rock mass loading failure. This method uses water as the sound wave transmission medium to achieve non-contact coupling between the ultrasonic probe and the rock sample, fundamentally solving the signal failure problem of traditional contact sensors during the rock mass failure stage. Through a multi-channel, multi-angle sound wave transmitting and receiving array, combined with an inversion algorithm, the spatial distribution image of the wave velocity field is reconstructed. Based on the temporal anomaly changes of the wave velocity field and the identification of spatial low-velocity anomaly areas, a multi-index fusion early warning system is constructed to identify precursors of rock mass instability and locate dangerous areas. This invention provides a scientific and reliable testing and characterization method for early warning of disasters in deep rock mass engineering. Summary of the Invention

[0005] (I) Purpose of the Invention

[0006] This invention addresses the technical problems in deep rock mass disaster monitoring, such as the failure of traditional contact-based wave velocity testing due to rock surface deformation, the inability of single-path testing to invert the spatial distribution of wave velocity, and the difficulty of accurately predicting disasters and locating dangerous areas using single indicators. It proposes a non-contact wave velocity testing and inversion characterization method for early warning of rock mass loading and failure disasters. This invention achieves non-contact acoustic coupling through a water medium, ensuring continuous signal acquisition throughout the entire rock mass failure process. It reconstructs the spatial distribution of the wave velocity field using multi-angle acoustic travel time data and inversion algorithms. Based on this, it integrates wave velocity temporal anomalies and spatial low-velocity zone identification to identify precursors of rock mass instability and locate dangerous areas.

[0007] (II) Technical Content

[0008] To solve the above-mentioned technical problems, the technical solution of the present invention is: a non-contact wave velocity testing and inversion characterization method for early warning of rock mass loading failure disasters, comprising the following steps:

[0009] S1. Construct a non-contact ultrasonic testing experimental system: After waterproofing the rock specimen, immerse it in a water tank, using water as the coupling medium; arrange an ultrasonic transducer array in a ring on the outer wall of the water tank to form a non-contact sound wave transmission path; connect the transducer array to an ultrasonic pulse transmitter and receiver and a data processing terminal to establish a synchronous triggering and multi-channel parallel acquisition link.

[0010] S2. Multi-angle acoustic travel time data acquisition: Each transmitting transducer is excited in sequence, and all receiving transducers synchronously acquire the transmitted waveform to form a multi-angle intersecting ray network covering the interior of the rock mass and pick up the first arrival travel time of each waveform.

[0011] S3. Iterative solution and travel time correction of the refraction path at the water-rock interface: In order to eliminate the systematic travel time error caused by strong refraction at the water-rock interface, an iterative solution method for the refraction path based on Snell's law and Fermat's principle is established to obtain the true travel time inside the rock.

[0012] S4. Wave velocity field tomographic inversion based on regularized joint iterative reconstruction: The rock mass section is discretized into a grid model, and a linear forward modeling equation for travel time and slowness is established; the actual travel time inside the rock is used as the observation value, and the slowness is updated by joint iterative reconstruction; Tikhonov regularization and first-order difference smoothing constraints are introduced in the iteration process, and spatial smoothing constraints are applied to the slowness field through explicit regularization steps; after the iteration converges, the spatial distribution of wave velocity inside the rock mass is obtained and visualized in the form of pseudo-color cloud map.

[0013] S5. Multi-indicator fusion disaster early warning based on spatiotemporal sequence deep learning prediction: Repeat steps S2 to S4 during the loading process to obtain a dynamic evolution image sequence of the wave velocity field; construct a spatiotemporal sequence prediction model based on a convolutional long short-term memory network, using historical wave velocity field image frames as input, to predict the evolution trend of the wave velocity field in future multiple frames; based on the predicted wave velocity field sequence, extract multi-dimensional early warning indicators such as average wave velocity attenuation rate, low wave velocity area growth rate, and wave velocity field spatial entropy, and establish multi-level early warning criteria based on the combination of indicator thresholds to achieve advanced identification of precursors to rock mass instability;

[0014] S6. Spatial localization and damage mode identification of hazard zones based on damage tensor evolution: The wave velocity field obtained from inversion and prediction is converted into a second-order damage tensor field; ellipsoid fitting is performed based on the wave velocity square ratio of ray clusters in different directions to obtain the principal values ​​and principal directions of the damage tensor; connected component labeling and multi-frame correlation tracking are performed on the damage field sequence to extract the geometric evolution features of the damage region; based on the change law of damage features, the three stages of damage initiation, stable expansion and accelerated localization are quantitatively divided; the damage mode is automatically identified based on the principal value ratio and principal direction features of the damage tensor; and a hazard probability map is generated by multiple random forward propagation of ConvLSTM to delineate high and medium risk areas and output hazard zone information.

[0015] Furthermore, in step S2, the initial arrival time is automatically picked up using a combination of the long short window ratio method (STA / LTA) and the Akaike Information Criterion (AIC), supplemented by manual correction.

[0016] Furthermore, in step S3, the iterative solution method for the refraction path includes: using the locations of the water-rock incident point and the exit point as variables, constructing a three-segment propagation path model that includes a straight segment in the water and a curved segment within the rock; using the shortest path method based on Dijkstra's algorithm to trace the curved ray within the rock; using the interface Snell's law residual as the objective function, iteratively solving the coordinates of the incident point and the exit point using the Newton-Raphson method, and obtaining the true sound wave travel time inside the rock after deducting the propagation time of the water medium, as the observation data for wave velocity field inversion.

[0017] Furthermore, in step S3, the complete sound wave propagation path is divided into three segments: the incident segment in water, the transmission segment inside the rock, and the exit segment in water. The method for establishing the three-segment refraction path model is as follows:

[0018] A horizontal section at the midpoint of the rock specimen's height is selected as the analysis plane. A two-dimensional rectangular coordinate system is established with the center of the specimen's cross-section as the origin. The radius of the rock specimen is set as... The boundary equation of the cross section is Let the coordinates of the transmitting probe be... The coordinates of the receiving probe are The wave velocity in the water medium is The current wave velocity field inside the rock is ;

[0019] The unknown quantity to be solved is the incident point of the sound wave at the water-rock interface. and the ejection point at the rock-water interface Both points are located on the circumference of the specimen, satisfying the condition that... , ;

[0020] Given and After setting the coordinates, calculate the travel time in the water for the two segments:

[0021]

[0022]

[0023] In the formula, For the travel time of the sound wave incident in water, The travel time of the sound wave emitted from the water segment.

[0024] Furthermore, step S3 uses the shortest path method combined with Dijkstra's algorithm to solve for the curved rays within the rock, specifically including:

[0025] The wave velocity field inside the rock is non-uniformly distributed, and the sound wave propagation satisfies the equation of process function. The rock cross section is discretized into a square grid and each grid is assigned a slowness. The grid node topology is constructed, and the edge propagation travel time is defined by the node distance and the average slowness of the grid.

[0026] The Dijkstra algorithm is run with the water-rock incident point as the source point. The grid nodes are traversed and the shortest propagation time is updated. The calculation is terminated early when the corresponding exit node is reached. The discrete ray path nodes are obtained by backtracking in reverse.

[0027] A continuous and smooth acoustic curve within the rock is obtained by using cubic spline interpolation to smooth the discrete node sequence. The interpolation process satisfies the conditions of coordinate continuity, second derivative continuity and natural boundary constraints.

[0028] Furthermore, in step S4, the slowness update formula of the joint iterative reconstruction method is:

[0029]

[0030]

[0031] in, The total number of rays passing through the j-th grid. Let be the path length of the i-th ray within the j-th grid.

[0032] Furthermore, in step S4, the explicit regularization smoothing step is as follows:

[0033] The regularization problem is transformed into an iterative form: in each iteration, the result is first updated using the joint iterative reconstruction method. Then solve for a smoothing step:

[0034]

[0035] The problem has an explicit solution:

[0036]

[0037] Construct sparse matrices :

[0038]

[0039] Construct the solution equation:

[0040]

[0041] It is a sparse symmetric positive definite matrix, which can be solved quickly using the conjugate gradient method without explicit inversion. The conjugate gradient method (CG) solver is called once after each iteration of the joint iterative reconstruction method to obtain the regularized slowness. In the above formula, λ is the regularization parameter, and L is the first-order difference regularization matrix.

[0042] Furthermore, in step S5, a spatiotemporal sequence prediction model based on a convolutional long short-term memory network is constructed, wherein the convolutional long short-term memory network adopts a symmetrical structure composed of an encoder and a decoder.

[0043] The encoder is composed of multiple layers of ConvLSTM (Convolutional Long Short-Term Memory) network, which extract the spatiotemporal features of the input wave velocity field sequence layer by layer. The input of the first ConvLSTM layer is the normalized wave velocity field image frame, and the input of each subsequent layer is the hidden state of the previous layer. The encoder finally outputs the hidden state and cell state of the last time step, which serve as the initial spatiotemporal features of the decoder.

[0044] The decoder stacks ConvLSTM layers of the same number, with a structure symmetrical to the encoder. The initial state directly copies the output features of the encoder. In the prediction stage, the decoder iterates to generate the wave velocity field for future multiple frames. The prediction result of each step is fed into the next layer of the network, and finally, a single-channel wave velocity field prediction image is obtained through mapping.

[0045] During offline training, the model employs a composite loss function that integrates mean squared error, image structure similarity, and gradient difference. The Adam optimizer is used for iterative convergence, and the trained network is used for real-time wave velocity field advance prediction.

[0046] Furthermore, in step S6, the wave velocity field is converted into a damage variable field, including:

[0047] Define wave velocity damage variable :

[0048]

[0049] In the formula: For a moment In position The inversion wave velocity value at that location; This represents the wave velocity at that location in its initial, undamaged state.

[0050] (III) Technical Effects

[0051] The advantages of this invention compared to the prior art are:

[0052] 1. This invention achieves continuous and stable acquisition of wave velocity data throughout the entire loading and destruction process through a non-contact testing system.

[0053] 2. This invention develops a regularized joint iterative reconstruction tomographic imaging inversion method, which realizes the visual characterization of the spatial distribution of wave velocity field inside the rock mass. It expands the wave velocity information of the traditional method from "point" or "line" to the distribution image of "surface" or "volume", revealing the spatial evolution process of non-uniform damage inside the rock mass.

[0054] 3. This invention constructs a spatiotemporal prediction model of a convolutional long short-term memory network that integrates physical constraints, realizing advanced prediction and multi-level early warning of the evolution trend of wave velocity field.

[0055] 4. This invention proposes a method for spatial localization of hazardous areas and identification of damage modes based on the evolution of second-order damage tensors, which realizes the delineation of high-risk disaster areas and the automatic identification of damage modes.

[0056] 5. This invention organically integrates non-contact acoustic testing, travel time extraction considering interface refraction, regularized tomographic inversion, spatiotemporal sequence prediction constrained by physical information, and damage tensor spatial localization, forming a complete closed-loop technology system of "sensing-inversion-prediction-localization", providing a scientific, systematic, and reliable technical solution for disaster early warning in deep rock mass engineering. Attached Figure Description

[0057] Figure 1 This is a schematic diagram of the wave velocity measurement point arrangement in the system of the present invention.

[0058] Figure 2 This is a schematic diagram of the sound wave propagation path of the present invention.

[0059] Figure 3 This is a flowchart of the wave velocity field inversion and spatial distribution characterization of the present invention.

[0060] Figure label:

[0061] 11. Axial pressure; 22. Rock mass sample; 33. Liquid region; 44. Water tank; 55. Wave velocity transducer; 1. Emission point T; 2. Rock mass incident point P1; 3. Center point O; 4. Rock mass ejection point P2; 5. Receiver point Rr. Specific Implementation

[0062] The present invention will now be described in further detail with reference to the accompanying drawings and embodiments.

[0063] Reference Appendix Figure 1 To be continued Figure 3 A non-contact wave velocity testing and inversion characterization method for early warning of rock mass loading failure disasters includes the following steps:

[0064] Step S1: Constructing a non-contact ultrasonic testing experimental system

[0065] The rock to be tested was processed into a standard specimen, and after the sides were coated with a flexible waterproof sealing layer, it was completely immersed in a transparent water tank filled with degassed constant-temperature water, using water as the acoustic coupling medium. Multiple ultrasonic transducers were arranged in a ring on the outer wall of the water tank, and an acoustic impedance matching layer was coated between the transducers and the tank wall to enhance the transmission efficiency, forming a non-contact acoustic wave transmission path of "transducer-tank wall-water-rock-water-tank wall-transducer". The transducer array was connected to an ultrasonic pulse transmitter and receiver and a data processing terminal through a multi-channel switching device to establish a synchronous triggering and multi-channel parallel acquisition link.

[0066] Step S2: Multi-angle acoustic travel time data acquisition

[0067] Each transmitting transducer is excited sequentially, and all receiving transducers synchronously acquire the transmitted waveforms, forming a multi-angle intersecting ray network covering the interior of the rock mass. The first arrival time of each waveform is automatically picked up using the long-short time window ratio method and the Akaike information criterion, and then manually corrected to obtain the measured total travel time.

[0068] Step S3: Iterative solution and travel time correction of the refraction path at the water-rock interface

[0069] To eliminate the system travel time error caused by strong refraction at the water-rock interface, this invention establishes an iterative and accurate solution method for the refraction path based on Snell's law and Fermat's principle: using the locations of the water-rock incident and exit points as variables, a three-segment propagation path model is constructed, including a straight segment in the water and a curved segment within the rock; the path within the rock is traced using the shortest path method based on Dijkstra's algorithm; using the residual of Snell's law at the interface as the objective function, the Newton-Raphson method is used to iteratively solve for the accurate coordinates of the incident and exit points, and after deducting the propagation time in the water medium, the true travel time of the sound wave inside the rock is obtained, which serves as the observation data for wave velocity field inversion.

[0070] Step S4: Wave velocity field tomography inversion based on regularized joint iterative reconstruction

[0071] The rock mass cross section is discretized into a grid model, and linear forward equations for travel time and slowness are established. The accurate internal travel time of the rock obtained in steps S2 and S3 is used as the observation value, and the slowness is updated by joint iterative reconstruction method. In order to solve the ill-conditioned nature of the inverse problem, Tikhonov regularization and first-order difference smoothing constraint are introduced in the iteration process. The slowness field is subjected to spatial smoothing constraint through explicit regularization step, which suppresses false high-frequency oscillations while ensuring convergence stability. After the iteration converges, the spatial distribution of wave velocity inside the rock mass is obtained and visualized in the form of pseudo-color cloud map.

[0072] Step S5: Multi-indicator fusion disaster early warning based on spatiotemporal sequence deep learning prediction

[0073] Throughout the loading process, steps S2 to S4 are repeated at fixed time intervals to obtain a sequence of dynamic evolution images of the wave velocity field. A spatiotemporal sequence prediction model based on a convolutional long short-term memory network is constructed, using historical wave velocity field image frames as input to predict the evolution trend of the wave velocity field in future multiple frames. The network loss function introduces structural similarity loss and gradient difference loss on the basis of mean square error, and embeds physical constraint terms based on the damage mechanics evolution equation to enhance the prediction generalization ability under small sample conditions.

[0074] Based on the predicted wave velocity field sequence, multi-dimensional early warning indicators such as average wave velocity attenuation rate, low wave velocity area growth rate, and wave velocity field spatial entropy are extracted and input into a lightweight classification network to output macroscopic damage probability. Based on the combination of damage probability and indicator thresholds, a three-level early warning criterion of yellow, orange, and red is established to achieve accurate and advanced identification of precursors to rock mass instability.

[0075] Step S6: Spatial localization of hazard zones and identification of damage patterns based on damage tensor evolution

[0076] The wave velocity field obtained from inversion and prediction is converted into a second-order damage tensor field: ellipsoid fitting is performed based on the wave velocity square ratio of ray clusters in different directions to obtain the principal value and principal direction of the damage tensor, overcoming the limitation that scalar damage models cannot describe anisotropic damage; connected component labeling and multi-frame correlation tracking are performed on the damage field sequence to extract geometric evolution features such as damage area, center coordinates, aspect ratio, principal direction angle, and localization index.

[0077] According to the variation laws of the damage area growth rate, localization index and aspect ratio, the three stages of damage initiation, stable propagation and accelerated localization are quantitatively divided; based on the principal value ratio and principal direction characteristics of the damage tensor, the shear failure, tensile splitting failure or mixed failure mode is automatically identified; combined with multiple random forward propagations of ConvLSTM to generate a catastrophe probability map, high and medium risk areas are delineated, and the center coordinates, range and predicted fracture arrival time of the dangerous area are output, so as to provide a clear spatial target area and time window for precise engineering prevention and control.

[0078] To make the objectives, technical solutions and advantages of the present invention clearer, the present invention is further described in detail below with reference to specific embodiments and the accompanying drawings.

[0079] Step S1: Construct a non-contact ultrasonic testing experimental system

[0080] This step aims to build an experimental platform capable of realizing non-contact wave velocity testing throughout the whole rock loading process, so as to fundamentally avoid the signal failure problem of traditional contact sensors in the rock mass failure stage.

[0081] (1) Prepare rock specimens

[0082] Select sandstone with good integrity, process the rock into standard cylindrical specimens with a diameter of 50mm and a height of 100mm in accordance with the standards of the International Society for Rock Mechanics, grind both end faces of the specimens flat with sandpaper, with a parallelism error not exceeding ±0.05mm to ensure uniform force during loading, place the specimens in an oven at 105°C for 24 hours for drying, take them out and cool them to room temperature, and prevent the specimens from absorbing moisture in the air during the cooling process.

[0083] Prepare a flexible waterproof sealing material: use polyurethane sealant, which still maintains a certain elasticity after curing, can adapt to the lateral deformation of the rock during loading, and avoid cracking of the sealing layer. Use a scraper to uniformly coat the sealant on the side surface of the specimen to a thickness of 1mm, and cure it at room temperature for 24 hours. Sealing effect inspection: fully immerse the sealed specimen in water for 24 hours, and it is qualified if the mass increase is less than 0.1%. The side surface of the specimen is sealed with a waterproof adhesive film to prevent water from seeping into the rock pores during the water immersion process and affecting the test results.

[0084] (2) Prepare the water tank and coupling medium

[0085] Fabricate a water tank made of transparent acrylic, with a diameter of 100mm, a height of 200mm, and a wall thickness of 10mm. The top of the water tank is open. To avoid direct contact between the press indenter and water, a layer of flexible waterproof gasket with a thickness of 2mm is arranged between the indenter and the end face of the specimen, and a retractable corrugated sealing cover is installed on the top of the water tank, so that the indenter is always isolated from water during the loading process, and the liquid level in the water tank is kept stable at the same time. The bottom of the water tank is flat and placed on the lower indenter of the press.

[0086] Pour deaerated purified water into the water tank to a height of 150 mm, ensuring that the rock specimen is completely submerged. Maintain the water temperature at 20±1℃. Replace the water before each experiment to reduce the effect of temperature on wave velocity.

[0087] (3) Arrange the ultrasonic probe array

[0088] An array of ultrasonic probes is arranged on the outer wall surface of the water tank. In this embodiment, a total of 16 probes are used, all of which can transmit and receive, and all of them are broadband ultrasonic probes with a center frequency of 500kHz. The probes are tightly attached to the tank wall by acrylic fixing brackets. A thin layer of ultrasonic coupling agent is coated on the probe end face, and the probe end face is tightly pressed against the outer wall of the water tank by the brackets to ensure that there is no air gap between the probe and the water tank wall.

[0089] The probes are arranged in a ring at equal angles around the rock specimen. With the center of the specimen as the center, 16 transmitting probes are evenly distributed on the circumference in the horizontal plane (the angle between adjacent probes is 22.5°). None of the probes touch the surface of the rock specimen. The center height of the probes is aligned with the midpoint of the height of the rock specimen to ensure that the sound waves mainly pass through the middle section of the specimen. Each probe is numbered T1-T16, and the spatial coordinates of each probe are recorded for subsequent ray tracing.

[0090] (4) Connect to the data acquisition system

[0091] Connect the probe to the transmitter output of the multi-channel switcher via a coaxial cable, and simultaneously connect the probe to the receiver input of the switcher. Connect the control line of the switcher to the USB interface of the computer. Connect the transmitter trigger line and receiver signal line of the switcher to the ultrasonic pulse transmitter and receiver. Connect the waveform output of the transmitter and receiver to the high-speed data acquisition card via a BNC cable, and then connect it to the computer via USB. Ensure that all electrical connections are reliable and set the synchronous trigger signal.

[0092] Set the transmitter and receiver to "external trigger" mode, and have the computer send the trigger signal through a switching device.

[0093] Triggering sequence: The computer sends a "channel switching" command → after switching is completed, it sends a "transmission trigger" pulse → the transmitter and receiver generate a high-voltage pulse to excite the transmitter probe → all receiver probes acquire waveforms → the data is transmitted to the computer → the cycle continues to the next transmitter probe.

[0094] Step S2: Multi-angle testing of sound waves and acquisition of travel time data

[0095] This step aims to obtain travel time data for multiple acoustic rays covering the interior of the rock specimen.

[0096] (1) Set the acquisition parameters

[0097] In the computer control software, set the sampling frequency to 10MHz, the transmission voltage to 400 to 800V (the specific value should be selected based on the signal-to-noise ratio of the received signal), the pulse repetition frequency to 100Hz, and the signal gain to 40dB. Select the one-to-many transmit mode.

[0098] (2) Perform automatic scanning

[0099] A computer-controlled multi-channel switching device sequentially transmits the i-th probe (i=1,2,…,16). For each transmission, all receiving probes synchronously start acquiring and recording the waveform data of the direct wave or transmitted wave. The 16 waveform data, along with the current transmitting probe number, receiving probe number, and timestamp, are saved to the hard disk. Then, the device switches to the next transmitting probe and repeats the above process. Each sound wave propagation path corresponds to a "transmit-receive" pair, resulting in a total of 16×16=256 waveform data points. Paths that deviate significantly or have excessive signal attenuation can be removed in subsequent processing.

[0100] (3) Pick up the items when you first arrive and leave

[0101] On each acquired waveform curve, the first arrival point is automatically selected using an automatic algorithm based on the short-term window ratio: short-term window length: 5 sampling points; long-term window length: 50 sampling points; the STA / LTA ratio is calculated for each sampling point, and when the ratio exceeds a preset threshold (3.5 in this embodiment), the point is determined to be the first arrival point. The Akaike Information Criterion (AIC) is then used for refined selection: within a range of 50 sampling points before and after the initially estimated first arrival point based on STA / LTA, the first arrival point is calculated... function:

[0102]

[0103] In the formula, This represents the number of sampling points within the window. From the 1st to the 2nd A sequence of observations; For from the first One to the first A sequence of observations; This represents the sample variance. Minimum value corresponding to This is the precise first arrival point.

[0104] The automatically acquired results are marked with vertical lines on the waveform diagram. Combined with manual correction, waveforms with low signal-to-noise ratios or obviously incorrect automatically acquired waveforms are manually re-acquired or marked as invalid. The system accurately captures the initial arrival time of the ultrasonic waves from the transmitting probe through the water and rock specimens, back to the receiving probe. .

[0105] Step S3: Iterative and accurate solution of the refraction path at the water-rock interface and travel time correction

[0106] Since both the transmitting and receiving probes are placed on the outer wall of the water tank, the sound wave needs to penetrate multiple layers of medium in sequence: "water tank wall-water-rock-water-water tank wall". The sound velocity difference between water and rock is significant (the P-wave velocity of water is approximately 1480 m / s, while the P-wave velocity of rock is approximately 3500 m / s), causing strong refraction of the sound wave at the water-rock interface. If the sound wave path is simplified to a straight line, a non-negligible system travel time error will be introduced, severely affecting the accuracy of subsequent wave velocity field inversion. Therefore, this invention establishes an iterative and accurate solution method for the refraction path based on Snell's law and Fermat's principle. The specific steps are as follows:

[0107] (1) Establishing a three-layer medium geometric model and defining variables

[0108] A two-dimensional rectangular coordinate system is established using the horizontal section at the midpoint of the rock specimen's height as the analysis plane, with the origin at... Located at the center of the cross-section of the rock specimen, let the radius of the rock specimen be... (Unit: mm), boundary equation is The coordinates of the transmitting probe are: The coordinates of the receiving probe are The wave velocity in the water medium is 1480 m / s; the current wave velocity field inside the rock is (Initially a uniform model with a wave velocity of 3500 m / s, subsequently updated by inversion iteration).

[0109] The unknown to be solved is the precise incident point of the sound wave at the water-rock interface. and the precise exit point at the rock-water interface Both points lie on the circumference and satisfy the condition that... , .

[0110] The total propagation path of sound waves consists of three segments: the incident segment in water. Internal transmission section of rock underwater launch section This invention focuses on solving the problem of accurate tracking of these three paths.

[0111] (2) Calculation of travel time on the water path

[0112] For a given and The travel times of the two straight-line propagation paths in the water are as follows:

[0113]

[0114]

[0115] (3) Tracing and travel time calculation of bending rays inside the rock

[0116] Internal wave velocity field of rock It is non-uniform, sound waves from spread to The actual path is a curved curve with the shortest travel time, satisfying the equation:

[0117]

[0118] In the formula, This represents the wavefront arrival time field.

[0119] This invention uses the Shortest Path Method (SPM) combined with Dijkstra's algorithm to solve the above equations. The specific steps are as follows:

[0120] 1) Mesh discretization and slow field construction

[0121] The square area containing the rock cross-section Divided into A square grid cell, grid side length To form a set of grid nodes ,in:

[0122] , ,

[0123] Each grid cell Give slowness Its value is determined by the wave velocity field in the current inversion iteration:

[0124]

[0125] 2) Node connection graph construction

[0126] To apply the shortest path method, connecting edges need to be established between nodes. This invention adopts a forward star topology, where each node is connected to multiple neighboring nodes to cover multi-directional propagation.

[0127] Establish connecting edges between grid nodes to form a directed graph and node set. Includes all grid intersections located within and on the boundaries of the rock region, for nodes and The traversal weight of an edge between two points is defined as follows:

[0128]

[0129]

[0130] In the formula: For nodes and Spacing; For the edge The average slowness of the traversed grid is calculated by taking the slowness value of the grid containing the midpoint of the edge. .

[0131] 3) Dijkstra's algorithm for finding the shortest travel time

[0132] At the point of incidence The source point (approximately bound to the nearest mesh node) ), execute Dijkstra's algorithm:

[0133] Initialization: Set source point Time of travel The time travel of all other nodes All nodes Establish a priority queue ,Will Join the team.

[0134] Loop: when If the node is not empty, pop the node with the smallest running time. ,like If the node has already been processed, skip it and proceed with the next iteration. All adjacent nodes Calculate candidate travel time In the formula Calculated from the above definition of travel time weight, if Then update ,Record predecessor node and will Join the team.

[0135] Relaxation: If That is, the target node (and If the nearest grid node is not found, the algorithm terminates prematurely. After the algorithm ends, the result is obtained from... arrive Shortest travel time inside the rock:

[0136]

[0137] 4) Ray path tracing and emission direction extraction

[0138] from Begin, along Array backtracking to The node sequence is obtained as follows:

[0139]

[0140] A smooth ray curve is obtained by performing cubic spline interpolation on the node sequence. .

[0141] Let the cumulative arc length parameter be... for:

[0142]

[0143] To each and Coordinates about Constructing a cubic spline function ,satisfy:

[0144] Interpolation conditions: , = ;

[0145] Second derivative continuity: , Continuous within internal nodes;

[0146] Natural boundary conditions: , ;

[0147] The smoothed path is:

[0148]

[0149] (4) Snell's law of refraction and residual function at the water-rock interface

[0150] To calculate the refraction angle at the interface, it is necessary to extract the tangent direction of the ray at its endpoint and take the ray path. Up close interpolation points Then the unit vector of the ejection direction inside the rock is:

[0151]

[0152] Similarly, take the ray path Up close interpolation points The unit vector of the incident direction inside the rock is obtained as follows:

[0153]

[0154] At the point of incidence At this point, the unit normal vector of the interface (pointing to the center of the circle) is:

[0155]

[0156] At the point of incidence At that location, the unit vector of the incident direction in the water is:

[0157]

[0158] Then the angle of incidence ( The angle between the ray and the normal (at a given point) satisfies:

[0159]

[0160]

[0161] Angle of refraction ( The angle between the direction of the projectile from the rock side and the normal satisfies:

[0162]

[0163]

[0164] According to Snell's law, we should have:

[0165]

[0166] in for Rock wave velocity within the grid. Define the Snell residual at the incident end. for:

[0167]

[0168] Similarly, at the point of departure At that location, the normal vector is The angle of incidence on the rock side is The water-side refraction angle is Snell residual at the launch end for:

[0169]

[0170] (5) Iterative solution of refraction point based on Newton-Raphson method

[0171] and The positions of these elements are interdependent and cannot be directly solved. This invention parameterizes them as angular coordinates on a circle:

[0172]

[0173] The problem is transformed into finding the angle. and This makes the residual vector and Approaching zero, the Newton-Raphson method is used for iterative solution:

[0174] 1) Initial guess

[0175] With the transmitting probe The line connecting the center of the circle intersects the rock boundary at... To receive probe The line connecting the center of the circle intersects the rock boundary at... The corresponding angle is:

[0176]

[0177]

[0178] 2) Iterative update formula

[0179] In the In the next iteration, the following steps are performed:

[0180] The first step is to determine the current situation. , calculate , ;

[0181] The second step is to call the shortest path method described above to calculate... , , ;

[0182] The third step is to calculate the residuals. , ;

[0183] Fourth step, if If the iteration converges, stop iterating; otherwise, proceed to the next step.

[0184] Step 5: Calculate the Jacobian matrix. The numerical approximation is for Apply small perturbation Recalculate the internal ray tracing to obtain the perturbed residual. , ,thereby:

[0185]

[0186] Similarly, Apply a perturbation to obtain , and , The Jacobian matrix G is formed as follows:

[0187]

[0188] Step 6: Solve the linear equation system G to obtain the angle correction. , ;

[0189] Step 7, update the angle:

[0190]

[0191] Step 8, Order Go back to step one.

[0192] 3) Post-convergence processing

[0193] After iterative convergence, the precise incident point is obtained. and launch point and the actual travel time inside the rock. The travel time inside the rock after deducting the water propagation time is:

[0194]

[0195] Total arrival time After deducting the time spent traveling in water, this is the result. This travel time will be taken as the final observation travel time of this ray. It is used for subsequent wave velocity field inversion.

[0196] (6) Algorithm effect description

[0197] Through the iterative solution described above, this invention reconstructs the actual propagation path of each acoustic ray in the water-rock-water multilayer medium, eliminating the systematic travel time error caused by simplifying the ray to a straight line due to neglecting refraction in traditional methods. Using this method, the root mean square value of the travel time picking residuals is less erroneous than the straight-line assumption, providing a high-precision data foundation for subsequent tomographic imaging inversion.

[0198] Step S4: Wave velocity field inversion and spatial distribution pattern characterization

[0199] This step utilizes the precise rock internal travel time data obtained in the previous steps. The spatial distribution of wave velocity field inside the rock specimen was reconstructed using a tomographic inversion algorithm.

[0200] (1) Establish a discretized forward model

[0201] In this embodiment, the horizontal section at the midpoint of the rock specimen (at a height of 50 mm) is used as the inversion plane. The discretized grid of the inversion region is consistent with the grid used for internal ray tracing in step S3 to ensure seamless data transfer.

[0202] Based on the actual dimensions of the rock specimen (φ50mm×100mm), a two-dimensional initial wave velocity model was established, and the cross-section at the midpoint of the rock specimen's height was discretized as follows: A regularized grid cell, which is taken in this embodiment. Each grid cell is 1mm x 1mm in size, and each grid cell is numbered. .in This represents the total number of effective cell grids. Based on rock type and experience, an initial wave velocity value is set, with the wave velocity within each grid cell being uniformly distributed. In this embodiment, the initial P-wave velocity for sandstone is set to 3500 m / s.

[0203] In travel-time tomography, travel time With slow It exhibits a linear relationship with wave speed, but a non-linear relationship with wave velocity; therefore, a slower wave speed is chosen. As unknown parameters in the inversion, the forward modeling problem is linearized.

[0204]

[0205] In the formula, It represents the slowness, which is the reciprocal of the wave speed.

[0206] Combine the slowness of all grids into a vector. , where is the unknown quantity to be solved, and the initial slowness is 0.2857 μs / mm.

[0207] Establish the forward equation, for the th ray ( , (Total effective rays), theoretical travel time Equal to the line integral of the slowness along the ray path:

[0208]

[0209] In the formula, For the first The first ray Path length within each grid cell This is the ray path.

[0210] Write the equations for all rays in matrix form:

[0211]

[0212] In the formula, for The Jacobian matrix, whose elements S is the slowness matrix; for The theoretical travel vector.

[0213] Actual timekeeping Keeping up with the theory residuals between for:

[0214]

[0215] The purpose of inversion is to find the slowness vector matrix. This makes the residual vector Minimum.

[0216] (2) Ray path tracing

[0217] This embodiment employs the Shortest Path Method (SPM) to solve for the ray path. To mitigate the "sawtooth effect" of SPM on a regular grid, linear interpolation is performed on the slowness field before each ray tracing, ensuring continuous variation of the slowness within the grid. Simultaneously, after tracing the path back, cubic spline interpolation is used to smooth the polyline path, resulting in a more realistic curved ray path. The Shortest Path Method discretizes the continuous medium into a node network and calculates the shortest path from the emission point to the receiving point using Dijkstra's algorithm in graph theory, serving as an approximate path for the acoustic ray.

[0218] For each ray whose refraction point has been accurately determined, its smooth path within the rock is... The path length of the ray within each grid cell was obtained in step three using cubic spline interpolation. The following steps are adopted:

[0219] Mesh nodes are set based on the four points around the perimeter of the mesh in the forward model. This embodiment uses the same nodes as the slow-speed mesh: X direction -25mm to 25mm, step size 1mm; Y direction -25mm to 25mm, step size 1mm. A total of 51 × 51 = 2601 nodes are used. Each node... Spatial coordinates are .

[0220] The travel time between adjacent nodes (horizontally, vertically, and diagonally) is equal to the distance multiplied by the average slowness of the grid traversed by that path. For connecting nodes... and The grid it passes through can be determined by linear interpolation. To simplify the calculation, the slowness of the grid where the midpoint of the line segment is located is usually taken as the slowness of the line segment. Then the travel time is:

[0221]

[0222]

[0223] In the formula, For nodes and Distance between The slowness of the grid containing the midpoint of the line segment.

[0224] As shown in the ray path tracing method in step S3, the Dijkstra algorithm is used to calculate the shortest travel time field. Taking the position of the transmitting probe as the source point (travel time initialized to 0), the travel time of all other nodes is initialized to infinity. A priority queue (min-heap) is used to maintain the nodes to be processed. Each time, the node with the smallest travel time is taken from the queue, and the travel time of all its adjacent nodes is updated.

[0225]

[0226] Repeat until all nodes have been processed or the receiving point converges.

[0227] Ray path backtracking starts from the receiver point and traces back to the source point along the direction of the time gradient descent: At the current node, check all adjacent nodes and select those that satisfy the following conditions: The adjacent node is taken as the previous node, and the process is repeated until the source point is reached. The coordinates of all nodes passed on the path are recorded and connected in order to form a polyline, which is the ray path.

[0228] For each ray After obtaining the path, it is necessary to calculate the length of the path within each grid cell. .

[0229] The ray path consists of a series of line segments (straight line segments between adjacent nodes). Each line segment may traverse multiple grids. For each line segment, a grid traversal algorithm is used to calculate the length of all grids traversed by the line segment and their interiors, given the endpoints of the line segments. and Calculate the intersection parameters of the line segment and the grid lines, divide the line segment into segments according to the intersection points, ensuring each segment lies entirely within a grid, calculate the length of each segment, and sum the total lengths within each grid to obtain the result. .

[0230] initialization sparse matrix All zeros for each ray Iterate through all the traversed grids. , calculate deposit .

[0231] (3) Inversion Iterative Solution

[0232] Because the number of rays is much smaller than the number of grids, and the ray distribution is uneven, the sensitivity matrix... The problem is highly ill-conditioned, and directly solving the least squares problem leads to non-unique solutions and extreme sensitivity to noise. To address this, this invention introduces Tikhonov regularization, which transforms the ill-conditioned problem into a well-posed problem by applying model smoothing constraints.

[0233] Constructing the least squares objective function Find the slowness vector This minimizes the L2 norm of the travel time residuals:

[0234]

[0235] Its solution satisfies the equation ,but It is usually pathological or singular, and difficult to solve directly.

[0236] This invention employs Tikhonov regularization in the objective function. Add model smoothness constraints:

[0237]

[0238] In the formula, The regularization parameter controls the balance between data fitting and model smoothing. This is the regularization matrix. Regularization terms. In effect, it is a penalty for the slow rate of change in space, forcing the inversion results to be smooth.

[0239] This embodiment employs a first-order difference operator to promote a smooth, gradual change in the speed of adjacent grids, suppressing spurious high-frequency oscillations. Regularization parameter. The choice of [aspect / factor] affects the inversion quality; this embodiment uses the L-curve method for automatic selection. For the data size of this embodiment, The initial value is set to 0.05, and then adjusted according to the inversion results: if the inverted wave velocity field shows violent oscillations, increase the value. If the residual is too large and the smoothness is too strong, reduce... In a series Solve the inverse problem under the given conditions and plot the results. and The relationship curve is selected based on the point corresponding to the maximum curvature. .

[0240] The complete Reg-SIRT algorithm flow is demonstrated by inputting measured time steps. Sensitivity matrix Initial slowness Regularization parameters Maximum number of iterations Convergence threshold (RMS residual). Cycle This process continues until convergence or the maximum number of iterations is reached. The inversion iteration process is as follows:

[0241] The first step is forward modeling.

[0242]

[0243] The theoretical time travel is calculated based on the initial or iterative input values.

[0244] The second step is to calculate the residuals.

[0245] Calculate measured time Compared with theoretical calculations residual :

[0246]

[0247] In the formula, The measured internal travel time of the rock obtained in step S3 .

[0248] Step 3: SIRT Update

[0249] Joint Iterative Reconstruction Method (SIRT) is a classic algorithm for solving the inverse problem of tomography. It has the advantages of convergence stability and insensitivity to noise. Based on the standard SIRT, this invention proposes Regularized SIRT (Reg-SIRT) by combining Tikhonov regularization.

[0250] For the In the next iteration, the temporary slow speed is updated to... :

[0251]

[0252]

[0253] In the formula, To cross the first The total number of rays in each grid;

[0254] The fourth step is regularization smoothing.

[0255] This invention employs explicit regularization to impose smoothing constraints on the solution, modifies the iterative formula by adding a regularization gradient term, and transforms the regularization problem into an iterative form: in each iteration, the solution is first updated according to SIRT. Then solve for a smoothing step:

[0256]

[0257] The problem has an explicit solution:

[0258]

[0259] Constructing sparse matrices :

[0260]

[0261] Construct the solution equation:

[0262]

[0263] because It is a sparse symmetric positive definite matrix, which can be solved quickly using the conjugate gradient method without explicit inversion. The CG solver is called once after each SIRT iteration to obtain the regularized slowness. .

[0264] Fifth step, convergence judgment criteria

[0265] Calculate the current residual :

[0266]

[0267] like ( )or Or reach the maximum number of iterations If the iteration stops, then stop.

[0268] Step 6: Update the iteration counter

[0269]

[0270] After the final iteration, the final slowness vector is obtained. The wave velocity within each grid is:

[0271]

[0272] (5) Visualization of wave velocity field

[0273] The wave velocity values ​​obtained from the inversion are mapped onto the grid to create a pseudo-color cloud map.

[0274] The inverted wave velocity values ​​were plotted as pseudo-color cloud maps using MATLAB's matplotlib library. The horizontal axis represents the width of the rock specimen, and the vertical axis represents the height. The colors gradually change from blue (low wave velocity) to red (high wave velocity), thus obtaining the spatial distribution map of the wave velocity inside the rock under the current stress level.

[0275] The stress level (MPa), axial strain (%), loading time (s), and wave velocity scale are labeled on the image. At the same time, the numerical data are output: the grid center coordinates and the corresponding inverted wave velocity values, and saved as a CSV file for subsequent index calculations.

[0276] Step S5: Multi-index rock mass disaster monitoring and early warning based on temporal anomalies

[0277] This step utilizes the wave velocity field image sequence obtained from the continuous inversion in the previous step to construct a spatiotemporal sequence prediction model based on Convolutional Long Short-Term Memory Network (ConvLSTM), thereby enabling advanced prediction of the wave velocity field evolution trend before rock mass failure, and performing multi-indicator fusion early warning based on the prediction results.

[0278] (1) Construction of wave velocity field image sequence dataset

[0279] During the rock loading process, steps S2 to S4 are repeated every 30 seconds to continuously acquire a series of wave velocity field distribution images, forming a dynamic evolution sequence of the wave velocity field with loading time.

[0280] Collected The wave velocity distribution images are denoted as a sequence. , For pixel height, The width is pixels, and each frame of the image is [pixel width]. for The matrix, where each element represents the inverted wave velocity value of the corresponding grid cell. .

[0281] To eliminate the influence of lithological differences on absolute wave velocity values, each image frame is normalized:

[0282]

[0283] In the formula, and These are the maximum and minimum values ​​of the wave velocity of all pixels in the entire sequence, respectively.

[0284] Training samples are generated using the sliding window method: with a length of... The consecutive frames are used as the input sequence, with the length of the next frame being... The continuous frames are used as the prediction target. Assuming a sliding step size of 1, the sample set can be represented as:

[0285]

[0286] in:

[0287]

[0288]

[0289] (2) ConvLSTM Spatiotemporal Prediction Network Architecture

[0290] This invention uses ConvLSTM as the core prediction module. This network organically combines the spatial feature extraction capability of convolutional neural networks (CNN) with the temporal modeling capability of long short-term memory networks (LSTM), which can effectively capture the local damage features and long-range dependencies in the evolution of wave velocity fields.

[0291] The network adopts an encoder-decoder structure, with the following architecture:

[0292] 1) Encoder

[0293] The encoder is made by It consists of stacked ConvLSTM layers, extracting spatiotemporal features of the input sequence layer by layer. For the th... layer( The calculation process is as follows:

[0294]

[0295] In the formula: This represents the convolution operation. Represents element-wise product; Use the Sigmoid activation function; It is the hyperbolic tangent activation function; For input gate, For the Gate of Oblivion For output gate; In cellular state, Candidate cell state, It is in a hidden state; Learnable convolution kernel parameters; This is the bias vector.

[0296] Input to the first layer of ConvLSTM The image frame is a normalized wave velocity field; the input for subsequent layers. The previous layer is in a hidden state. .

[0297] The encoder ultimately outputs the hidden state at the last time step. and cell state , serving as the initial spatiotemporal features of the decoder.

[0298] 2) Decoder

[0299] The decoder is also made by It consists of layers of ConvLSTM, with a structure symmetrical to the encoder. The initial state of the decoder is copied from the final state of the encoder.

[0300] During the prediction phase, the decoder iteratively generates the future, starting from the final state output by the encoder. The prediction result for the first frame One prediction step ( The decoder receives the output of the previous prediction step as the current input, processes it through a ConvLSTM layer, and then maps it to a single-channel wave velocity field prediction image by the final output layer.

[0301]

[0302] 3) Network parameter configuration

[0303] In this embodiment, the number of input frames Output frame count Image size 0, ConvLSTM layer number The number of convolutional kernels in each layer are 64, 128, and 64, respectively, and the kernel size is [missing information]. The step size is 1, and the padding is 1 to keep the feature map size unchanged.

[0304] (3) Loss function and model training

[0305] To simultaneously ensure pixel-level accuracy and structural similarity of the predicted image, this invention designs a combined loss function:

[0306]

[0307] The definitions are as follows:

[0308] 1) Mean squared error loss:

[0309]

[0310] 2) Structural Similarity Loss: The Structural Similarity Index (SSIM) measures the similarity between two images from three dimensions: brightness, contrast, and structure. The SSIM loss is defined as:

[0311]

[0312] 3) Gradient difference loss: The spatial gradient of the wave velocity field reflects the degree of damage localization. Introducing gradient difference loss can enhance the model's ability to predict the evolution of damage boundaries.

[0313]

[0314] In the formula, The spatial gradient is calculated for the Sobel operator. and The balance coefficient (taken in this embodiment) , ).

[0315] Training uses the Adam optimizer with an initial learning rate of The batch size is 8, and the training epochs are 200. An early stopping strategy is used to prevent overfitting: if the loss on the validation set does not decrease for 20 consecutive epochs, training is stopped.

[0316] (4) Multi-level early warning criteria based on predicted wave velocity field

[0317] After training is completed, during the real-time monitoring phase, the most recent Frame velocity field image input model, obtain future Predicted sequence of frames .

[0318] For each frame in the prediction sequence, extract the following warning-sensitive features:

[0319] 1) Predict the average wave velocity attenuation rate Average wave velocity of the predicted frame Perform a linear fit, and its slope Indicates the rate of decrease in wave speed:

[0320]

[0321] 2) Predict the growth rate of the low-wave velocity region area The least squares method is used to analyze the data points. Performing a univariate linear regression, the fitted straight line equation is:

[0322]

[0323] According to the least squares principle, the slope and intercept The calculation formula is:

[0324]

[0325]

[0326]

[0327]

[0328] In the formula: This is the average of the time sequence numbers; for Average of the frame prediction area percentage; slope The unit is % / frame, which represents the average percentage increase in the area of ​​the low beam velocity region after each prediction time step.

[0329] 3) Predicting the spatial entropy of the wave velocity field The spatial entropy of the wave velocity field reflects the degree of non-uniformity in damage distribution, and is defined as:

[0330]

[0331] In the formula, the wave velocity value is divided into equal parts. Each interval For wave speed values ​​falling into the first The pixel ratio of each interval. The sharp increase indicates the formation of localized damage.

[0332] 4) Predicting the probability of destruction : Predict the frame , and Input a lightweight fully connected classification network and output the future... The probability of macroscopic destruction occurring within a given time period .

[0333] The criteria for determining the warning level are as follows:

[0334] Yellow alert: or < ( (The value is the historical average minus two standard deviations), indicating that the damage has entered an accelerated phase.

[0335] Orange alert: and ( (The value is the historical average minus two standard deviations), indicating that damage localization has been significantly formed.

[0336] Red Alert: Furthermore, the predicted low-velocity region has already connected into a continuous strip, suggesting that a macroscopic rupture is imminent.

[0337] Step S6: Detection of rock hazard zones based on spatial distribution patterns

[0338] This step transforms the historical wave velocity field sequence obtained from step S4 and the future wave velocity field sequence predicted from step S5 into a damage variable field sequence. From the perspective of damage mechanics, it dynamically tracks the entire process of the initiation, expansion, localization and penetration of damage zones inside the rock mass. Combined with the prediction capabilities of ConvLSTM, it achieves advanced prediction and accurate location of future disaster areas.

[0339] (1) Quantitative conversion of wave velocity field to damage variable field

[0340] The decrease in wave velocity is essentially a stiffness degradation caused by the initiation and propagation of microcracks within the rock mass. To establish a direct link between wave velocity changes and damage mechanics, this invention defines a wave velocity damage variable. :

[0341]

[0342] In the formula: For a moment In position The inversion wave velocity value at that location; The wave velocity value at this location in the initial lossless state is taken from the first inversion result before loading begins; , Indicates a lossless state. Indicates complete damage.

[0343] The wave velocity fields for all historical moments obtained in step four, and the future moments predicted in step five. The wave velocity field is converted into the corresponding damage variable field sequence:

[0344]

[0345]

[0346] (2) Identification of damage localization zones and multi-frame correlation tracking

[0347] During the damage evolution process, the aggregation and penetration of microcracks will form a spatially continuous damage localization zone. Its geometric characteristics and evolution trajectory are the key basis for determining the final failure location.

[0348] 1) Single-frame identification of damage localization zones

[0349] For each frame of damage field Perform the following identification process:

[0350] Binarization: Setting a damage threshold Generate a binary mask:

[0351] (1)

[0352] Morphological filtering: 3×3 median filtering is used to remove isolated single-pixel noise, and then morphological closing operation is used to fill small holes to make the damaged area continuous.

[0353] Connectivity labeling: All connected regions are labeled using the 8-neighborhood connectivity criterion. Remove areas smaller than the threshold In small areas, the main localized damage zone is preserved.

[0354] 2) Extraction of geometric feature parameters of the damaged area

[0355] For each effective damage area Extract the following quantitative geometric features:

[0356] ① Area $: The total number of valid grids in the region multiplied by the area of ​​a single grid.

[0357] ② Equivalent center coordinates , :

[0358] (2)

[0359] (3)

[0360] In the formula, For the region The number of pixels within.

[0361] ③ Major axis length With minor axis length The calculation process is as follows: Construct a covariance matrix for the pixel coordinates within each damaged area. :

[0362] (4)

[0363] Eigenvalue decomposition is performed on the covariance matrix to measure the dispersion and directionality of pixels within the damaged region relative to the center point, and the results are obtained. eigenvalues , And the corresponding eigenvectors, the eigenvector corresponding to the largest eigenvalue represents the main extension direction of the damaged area, and the eigenvector corresponding to the smallest eigenvalue represents the secondary direction orthogonal to the main extension direction;

[0364] Based on two eigenvalues and The equivalent major axis length and equivalent minor axis length of the damaged region are calculated separately to obtain the morphological scale characteristics of the damaged region:

[0365] , (5)

[0366] ④ Principal direction angle : Major axis direction and The angle in the positive direction of the axis (counterclockwise is positive).

[0367] ⑤ Average damage value :

[0368] (6)

[0369] 3) Temporal correlation and evolution trajectory construction of multi-frame damage regions

[0370] To track the evolution of the same damaged area, it is necessary to match and correlate damaged areas at adjacent time points. A joint matching criterion based on location and area similarity is adopted:

[0371] For time area and time area Define the matching cost function:

[0372]

[0373] In the formula, Where is the radius of the rock specimen. The weighting coefficient (in this embodiment, we take...) =0.6). The Hungarian algorithm is used to solve for the optimal match, and the successfully matched regions are assigned the same tracking ID to construct the evolution trajectory of the damaged region:

[0374]

[0375] (3) Quantitative division of damage evolution stages and identification of damage precursors

[0376] Based on the evolution trajectory of the damaged area, three key stages of damage evolution can be quantitatively divided:

[0377] 1) Damage initiation stage

[0378] Feature: First appearance satisfies A continuous region with area S < 0.1 ( (The cross-sectional area of ​​the specimen), the center position is relatively stable, and the average damage value is... Slow growth. This stage corresponds to the linear elastic segment of the stress-strain curve.

[0379] 2) Damage Stabilization and Propagation Phase

[0380] Characteristics: The area S grows at an approximately constant rate, and the principal direction angle... It gradually tends to form a specific angle with the direction of the maximum principal stress (consistent with the Mohr-Coulomb fracture criterion). Average damage value. It continues to increase.

[0381] Define the area growth rate:

[0382]

[0383] when When the time interval is positive for N consecutive time intervals and the variance is small, it is determined to be in the stable expansion phase.

[0384] 3) Accelerated localization of damage (critical stage of precursory disaster)

[0385] Feature: The area S exhibits exponentially accelerating growth. Rapidly increasing; damage localization index A significant increase, defined as the ratio of the area of ​​the largest damaged region to the total area of ​​all damaged regions:

[0386]

[0387] This indicates that the damage is highly concentrated in a single main rupture zone, which is a direct precursor to impending macroscopic destruction.

[0388] In addition, the aspect ratio of the damaged area is defined. :

[0389]

[0390] when A rapid increase indicates that the damage zone is extending rapidly and tending to become continuous.

[0391] Comprehensive criteria for identifying precursors to destruction: When the following conditions are met simultaneously, the rock mass is determined to have entered a critical state of catastrophic event:

[0392] , To stabilize the average rate during the expansion phase, It is the instantaneous rate of change at the current moment; >0.8; >3.0 and principal direction angle The deviation from the predicted fracture direction is less than 15°.

[0393] (4) Advance prediction of future danger zones based on ConvLSTM prediction results

[0394] The ConvLSTM network output in step S5 yields the future... The predicted wave velocity field sequence of the frame. This invention further transforms this predicted sequence into a predicted damage field sequence. The aforementioned damage identification and tracking algorithms are then applied to the prediction frame to achieve advanced prediction of future danger zones.

[0395] 1) Identification and feature extraction of predicted damage areas

[0396] For each predicted frame The geometric features of the damaged area are extracted according to equations (1)-(6) to obtain the predicted center location. ,area Main direction angle .

[0397] 2) Prediction of future fracture surface location

[0398] Based on the predicted evolution trend of the damage center, a quadratic polynomial is used to fit and predict the center trajectory:

[0399]

[0400] Extrapolation This yields the expected location of the fracture surface.

[0401] 3) Spatial probability distribution of future low-velocity anomaly regions

[0402] Multiple predicted trajectories are generated using multiple random forward propagations of a ConvLSTM network, and the predicted damage area for each spatial location at a future time is calculated. Construct a catastrophe probability graph based on the probability of ( ).

[0403]

[0404] In the formula, The number of random forward propagations (in this embodiment, we take...) ), For the first Predicted damage field from subsampling.

[0405] Probability values ​​in a disaster probability graph The area was designated as a high-risk disaster zone. It is a medium-risk area.

[0406] (5) Automatic identification of destruction mode

[0407] Based on the evolution characteristics of the damaged area, this invention achieves automatic classification of rock mass failure modes:

[0408] 1) Shear failure mode

[0409] Characteristics: The localized lesion appears as a single, continuous band with a length-to-width ratio of [missing information]. Main direction angle When the angle between the direction of the maximum principal stress and the direction of the maximum principal stress is in the range of 20°-40°, the area growth exhibits a typical three-stage characteristic of "slow → accelerated → rapid penetration".

[0410] 2) Tension splitting failure mode

[0411] Characteristics: The damaged area is distributed in multiple stripes approximately parallel to the principal stress direction, with the principal stress angle... 0° or 180°, with a large length-to-width ratio but the strips are independent of each other. The value is relatively low ( <0.7).

[0412] 3) Hybrid destruction mode

[0413] Characteristics: It has the morphological characteristics of both of the above modes. The damaged area has both shear bands that are oblique to the principal stress and tensile cracks that are parallel to the principal stress direction, and the two intersect or are close to each other in space.

[0414] The failure mode identification results are marked with text labels on the output report to provide a basis for subsequent support design.

[0415] (6) Visual location and report output of hazardous areas

[0416] Based on the above analysis, the following information will be integrated and visualized:

[0417] 1) Current damage state diagram: The damage field at the current moment is displayed on the cross-sectional profile of the rock specimen using a pseudo-color overlay. And marked with red contour lines. The boundary of the damaged area is marked with an arrow indicating the evolution trajectory of the damage center.

[0418] 2) Catastrophic probability prediction map: displayed in the form of a heat map. The medium-risk area is marked with a yellow dashed circle, and the high-risk disaster area is marked with a red solid circle. The predicted damage mode is also marked in the legend.

[0419] 3) Quantitative characteristic curve: Outputs the area of ​​the main damage region. Aspect Ratio Localization index The curve shows the change over time and indicates the current stage of damage evolution.

[0420] 4) Spatial positioning data table: Outputs the geometric center coordinates, equivalent radius, major axis direction angle and predicted arrival time of high-risk disaster areas in tabular form.

[0421] The visualization results and location data are saved in report form, providing clear spatial target areas and time windows for on-site pressure relief, anchor bolt support, or grouting reinforcement. Through the above steps, this invention achieves a complete process of rock mass loading failure disaster early warning and characterization, from non-contact data acquisition, wave velocity field inversion, time-series early warning to hazard zone location.

[0422] The present invention and its embodiments have been described above. This description is not restrictive. The accompanying drawings are only one embodiment of the present invention. In short, if a person skilled in the art is inspired by this description and designs a similar structure and embodiment without departing from the spirit of the present invention, such design should fall within the scope of protection of the invention.

Claims

1. A non-contact wave velocity testing and inversion characterization method for early warning of rock mass loading failure disasters, characterized in that, Includes the following steps: S1. Construct a non-contact ultrasonic testing experimental system: After waterproofing the rock specimen, immerse it in a water tank, using water as the coupling medium; arrange an ultrasonic transducer array in a ring on the outer wall of the water tank to form a non-contact sound wave transmission path; connect the transducer array to an ultrasonic pulse transmitter and receiver and a data processing terminal to establish a synchronous triggering and multi-channel parallel acquisition link. S2. Multi-angle acoustic travel time data acquisition: Each transmitting transducer is excited in sequence, and all receiving transducers synchronously acquire the transmitted waveform to form a multi-angle intersecting ray network covering the interior of the rock mass and pick up the first arrival travel time of each waveform. S3. Iterative solution and travel time correction of the refraction path at the water-rock interface: In order to eliminate the systematic travel time error caused by strong refraction at the water-rock interface, an iterative solution method for the refraction path based on Snell's law and Fermat's principle is established to obtain the true travel time inside the rock. S4. Wave velocity field tomographic inversion based on regularized joint iterative reconstruction: The rock mass section is discretized into a grid model, and a linear forward modeling equation for travel time and slowness is established; the actual travel time inside the rock is used as the observation value, and the slowness is updated by joint iterative reconstruction; Tikhonov regularization and first-order difference smoothing constraints are introduced in the iteration process, and spatial smoothing constraints are applied to the slowness field through explicit regularization steps; after the iteration converges, the spatial distribution of wave velocity inside the rock mass is obtained and visualized in the form of pseudo-color cloud map. S5. Multi-indicator fusion disaster early warning based on spatiotemporal sequence deep learning prediction: Repeat steps S2 to S4 during the loading process to obtain a dynamic evolution image sequence of the wave velocity field; construct a spatiotemporal sequence prediction model based on a convolutional long short-term memory network, using historical wave velocity field image frames as input, to predict the evolution trend of the wave velocity field in future multiple frames; based on the predicted wave velocity field sequence, extract multi-dimensional early warning indicators such as average wave velocity attenuation rate, low wave velocity area growth rate, and wave velocity field spatial entropy, and establish multi-level early warning criteria based on the combination of indicator thresholds to achieve advanced identification of precursors to rock mass instability; S6. Spatial localization and damage mode identification of hazard zones based on damage tensor evolution: The wave velocity field obtained from inversion and prediction is converted into a second-order damage tensor field; ellipsoid fitting is performed based on the wave velocity square ratio of ray clusters in different directions to obtain the principal values ​​and principal directions of the damage tensor; connected component labeling and multi-frame correlation tracking are performed on the damage field sequence to extract the geometric evolution features of the damage region; based on the change law of damage features, the three stages of damage initiation, stable expansion and accelerated localization are quantitatively divided; the damage mode is automatically identified based on the principal value ratio and principal direction features of the damage tensor; and a hazard probability map is generated by multiple random forward propagation of ConvLSTM to delineate high and medium risk areas and output hazard zone information.

2. The non-contact wave velocity testing and inversion characterization method for early warning of rock mass loading failure disasters according to claim 1, characterized in that, In step S2, the initial arrival time is automatically picked up using a combination of the long and short time window ratio method and the Akaike information criterion, supplemented by manual correction.

3. The non-contact wave velocity testing and inversion characterization method for early warning of rock mass loading failure disasters according to claim 1, characterized in that, In step S3, the iterative solution method for the refraction path includes: using the locations of the water-rock incident point and the exit point as variables, constructing a three-segment propagation path model that includes a straight segment in the water and a curved segment within the rock; using the shortest path method based on Dijkstra's algorithm to trace the curved ray within the rock; using the interface Snell's law residual as the objective function, using the Newton-Raphson method to iteratively solve for the coordinates of the incident point and the exit point, and obtaining the true sound wave travel time inside the rock after deducting the propagation time in the water medium, which serves as the observation data for wave velocity field inversion.

4. The non-contact wave velocity testing and inversion characterization method for early warning of rock mass loading failure disasters according to claim 3, characterized in that, In step S3, the complete sound wave propagation path is divided into three segments: the incident segment in water, the transmission segment inside the rock, and the exit segment in water. The method for establishing the three-segment refraction path model is as follows: A horizontal section at the midpoint of the rock specimen's height is selected as the analysis plane. A two-dimensional rectangular coordinate system is established with the center of the specimen's cross-section as the origin. The radius of the rock specimen is set as... The boundary equation of the cross section is Let the coordinates of the transmitting probe be... The coordinates of the receiving probe are The wave velocity in the water medium is The current wave velocity field inside the rock is ; The unknown quantity to be solved is the incident point of the sound wave at the water-rock interface. and the ejection point at the rock-water interface Both points are located on the circumference of the specimen, satisfying the condition that... , ; Given and After setting the coordinates, calculate the travel time in the water for the two segments: In the formula, For the travel time of the sound wave incident in water, The travel time of the sound wave emitted from the water segment.

5. The non-contact wave velocity testing and inversion characterization method for early warning of rock mass loading failure disasters according to claim 3, characterized in that, Step S3 uses the shortest path method combined with Dijkstra's algorithm to solve for the curved rays within the rock, specifically including: The wave velocity field inside the rock is non-uniformly distributed, and the sound wave propagation satisfies the equation of process function. The rock cross section is discretized into a square grid and each grid is assigned a slowness. The grid node topology is constructed, and the edge propagation travel time is defined by the node distance and the average slowness of the grid. The Dijkstra algorithm is run with the water-rock incident point as the source point. The grid nodes are traversed and the shortest propagation time is updated. The calculation is terminated early when the corresponding exit node is reached. The discrete ray path nodes are obtained by backtracking in reverse. A continuous and smooth acoustic curve within the rock is obtained by using cubic spline interpolation to smooth the discrete node sequence. The interpolation process satisfies the conditions of coordinate continuity, second derivative continuity and natural boundary constraints.

6. The non-contact wave velocity testing and inversion characterization method for early warning of rock mass loading failure disasters according to claim 1, characterized in that, In step S4, the slowness update formula of the joint iterative reconstruction method is: in, The total number of rays passing through the j-th grid. Let be the path length of the i-th ray within the j-th grid.

7. The non-contact wave velocity testing and inversion characterization method for early warning of rock mass loading failure disasters according to claim 6, characterized in that, In step S4, the explicit regularization smoothing step is as follows: The regularization problem is transformed into an iterative form: in each iteration, the result is first updated using the joint iterative reconstruction method. Then solve for a smoothing step: The problem has an explicit solution: Constructing sparse matrices : Construct the solution equation: It is a sparse symmetric positive definite matrix, which can be solved quickly using the conjugate gradient method without explicit inversion. The conjugate gradient method solver is called once after each iteration of the joint iterative reconstruction method to obtain the regularized slowness. In the above formula, λ is the regularization parameter, and L is the first-order difference regularization matrix.

8. The non-contact wave velocity testing and inversion characterization method for early warning of rock mass loading failure disasters according to claim 1, characterized in that, In step S5, a spatiotemporal sequence prediction model based on a convolutional long short-term memory network is constructed, wherein the convolutional long short-term memory network adopts a symmetrical structure composed of an encoder and a decoder. The encoder is composed of multiple stacked ConvLSTM layers, which extract the spatiotemporal features of the input wave velocity field sequence layer by layer; The first ConvLSTM layer takes a normalized wave velocity field image frame as input, and the subsequent layers take the hidden state output from the previous layer as input. The encoder outputs the hidden state and cell state of the last time step as the initial spatiotemporal features of the decoder. The decoder stacks the same number of ConvLSTM layers, and its structure is symmetrical with that of the encoder. The initial state directly copies the output features of the encoder. In the prediction phase, the decoder iterates to generate future multi-frame wave velocity fields. The prediction results of each step are fed into the next layer of the network, and finally a single-channel wave velocity field prediction image is obtained through mapping. During offline training, the model employs a composite loss function that integrates mean squared error, image structure similarity, and gradient difference. The Adam optimizer is used for iterative convergence, and the trained network is used for real-time wave velocity field advance prediction.

9. The non-contact wave velocity testing and inversion characterization method for early warning of rock mass loading failure disasters according to claim 1, characterized in that, In step S6, the wave velocity field is converted into a damage variable field, including: Define wave velocity damage variable : In the formula: For a moment In position The inversion wave velocity value at that location; This represents the wave velocity at that location in its initial, undamaged state.

Images