TBM crossing fault fracture zone disturbance area rock instability dense array acoustic wave dynamic detection method

Abstract

The present application belongs to the field of tunnel advanced geological prediction, and provides a TBM passing through fault fracture zone disturbance area rock mass instability intensive array acoustic dynamic detection method, obtains multi-period observation data, pre-processes each period of observation data, carries out time picking on the pre-processed data, and obtains multi-period travel time; based on the multi-period travel time, carries out benchmark inversion and monitoring inversion, the benchmark inversion obtains the benchmark wave velocity inversion result by iteratively fitting the observation data time, which is input into the monitoring inversion as the initial model, the monitoring inversion reconstructs the TBM disturbance area wave velocity distribution at different time by fitting the double difference of the time data before and after the two periods, and obtains the monitoring wave velocity inversion result; the benchmark wave velocity inversion result and the monitoring wave velocity inversion result are used as the background model, and the reverse time migration imaging method is used for dynamic imaging of the TBM excavation disturbance area development process. The present application realizes the dynamic detection and description of the TBM excavation disturbance area development range and its development state.

Description

Technical Field

[0001] This invention belongs to the field of tunnel advanced geological prediction, specifically involving a method for dynamic detection of rock mass instability by dense array acoustic waves when a TBM traverses a fault fracture zone disturbance area. Background Technology

[0002] The statements in this section are merely background information related to the present invention and do not necessarily constitute prior art.

[0003] The excavation and unloading of rock mass during TBM tunnel construction disturbs the in-situ stress field, leading to stress redistribution. Because energy is easily released into the tunnel cavities, the surrounding rock mass fractures. As the depth away from the tunnel walls increases, the rock mass's bearing capacity increases, and the degree of disturbance weakens, forming a disturbed rock mass zone encompassing fractured, plastic, and elastic zones. Since stress redistribution is a slow process, the rock fractures within the disturbed zone can further expand, potentially leading to rockfalls and collapses.

[0004] Existing methods for detecting disturbance zones focus on predicting the loosening zone of the tunnel arch after excavation, and can be divided into active detection methods and passive monitoring methods.

[0005] Active detection methods primarily utilize geophysical fields to observe disturbed areas. By processing and inverting the observation data into images, the extent of the disturbed area can be delineated. However, the application of active detection methods in TBM tunnels faces several challenges: First, the massive metal frame and electrical equipment of the TBM cause severe electromagnetic interference to resistivity and ground-penetrating radar methods, leading to effective signal attenuation and abnormal response distortion, making it difficult to guarantee the stability and reliability of the detection results. Second, acoustic testing and cross-hole seismic methods rely on drilling, resulting in long operation cycles and high economic costs, which cannot meet the needs of rapid TBM construction.

[0006] Passive monitoring methods primarily infer the condition of the disturbance zone in the tunnel arch through monitoring data, mainly including monitoring and measurement methods and microseismic monitoring methods. Monitoring and measurement methods generally use multi-point displacement gauges to monitor the convergence deformation of the surrounding rock in real time and determine the development status of the disturbance zone. However, because it only uses apparent convergence deformation data, it is difficult to effectively identify the spatial extent of the disturbance zone and usually needs to be used in conjunction with other methods. Microseismic monitoring or acoustic emission methods receive wave signals released by the fracturing of the surrounding rock through acoustic sensors, and combine this with the source location results to invert the location of the rock mass fracturing, thereby analyzing the development degree and spatial extent of the disturbance zone. In field applications, monitoring and measurement and microseismic monitoring are often used in combination to fully obtain the development characteristics of the disturbance zone. However, these methods usually rely on the already occurring response information such as surrounding rock deformation, stress anomalies, and microseismic activity to identify the rock mass disturbance state. Their monitoring results have a significant lag, making it difficult to effectively predict potential rock mass disturbance areas and their development trends before a disaster occurs, and thus failing to effectively support the prevention and pre-treatment of TBM jamming disasters.

[0007] In summary, existing methods are not effectively applicable to TBM-constructed tunnels, and they mainly detect the disturbance of the surrounding rock at the top of the excavated tunnel arch, which is not effective in reflecting the dynamic changes in the rock mass in front of the TBM cutterhead. As a result, the lack of geological information to guide TBM jamming disasters leads to false alarms and missed alarms. Summary of the Invention

[0008] To address the aforementioned problems, this invention proposes a dynamic detection method using dense array acoustic waves to detect rock mass instability in the disturbed zone of a TBM crossing a fault fracture zone. This invention enables dynamic detection and characterization of the development range and state of the TBM excavation disturbance zone.

[0009] According to some embodiments, the present invention adopts the following technical solution: A method for dynamic detection of rock mass instability in a fault fracture zone using dense array acoustic waves via TBM includes the following steps: Acquire observation data from multiple periods, preprocess the observation data from each period, and then pick up the time of arrival of the preprocessed data to obtain the travel times for multiple periods; Based on the multi-period travel times, benchmark inversion and monitoring inversion are performed. The benchmark inversion obtains the benchmark wave velocity inversion result by iteratively fitting the arrival time of the observation data, and uses it as the initial model input into the monitoring inversion. The monitoring inversion reconstructs the wave velocity distribution of the TBM disturbance area at different times by fitting the double difference of the arrival time data of the two periods before and after, and obtains the monitoring wave velocity inversion result. Using the reference wave velocity inversion results and the monitoring wave velocity inversion results as background models, the reverse time migration imaging method is used to perform dynamic imaging of the development process of the TBM excavation disturbance zone.

[0010] As an alternative implementation method, the preprocessing of observation data for each period includes: estimating the tilt information and structural coherence of acoustic signals in the observation data using the multi-scale structure tensor method, setting a coherence threshold using the correlation between local coherence and acoustic signals, and dividing the acoustic signals into signal regions and noise regions. A hybrid total variational denoising random noise suppression method based on structural tilt angle is adopted to adaptively denoise the signal region and the noise region, thereby suppressing the background noise. The hybrid total variational denoising random noise suppression method based on structural tilt angle guidance adaptively denoises the signal region and the noise region, including: For the signal region, where the tilt information is stable, a total variation denoising method guided by the data tilt angle is used for noise suppression. For the noise region, where there is no obvious tilt information, a total variation denoising method is directly used for noise suppression, thereby effectively balancing signal leakage and noise suppression.

[0011] As a further defined implementation, the tilt angle information and structural coherence of the acoustic signal are estimated using the multi-scale structure tensor method. The multi-scale structure tensor method quantifies the local structural features of a pixel by analyzing the gradient in its neighborhood, including: First, the acoustic data is filtered using a derivative Gaussian window function. Then, a ring filter is used to filter the structure tensor. This filter can align with the feature edges and accurately normalize the edge energy to the feature center, thus obtaining the filtered set of structure tensors. The filtered structure tensors at different scales are normalized, and the trace of the structure tensor at each scale is solved at each acoustic data location. The optimal filter scale at each pixel is obtained through one-dimensional maximum search. By extracting the structure tensor corresponding to the optimal filter scale pixel by pixel, the optimal structure tensor is obtained. The eigenvalues ​​of the optimal structure tensor are then solved to obtain the signal structure tilt angle information. The coherence of the local structure is calculated using the eigenvalues ​​of the structure tensor, yielding the results of the local structure coherence.

[0012] As a further defined implementation, the tilt angle information and structural coherence of the acoustic signal are estimated using the multi-scale structure tensor method. The multi-scale structure tensor method quantifies the local structural features of a pixel by analyzing the gradient in its neighborhood. First, the acoustic data is filtered using a derivative Gaussian window function. Then, a ring filter is used to filter the structure tensor. This filter can align with the feature edges and accurately normalize the edge energy to the feature center, thus obtaining the filtered set of structure tensors. The filtered structure tensors at different scales are normalized, and the trace of the structure tensor at each scale is solved at each acoustic data location. The optimal filter scale at each pixel is obtained through one-dimensional maximum search. By extracting the structure tensor corresponding to the optimal filter scale pixel by pixel, the optimal structure tensor is obtained. The eigenvalues ​​of the optimal structure tensor are then solved to obtain the signal structure tilt angle information. The coherence of the local structure is calculated using the eigenvalues ​​of the structure tensor, yielding the results of the local structure coherence.

[0013] As a further defined implementation method, the preprocessing process for observation data in each period also includes: obtaining a clean effective signal based on denoising preprocessing, selecting multiple wave data in the effective signal by pre-windowing strategy, extracting the primary reflection wave component by using generalized singular value decomposition and low-rank approximation, reconstructing the primary reflection wave component in the effective signal, and obtaining reconstructed data. The local correlation method is used to calculate the local similarity between the reconstructed data and the effective signal. A similarity threshold is set, a binary mask of the primary reflection wave data is constructed and applied to the original clean signal. At the same time, the reconstructed noise and the spatial location information of the residual multiple reflection waves are removed, so as to achieve effective separation of the primary reflection wave and the multiple waves.

[0014] As a further defined implementation, multiple wave data in the effective signal are selected by pre-windowing strategy, and the primary reflection wave component is extracted using generalized singular value decomposition and low-rank approximation. The primary reflection wave component in the effective signal is then reconstructed to obtain reconstructed data, including: First, the multiple wave data in the effective signal are selected by box selection, and the effective signal is divided into a region where the primary and multiple reflection waves overlap and a region where multiple reflection waves overlap, and data matrices for the two regions are constructed. The two data matrices are decomposed using the generalized singular value decomposition method to obtain the eigenvector matrix and the diagonal matrix, respectively. By comparing the magnitudes of the singular values ​​in the diagonal matrix, the dominant component of the primary reflected wave is obtained based on the comparison results, and the reconstructed data is obtained. By employing the local correlation method, the local similarity between the reconstructed data and the effective signal is calculated. A similarity threshold is set, and a binary mask for the primary reflection wave data is constructed and applied to the reconstructed data. This effectively suppresses reconstruction noise and residual multiple reflection waves, thereby achieving effective separation of the primary reflection wave and multiple reflection waves.

[0015] As an alternative implementation method, the process of acquiring multiple periods of observation data and picking up the arrival time of each period of observation data includes: using the kurtosis characteristics of the transmitted signal of the observation data as the arrival time indicator, and constructing a path cost function based on the kurtosis characteristics of the observation data; Apply hard constraints to the arrival candidate points between adjacent channels; Under the hard constraints, based on the path cost function, the problem of picking up the arrival time of the reflected wave is regarded as a decision problem of picking up the wave one by one along the phase axis of the reflected wave. The result of picking up the arrival time of the reflected wave is regarded as an arrival curve. The decision problem is solved recursively by dynamic programming to obtain the preliminary automatic picking result. The initial automatic picking results are fine-tuned based on the reliability to obtain the final automatic picking results.

[0016] As a further defined implementation, the process of constructing the path cost function includes: Path cost function ; Where, d K + Indicates the kurtosis characteristics of the observed data. λ As a scaling factor, τ i This indicates the pickup result when the reflected wave corresponding to channel i arrives. Let x represent time, x represent the signal, and the first term be the confidence cost function, calculated by observing the kurtosis characteristic d of the data. K + The first term characterizes the reliability of the candidate point as the arrival time of the reflected wave. The higher the energy, the higher the probability of it being the correct arrival time of the reflected wave, and the smaller the loss. The second term is a smoothing function, which describes whether the travel time changes between two adjacent candidate points are smooth.

[0017] As a further defined implementation, the process of solving the decision problem using dynamic programming recursion includes: using dynamic programming recursion, we obtain: ; in, Let Δ be the cost function of the optimal path. t max For the maximum time interval, τ i d represents the pickup result when the reflected wave corresponding to channel i arrives. K + Indicates the kurtosis characteristics of the observed data. λ This is a scaling factor; it records the optimal index for each path and obtains the entire path by backtracking.

[0018] As a further defined implementation, the process of fine-tuning the preliminary automatic picking results based on credibility to obtain the final automatic picking results includes: introducing expert experience to correct or rectify the automatic picking results; manually picking multiple arrival points; automatically fine-tuning based on the credibility cost function to ensure that the manually selected points fall at the positions with the strongest kurtosis; and using the manually picked arrival points as anchor points and adding them to the path cost function to obtain: ; Where, d K + Indicates the kurtosis characteristics of the observed data. λ As a scaling factor, τ i This indicates the pickup result when the reflected wave corresponding to channel i arrives. This indicates manually selected anchor points. N t Indicates the sampling point. γ As a scaling factor, when the complexity of the in-phase axis of the reflection data exceeds a predetermined level, γ >0, introduce soft constraints to make the path close to the prior, and realize accurate picking of the arrival time of complex reflected wave phase axis.

[0019] As an alternative implementation method, the process of benchmark inversion and monitoring inversion based on the multi-period travel times includes: referring to the data obtained from a certain detection as the benchmark data arrival time, referring to the subsequent detection data as the monitoring data arrival time, using the benchmark data arrival time as the benchmark state, and using the wave velocity and interface results obtained from its inversion as the background model; extracting the travel time change from the observation data of the subsequent monitoring period to minimize the impact of static structural errors and systematic time delays on the inversion; and solving for the time shift change of the disturbance region based on the difference between the inversion results of the two detection periods.

[0020] As an optional implementation, the process of benchmark inversion based on the multi-period travel times includes: obtaining basic inversion parameters, which include the deployment parameters of the TBM tunnel array acoustic observation system, the arrival time picking results of the direct wave and reflected wave of the observation data, and using the number of in-phase axes of the reflected acoustic waves to evaluate the stratification information of the surrounding rock medium; and establishing an initial model based on geological information, including the interface location and wave velocity distribution. Based on the inversion basic parameters, and combined with the equation, forward modeling models of incident wave travel time field and reflected wave travel time field are constructed respectively, and the simulated travel time of incident wave and simulated travel time of reflected wave are obtained by solving them. A hierarchical joint loss function is constructed that includes both absolute travel time loss term and differential travel time loss term. The average wave velocity update gradient and the B-spline reflection interface control coefficient update gradient are solved according to the hierarchical joint loss function. The average wave velocity of the medium and the control coefficient of the B-spline reflection interface are updated alternately by adopting a layer-by-layer alternating iterative strategy to complete the joint inversion of the wave velocity distribution of the surrounding rock and the geometric morphology of the reflection interface in the TBM excavation disturbance zone.

[0021] As an alternative implementation method, the monitoring and inversion process, which reconstructs the wave velocity distribution of the TBM disturbance zone at different times by fitting the double difference of arrival data from two consecutive periods, includes: using the tunnel reflection wave double difference travel-time shift inversion loss function, specifically: ; in, J f and J r Let these represent the event arrival loss and the receiver double difference loss function, respectively. For the receiving point residual, For the residual at that time, is a coefficient.

[0022] Based on the aforementioned inversion fundamental parameters, and in conjunction with the equation, forward modeling models of the incident wave travel time field and the reflected wave travel time field are constructed respectively, including: The equation is solved by taking the location of the seismic source as the global initial zero-to-time boundary condition, and the global solution yields the incident simulation travel time. Using the incident simulated travel time corresponding to the geological reflection interface as the global initial boundary condition for the reflected wave, the equation is solved again to obtain the global simulated travel time of the reflection. Based on the incident simulation travel time and the reflection simulation travel time, the theoretical arrival times of the corresponding direct waves and layered reflected waves of each detector are extracted, providing theoretical numerical support for the calculation of the loss function.

[0023] As a further defined implementation, the step of solving for the medium average wave velocity update gradient based on the hierarchical joint loss function includes: Based on the adjoint state method, the incident wave specific adjoint field equation and the reflected wave specific adjoint field equation are constructed respectively. Substituting them into the incident simulation travel time and the reflection simulation travel time, the incident adjoint field and the reflection adjoint field are obtained by solving. Combining the mathematical relationship that the slowness of sound waves and the velocity of medium waves are reciprocals, the original gradient of global wave velocity is generated by the global superposition of the incident accompanying field and the reflected accompanying field. The spatial medium partitioning is completed by retrieving the prior information of the surrounding rock medium layering. The original global wave velocity gradients corresponding to all grids within the same medium partition are accumulated and summed to generate the medium average wave velocity update gradient that is uniform across the partition.

[0024] As a further defined implementation, the step of updating the gradient by solving the B-spline reflection interface control coefficients based on the hierarchical joint loss function includes: Gradient transformation calculation is performed based on the wave velocity difference between two adjacent media layers, and the original gradient of the full-space interface geometry is generated from the original global wave velocity gradient. Multi-order, multi-segment B-spline basis functions are used to complete the parametric characterization of the geometry of all reflective interfaces, and a one-to-one mapping mathematical relationship between interface spatial coordinates and B-spline control coefficients is established. The original geometric gradient of the full-space interface is projected along the dimension of the B-spline basis function, and a second-order curvature regularization constraint term to suppress interface overfitting is superimposed to finally generate the updated gradient of the B-spline reflection interface control coefficient.

[0025] As a further defined implementation, the construction of a hierarchical joint loss function that simultaneously includes an absolute travel time loss term and a differential travel time loss term includes: For the observation data corresponding to each layer of reflected wave, a single-layer absolute travel time loss is generated. The absolute travel time loss is the weighted average of the simulated travel time of the reflection and the square of the original observed time difference of the reflected sound wave under the same source and the same detector combination. The differential travel loss is generated in layers. The differential travel loss of the first layer is the weighted average of the simulated travel time of direct incident and the square of the original direct arrival time difference. The differential travel loss of the remaining layers is the weighted average of the simulated travel time of reflection in this layer and the square of the original observed time difference of the reflected sound wave in the previous layer. Set a differential loss growth rate that increases linearly with the iteration rounds. In each iteration, add the corresponding weight to the differential time loss according to the growth rate. The absolute time loss and the weighted differential time loss of the same layer are weighted and fused to generate the total loss of a single layer. The total loss of each single layer corresponding to all reflection layers is combined and summarized to generate the complete joint loss function of the layers.

[0026] As a further defined implementation, the method of alternately updating the average wave velocity of the medium and the B-spline reflection interface control coefficient using a layer-by-layer alternating iterative strategy includes: According to the spatial layering order of the reflective interface from near to far and from shallow to deep, complete iterative calculations are performed on each reflective interface layer separately. Within a single layer, continuous iteration rounds are divided, and differential parameter update operations are performed between odd-numbered and even-numbered iteration rounds. In a single-layer even-numbered iteration cycle, the gradient of the B-spline reflection interface control coefficient is retrieved to correct and update the currently stored B-spline interface control parameters. In a single-layer odd-numbered iteration round, the average wave velocity update gradient of the medium is retrieved, and the currently stored average wave velocity of the layered medium is corrected and updated. After a single layer completes all preset iteration rounds of calculation, the average wave velocity of the medium and the control coefficients of the B-spline reflection interface obtained from the final update of this iteration are all saved as the initial input model for the next layer of reflection interface iteration calculation.

[0027] A TBM-based close-array acoustic dynamic detection system for rock mass instability in a fault fracture zone includes: The preprocessing module is configured to acquire observation data from multiple periods and preprocess the observation data from each period. The time-out pickup module is configured to pick up the preprocessed data at the time-out point to obtain multiple time-out periods; The inversion module is configured to perform baseline inversion and monitoring inversion based on the multi-period travel times. The baseline inversion obtains the baseline wave velocity inversion result by iteratively fitting the arrival time of the observation data, and uses it as the initial model input into the monitoring inversion. The monitoring inversion reconstructs the wave velocity distribution of the TBM disturbance area at different times by fitting the double difference of the arrival time data of the two periods before and after, and obtains the monitoring wave velocity inversion result. The reverse time migration imaging module is configured to use the reference wave velocity evolution results and the monitoring wave velocity inversion results as background models, and to perform dynamic imaging of the TBM excavation disturbance zone development process using the reverse time migration imaging method.

[0028] Compared with the prior art, the beneficial effects of the present invention are as follows: This invention uses a certain period of observation as the baseline state and the wave velocity and interface results obtained from its inversion as the background model; it extracts the travel time variation from the observation data of subsequent monitoring periods, thereby minimizing the impact of static structural errors and systematic time delays on the inversion, and making the inversion results more focused on the rock mass instability changes caused by the TBM excavation and unloading disturbance.

[0029] This invention utilizes a tunnel reflection wave double-difference time-shifting wave velocity inversion method based on the process function equation. It can eliminate the influence of static structural errors and systematic time delays on the inversion by subtracting data, thereby improving the inversion accuracy.

[0030] When picking up the detected wave, this invention uses signal kurtosis as an indicator, which can accurately pick up the arrival time of the reflected wave. At the same time, it takes into account the structural characteristics of the in-phase axis of the reflected wave and adds smoothness constraints to enhance the accuracy of the inter-channel picking results.

[0031] This invention utilizes the multi-scale structure tensor method to solve for the tilt information and structural coherence of acoustic wave data. Then, by leveraging the correlation between structural coherence and the acoustic signal and setting a coherence threshold, the acoustic wave data is divided into signal and noise regions. In the signal region, where tilt information is stable, a tilt-guided total variation denoising method is employed. In the noise region, where there is no significant tilt information, a direct total variation denoising method is used. Through this adaptive constraint denoising method, a balance can be effectively struck between signal fidelity and noise removal.

[0032] This invention, based on the characteristic of late arrival of multiple waves, selects the multiple wave development region, and then uses the generalized singular value decomposition method to project multiple reflections and primary reflections onto the feature space. Using the low-rank approximation, the primary reflection data in the acoustic data profile is reconstructed. Furthermore, by employing the local correlation method, the local similarity between the reconstructed data and the original acoustic data profile is calculated. By setting a similarity threshold, a binary mask for the primary reflection data is constructed, thereby achieving effective extraction of primary reflections and effective suppression of multiple reflections.

[0033] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, preferred embodiments are described below in detail with reference to the accompanying drawings. Attached Figure Description

[0034] The accompanying drawings, which form part of this invention, are used to provide a further understanding of the invention. The illustrative embodiments of the invention and their descriptions are used to explain the invention and do not constitute an improper limitation of the invention.

[0035] Figure 1 This is a schematic diagram of a TBM (Tunnel Boring Machine) method for dynamic detection of rock mass instability in a fault fracture zone disturbance area using dense array acoustic waves, according to one embodiment. Figure 2 This is a schematic diagram of a dynamic imaging process based on time-shifting wave velocity inversion using the process function equation, according to one embodiment. Figure 3(a) is a schematic diagram of a wave speed model of one embodiment, and Figure 3(b) is a schematic diagram of an initial model of one embodiment; Figure 4 This is a schematic diagram of the result of a dynamic imaging method according to one embodiment; Figure 5 This is a schematic diagram of dynamic imaging results obtained by reverse time-shift imaging according to one embodiment. Figure 6(a) shows a fault fracture zone model of one embodiment, and Figure 6(b) shows a schematic diagram of the baseline wave velocity inversion results of one embodiment. Figure 7 This is a schematic diagram of an inversion imaging method according to one embodiment; Figure 8 This is a schematic diagram of the dynamic imaging result of one embodiment; Figure 9 This is a schematic diagram of a control device according to one embodiment. Detailed Implementation

[0036] The present invention will be further described below with reference to the accompanying drawings and embodiments.

[0037] It should be noted that the following detailed description is illustrative and intended to provide further explanation of the invention. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains.

[0038] It should be noted that the terminology used herein is for the purpose of describing particular embodiments only and is not intended to limit the scope of exemplary embodiments according to the invention. As used herein, the singular form is intended to include the plural form as well, unless the context clearly indicates otherwise. Furthermore, it should be understood that when the terms "comprising" and / or "including" are used in this specification, they indicate the presence of features, steps, operations, devices, components, and / or combinations thereof.

[0039] Where there is no conflict, the embodiments and features described in this application may be combined with each other.

[0040] Example 1 A method for dynamic detection of rock mass instability in a fault fracture zone using dense array acoustic waves, such as... Figure 1 As shown, it includes the following steps: Acquire observation data from multiple periods, preprocess the observation data from each period, and then pick up the time of arrival of the preprocessed data to obtain the travel times for multiple periods; Based on the multi-period travel times, benchmark inversion and monitoring inversion are performed. The benchmark inversion obtains the benchmark wave velocity inversion result by iteratively fitting the arrival time of the observation data, and uses it as the initial model input into the monitoring inversion. The monitoring inversion reconstructs the wave velocity distribution of the TBM disturbance area at different times by fitting the double difference of the arrival time data of the two periods before and after, and obtains the monitoring wave velocity inversion result. Using the reference wave velocity inversion results and the monitoring wave velocity inversion results as background models, the reverse time migration imaging method is used to perform dynamic imaging of the development process of the TBM excavation disturbance zone.

[0041] The following is a detailed description.

[0042] First, the principle of time-shift inversion is introduced. It is assumed that a TBM (Tunnel Boring Machine) traverses a fault fracture zone, and the rock mass within the fault fracture zone is disturbed by the TBM excavation, creating a disturbed zone that further develops. Detection work is conducted at two different times, and the data obtained from the first detection is used to... T A,obs When the baseline data is reached, the second probe data will be obtained. T B,obs The time shift change in the model between two detections, referred to as the arrival of monitoring data, is Δ. m It can be represented as: ; in, mCharacterizing the rock mass properties within the target area at a certain moment, in this embodiment, it mainly refers to wave velocity. v As shown in the above equation, given the inversion results of the baseline data, it is only necessary to solve for the time shift change Δ. m This allows us to obtain the rock mass medium conditions during the monitoring period. In numerical simulations, the monitoring period data can be represented as follows: ; Where Φ represents the travel-time forward modeling operator. Based on this, the time shift characteristic is the arrival time difference between two observation periods: ; Time-shift inversion is achieved by fitting two time-shift points. T A,obs , T B,obs Or the time difference between the two periods δT obs The monitoring period model can then be solved. m B or Δ m This enables dynamic inversion of the TBM excavation disturbance zone. Depending on the loss function used, time-shift inversion strategies can be categorized into three types: serial inversion, parallel inversion, and differential inversion.

[0043] The double-difference inversion strategy also uses the reconstructed baseline model as the initial model, and selects the residuals of the two observations as the loss function. The residuals of the observations used can be expressed as: ; The above formula mainly includes two terms: the time-shift residual of the actual observed data and the time-shift residual of the simulated data, hence it is called the two-difference time-shift inversion strategy. The two-difference inversion strategy can be expressed as: ; This method can eliminate the influence of static structural errors and systematic time delays on the inversion by subtracting data, thereby improving the accuracy of the inversion.

[0044] The core idea of ​​the tunnel reflection wave double-difference travel time inversion method proposed in this embodiment is: take a certain period of observation as the baseline state, and use the wave velocity and interface results obtained from its inversion as the background model; extract the travel time change from the observation data of subsequent monitoring periods, thereby minimizing the impact of static structural errors and systematic time delays on the inversion, and making the inversion results more focused on the rock mass instability changes caused by the TBM excavation unloading disturbance.

[0045] Assume a TBM traverses a fault fracture zone, and the rock mass within the fault fracture zone undergoes instability under the disturbance of the TBM. Exploration work is conducted at two different times, and the data obtained from the first exploration is... TA,obs When the baseline data is reached, the second probe data will be obtained. T B,obs When the monitoring data is available, the results of the two-phase detection inversion can be obtained using the following formula: ; The time shift change of the disturbance region can then be solved by subtracting the above equation: ; Through the above process, this paper derives the double-difference loss function for the double-difference travel-time inversion method of tunnel reflected waves, which characterizes that the time shift transformation in the disturbance region can be solved by iteratively fitting the time difference between the two periods. For a specific reflected wave, assume that it is derived from the first... i The earthquake was triggered by the first source and occurred at the second... k A reflected wave is generated at the first reflecting interface and is then reflected by the second. j If multiple detectors receive the signal, the resulting double-difference loss can be expressed as: ; The above equation represents the time residual of a reflection event, which characterizes the changes in the ray path and travel time before and after a disturbance for a given reflection event. Therefore, it is sensitive to changes in wave velocity and interface position, which can lead to coupling crosstalk between wave velocity and interface inversion. To address this, this paper further introduces the receiver double-difference loss, defined as follows: ; in j 0 is the reference detector, and the corresponding residual is: ; The double-difference loss at the receiver points, as characterized by the above formula, emphasizes the travel time differences between different receiver points, which is more conducive to constraining interface morphology or local anomalies. Therefore, the loss function used in the tunnel reflection wave double-difference travel time shift inversion method proposed in this embodiment is: ; in, J f and J r These represent the event arrival time loss and the receiver point double difference loss, respectively. The time-shift inversion proposed in this embodiment requires first using the tunnel array acoustic path function equation wave velocity-interface joint travel-time inversion method to obtain the inversion result of the reference data. This result is then used as the initial model, and the above loss function is used to invert the wave velocity distribution state of the disturbed rock mass. In this embodiment, the time-shift wave velocity inversion can be regarded as a special inversion method, which mainly achieves the inversion of the wave velocity distribution state and its changes in the disturbed rock mass by fitting the travel-time difference of the two reflected waves.

[0046] By utilizing the tunnel reflected wave double-difference travel-time inversion method, this embodiment can obtain the wave velocity distribution of the medium in front of the tunnel at different observation times, and realize the dynamic characterization of wave velocity changes during the instability of rock mass in the TBM disturbance zone. However, the travel-time wave velocity inversion method can only obtain the wave velocity distribution. It needs to be further combined with the reverse time migration (RTM) method to realize the spatial structure information of the TBM disturbance zone and its dynamic change process in detail.

[0047] The acoustic field can be derived from a geological model using a linearized forward modeling operator. L The calculations show that the model can now be decomposed into a background model and a disturbance model: ; Among them, m0 and δ m represents the background model and the disturbance model, respectively. Among them, m0 mainly contains low-frequency background components, which characterize the overall velocity variation trend of the tunnel surrounding rock. δ m mainly contains high-frequency detail components, characterizing the detailed information of the interface structure in the perturbation region. Similarly, the background wavefield and the perturbation wavefield can be written as: ; Where U0 is the background wave field generated by the background model m0, representing the wave field in the ideal background medium; δU is the perturbation model. δ The scattered or reflected wave field caused by m is an effective signal for characterizing the spatial location of the disturbance region. Furthermore, since the total wave field U and model m, and the background wave field U0 and background model m0 all satisfy the above formula, we have: ; Furthermore, under the Born approximation, assuming the perturbation is weaker than the background model, the perturbation wavefield and the perturbation model satisfy an approximately linear relationship: ; The above equation shows that the observed reflected wavefield is essentially the response information of the perturbation model in the background medium. Therefore, as long as a relatively reliable background model is known, the reflection information in the observed data can be shifted and realigned to the true reflection position through wavefield backpropagation and imaging conditions, thereby achieving imaging of the spatial structure of the perturbation region. Therefore, the imaging results of the TBM perturbation region can be solved using the following formula: ; in, δ d represents the reflection data of the interface in the disturbed area, considering L The inverse computation is quite large, therefore it adopts... L T Approximate backpropagation operator. This is achieved by superimposing correlation results from different times and source conditions: ; in,u s and u R These are the forward propagation wavefield and the inverse propagation wavefield, respectively. RTM imaging results can reveal the spatial structure information of the perturbation region, complementing the time-shift velocity inversion results, thereby achieving a detailed characterization of the TBM perturbation region development process.

[0048] The flowchart of the dynamic imaging method based on time-shifting wave velocity inversion using the equation proposed in this embodiment is as follows: Figure 2 As shown, using multi-period observation data as input, the arrival times of each period's reflected wave data are first preprocessed to obtain multi-period travel times. Further benchmark inversion and monitoring inversion are then performed. The benchmark inversion obtains the inversion result by iteratively fitting the arrival times of the observation data. This inversion result can be used as the initial model input for the monitoring inversion. By fitting the double difference between two periods of arrival time data, the wave velocity distribution of the TBM disturbance zone at different times can be effectively reconstructed. Finally, using the inversion result as the background model, the RTM method can be used to achieve dynamic imaging of the TBM excavation disturbance zone development process.

[0049] This paper further establishes a wave velocity model after TBM excavation disturbance based on the wave velocity model shown in Figure 3(a), and verifies the effectiveness of the proposed dynamic imaging method based on the time-shifted reflected wave velocity inversion using the equation. Figure 4 As shown in (a), the fracture network within the fault fracture zone further develops due to the disturbance caused by TBM excavation, and the fault core extends towards the TBM, as... Figure 4 Figure (b) shows the wave velocity disturbance caused by TBM excavation. The time-shifted wave velocity inversion method proposed in this paper is used to invert the wave velocity model after the disturbance and the wave velocity model after the disturbance. The initial model adopts the inversion results shown in Figure 3(b). Figure 4 (c) and (d) respectively show the wave velocity model inversion results after the disturbance and the wave velocity disturbance change inversion results. Figure 4 Figures (e) and (f) show the corresponding interface inversion results. The inversion results can reflect the overall distribution of wave velocity and interface after the disturbance well, and the wave velocity inversion results are basically consistent with the real model. The wave velocity disturbance model inversion results are also relatively consistent with the real model. Similarly, in terms of the accuracy of the interface inversion results, overall, the interface inversion results are basically consistent with the real interface position and shape, with only a certain error in the wave velocity disturbance change inversion results. Figure 4 As shown in (g), the overall wave velocity distribution of the inversion results is consistent with the actual model and can well reflect the time-shift changes of the wave velocity model; similarly, Figure 4 The wave velocity disturbance change inversion results shown in (h) can well show the relative changes in the disturbance range and wave velocity.

[0050] Reverse-time migration imaging was performed based on the wave velocity inversion results before and after the disturbance, and the resulting dynamic imaging results are as follows: Figure 5 As shown, (a) is the imaging result before perturbation, and (b) is the imaging result after perturbation. Overall, the imaging results can reflect the relative position and tendency of the interface. However, due to the limitation of the observation aperture, the reflection information above the interface was not observed by the detector array, thus preventing effective imaging of the upper interface position. This is consistent with the phenomenon in the above inversion results where the slope of the upper interface of the model is difficult to update effectively. Furthermore, the imaging results show that the imaging range of the interface increases with the depth of the reflecting interface, which is due to the enhancement of the illumination range caused by the geometric diffusion of the wave field. Further comparison... Figure 5 In Figures 5(a) and 5(b), the second reflection interface shifted significantly forward due to the TBM excavation disturbance. However, due to the limitation of the observation aperture, the imaging results could not accurately depict the morphological changes of the disturbance zone, requiring joint interpretation with the aforementioned wave velocity inversion results. Therefore, in summary, the reverse time migration imaging method can effectively realize dynamic imaging of spatial structural changes during the development of the TBM disturbance zone.

[0051] To further test the applicability of the proposed method to the bending interface model, this paper further uses the fault fracture zone model shown in Figure 6(a) for testing. Figure 7 As shown in (a), after the TBM excavation disturbance, the low wave velocity region in the fault core extends towards the TBM tunnel. The wave velocity change caused by the TBM excavation is as follows: Figure 7 As shown in (b), the anomalies are mainly concentrated within the fracture zone and continuously distributed along the tortuous zonal structure. Using the baseline wave velocity inversion results shown in Figure 6(b) as the initial model, the monitoring period inversion results are obtained as follows: Figure 7 As shown in (c), the location and morphology of the overall fault fracture zone interface are basically consistent with the real model. Meanwhile, the wave velocity inversion results are quite similar to the actual wave velocity distribution. By subtracting from the initial model, further wave velocity perturbation change inversion results were obtained, as shown in (c). Figure 7 As shown in (d), the wave velocity disturbance region and morphology are basically consistent with the actual disturbance, with significant differences only in the upper region of the disturbance. This indicates that the method presented in this paper can accurately invert the spatial range and geometric morphology of the rock mass temporal migration caused by TBM excavation. The interface inversion results further support the above conclusions: Figure 7 In (e), the inverted interface, shown by the blue dashed line, generally follows the same orientation and relative position as the real interface. However, due to the limitation of a finite number of first-order B-splines, the shape of the inverted result differs somewhat from the real interface. Further analysis... Figure 7 The results of the perturbation model interface inversion shown in (f) demonstrate that the proposed method can accurately invert the interface of the perturbation model, with its relative position and shape largely consistent with the real model. Finally, this paper quantitatively compares the accuracy of the wave velocity inversion results. Figure 7 (g) shows that the inverted wave velocity curve converges significantly to the true wave velocity compared to the initial model, and the low wave velocity region can be located well. Figure 7 The result (h) indicates that the wave velocity disturbance change obtained by the inversion is consistent with the actual result in terms of the anomaly center location and the sign of the change. In summary, the time-shifted wave velocity inversion method proposed in this paper can effectively realize the dynamic characterization of wave velocity disturbance changes caused by TBM excavation disturbance and the wave velocity distribution after disturbance.

[0052] Based on the time-shift inversion results, this paper obtains dynamic imaging results such as Figure 8 As shown. Compared with the previous imaging result. Figure 5 In comparison, the imaging results of this model better represent the trend of the real boundary, especially above the model, where the interface position and trend are also well imaged. This is due to the larger boundary tilt angle of the model, resulting in sufficient effective reflection signals. However, because the reflection information below the model cannot be well received by the detector array, its imaging results deviate significantly from the reality. Consequently, the imaging results above and below the interface are not on the same curve, thus affecting the effective identification of the reflection interface position. Further comparison... Figure 8 The imaging results before and after the disturbance shown in Figure (a) and Figure (b) of Figure 8, the reverse time migration imaging method can better invert the development of the second reflection interface. The imaging results are basically consistent with the actual interface position and tendency. However, due to the small observation aperture, it is difficult to reflect the morphology of the upper interface. Nevertheless, the dynamic imaging results can provide sufficient support for the identification of the development state of the TBM disturbance zone.

[0053] In summary, the wave velocity dynamic imaging method based on time-shifted reflection tomography proposed in this paper can not only effectively reconstruct the wave velocity distribution after TBM excavation disturbance, but also, based on the inverse time migration imaging method, further realize the dynamic imaging of the spatial structure changes in the disturbance zone, accurately characterizing the range and morphology of the medium disturbance caused by TBM excavation. The effectiveness of this method in the dynamic detection of rock mass instability in the TBM excavation disturbance zone of fault fracture zone is verified by inversion results and comparative analysis.

[0054] As a further specific explanation, in this embodiment, the specific process of using the tunnel array acoustic path function equation wave velocity-interface joint travel time inversion method includes: The basic inversion parameters are obtained, including the deployment parameters of the TBM tunnel array acoustic observation system, the observation data, the arrival time of direct waves and reflected waves, and the number of in-phase axes of reflected acoustic waves. Based on the stratification information of the surrounding rock medium, an initial model is established, which includes the interface location and wave velocity distribution. Based on the inversion basic parameters and the initial model, the incident wave travel time field and the reflected wave travel time field are constructed using the equation of process function, and the simulated travel time of the direct wave and the simulated travel time of the reflected wave are obtained by solving. A hierarchical joint loss function is constructed that includes both absolute travel time loss term and differential travel time loss term. The average wave velocity update gradient and the B-spline reflection interface control coefficient update gradient are solved respectively based on the hierarchical joint loss function. A layer-by-layer alternating iterative strategy was adopted to alternately update the medium average wave velocity and the B-spline reflection interface control coefficient, thereby completing the joint inversion of the wave velocity distribution of the surrounding rock and the geometric morphology of the reflection interface in the TBM excavation disturbance zone.

[0055] In this embodiment, the incident wave travel time field and the reflected wave travel time field are constructed based on the inversion basic parameters and the equation, including: Using the location of the seismic source as the global initial zero-to-time boundary condition, the equation of the process function is solved, and the global solution yields the simulated travel-time field of the incident wave. Using the simulated travel time field of the incident wave at the geological reflection interface as the initial boundary condition for the global reflection wave, the equation of process function is solved again, and the simulated travel time field of the reflection wave is obtained by solving the global equation. Based on the simulated travel time fields of the incident wave and the simulated travel time fields of the reflected wave, the theoretical arrival times of the corresponding direct wave and layered reflected wave of each detector are extracted, providing theoretical numerical support for the calculation of the loss function.

[0056] In this embodiment, the gradient update method for solving the medium average wave velocity based on the hierarchical joint loss function includes: The specific adjoint field equation of the reflected wave is constructed based on the adjoint state method, and the adjoint field of the reflected wave is obtained by substituting it into the simulated travel time field of the reflected wave. Combining the mathematical relationship that the slowness of sound waves and the velocity of medium waves are reciprocals, the original gradient of global wave velocity is generated by the global superposition of the reflection-associated field. The spatial media partitioning is completed by retrieving the prior information of the surrounding rock medium layering. The original global wave velocity gradients corresponding to all grids within the same medium partition are accumulated and summed to generate the medium average wave velocity update gradient with uniformity for each partition.

[0057] In this embodiment, the gradient update method for solving the B-spline reflection interface control coefficients based on the hierarchical joint loss function includes: Gradient transformation calculation is performed based on the wave velocity difference between two adjacent media layers, and the original gradient of the full-space interface geometry is generated from the original global wave velocity gradient. Multi-order, multi-segment B-spline basis functions are used to complete the parametric characterization of the geometry of all reflective interfaces, and a one-to-one mapping mathematical relationship between interface spatial coordinates and B-spline control coefficients is established. The original geometric gradient of the full-space interface is projected along the dimension of the B-spline basis function, and then a second-order curvature regularization constraint term is superimposed to suppress interface overfitting, finally generating the B-spline reflection interface control coefficients to update the gradient.

[0058] In this embodiment, a hierarchical joint loss function is constructed that simultaneously includes an absolute travel time loss term and a differential travel time loss term, including: For the observation data corresponding to each layer of reflected wave, a single-layer absolute travel time loss is generated. The absolute travel time loss is the weighted average of the simulated travel time of the reflection and the square of the original observed time difference of the reflected sound wave under the same source and the same detector combination. Differential travel loss is generated in layers. The differential travel loss of the first layer is the weighted average of the squared difference between the difference between the simulated travel time of the reflected wave and the simulated arrival time of the direct wave and the original observed time difference corresponding to the simulated arrival time difference. The differential travel loss of the remaining layers is the weighted average of the squared difference between the simulated travel time of the reflected wave in the current layer and the simulated arrival time of the reflected wave in the previous layer and the original observed time difference corresponding to the simulated arrival time difference. Set the differential loss growth rate to increase linearly with the iteration rounds. In each iteration, add the corresponding weight according to the growth rate to the differential travel loss. The weighted fusion of the absolute travel loss of the same layer and the weighted differential travel loss generates the total loss of a single layer. The total loss of each single layer corresponding to all reflection layers is combined and summarized to form a complete layered joint loss function.

[0059] In this embodiment, a layer-by-layer alternating iterative strategy is used to alternately update the medium's average wave velocity and the B-spline reflection interface control coefficients, including: According to the spatial layering order of the reflection interface from near to far and from shallow to deep, complete iterative calculations are carried out on each reflection interface and wave velocity separately. Within a single layer, continuous iteration rounds are divided, and differential parameter update operations are performed between odd-numbered and even-numbered iteration rounds. In a single-layer even-numbered iteration, the gradient of the B-spline reflection interface control coefficient is retrieved to correct and update the currently stored B-spline interface control parameters. In a single-layer odd-numbered iteration cycle, the average wave velocity of the medium is retrieved to update the gradient, and the currently stored average wave velocity of the layered medium is corrected and updated. After a single layer completes all preset iteration rounds of calculations, the average wave velocity of the medium and the B-spline reflection interface control coefficients obtained from the final update of this iteration are all saved as the initial input model for the next layer's reflection interface iteration calculation.

[0060] In some embodiments, the preprocessing of observation data for each period includes: estimating the tilt information and structural coherence of the acoustic signal using the multi-scale structural tensor method, setting a coherence threshold using the correlation between structural coherence and the acoustic signal, and dividing the acoustic signal into a signal region and a noise region. A hybrid total variation denoising method based on data tilt angle is adopted to perform adaptive constraint denoising on the signal region and the noise region, thereby suppressing background noise and protecting the effective signal. The hybrid total variation denoising method based on data tilt angle guidance performs adaptive constraint denoising on the signal region and the noise region, including: For the signal region, where the tilt information is stable, a total variational denoising method guided by the data tilt angle is used for noise suppression. By reducing the denoising intensity in the signal direction, the signal can be effectively protected. For the noise region, where there is no obvious tilt information, a total variational denoising method is directly used for noise suppression, thereby effectively balancing signal leakage and noise suppression.

[0061] As a further step, the acoustic signal detected by the TBM disturbance zone array is obtained, wherein the acoustic signal contains valid signal and noise data; the acoustic signal denoising process is regarded as an inverse problem, that is, recovering the valid signal from the noisy acoustic signal.

[0062] Furthermore, the tilt angle information and structural coherence of the acoustic signal are estimated using the multi-scale structure tensor method. This multi-scale structure tensor method quantifies the local structural features of a pixel by analyzing the gradient in its neighborhood, including: First, the acoustic data is filtered using a derivative Gaussian window function. Then, a ring filter is used to filter the structure tensor. This filter can align with the feature edges and accurately normalize the edge energy to the feature center, thus obtaining the filtered set of structure tensors. The filtered structure tensors at different scales are normalized, and the trace of the structure tensor at each scale is solved at each acoustic data location. The optimal filter scale at each pixel is obtained through one-dimensional maximum search. By extracting the structure tensor corresponding to the optimal filter scale pixel by pixel, the optimal structure tensor is obtained. The eigenvalues ​​of the optimal structure tensor are then solved to obtain the signal structure tilt angle information. The coherence of the local structure is calculated using the eigenvalues ​​of the structure tensor, yielding the results of the local structure coherence.

[0063] Furthermore, by utilizing the correlation between structural coherence and acoustic signal to set a coherence threshold, the acoustic signal is divided into a signal region and a noise region, including: Coherence is an evaluation index distributed between [0, 1]. When a correlation threshold is set, when the local coherence C=0, the local structure exhibits obvious isotropy, that is, there is basically no change in the amplitude of the signal, so the region is judged as noise; when the local coherence C=1, the local structure exhibits obvious anisotropy, so it is judged as a valid signal. By employing a coherence threshold, the coherence results are converted into a binary mask matrix, thereby distinguishing between the in-phase region and the noise region of the acoustic data.

[0064] Furthermore, for the signal region where the tilt information is stable, a data tilt-guided total variation denoising method is used for noise suppression, including: The obtained acoustic signal structure tilt angle is introduced to reduce the noise reduction intensity along the phase axis direction, while increasing the noise reduction intensity in the phase axis normal direction; The scaling factor is set to adjust the denoising intensity in the tangential and normal directions of the in-phase axis. The smaller the scaling factor, the weaker the denoising intensity in the in-phase axis direction, which can protect the in-phase axis. When the scaling factor = 1, the denoising intensity in both directions is the same.

[0065] As a further step, a hybrid total variation denoising method based on data tilt angle is used to adaptively constrain denoising the signal and noise regions, as shown in the following formula: ; Where J is the loss function, s , These are the valid signal and the valid signal obtained after denoising, respectively. d The original data, λ It is a scaling factor; TV and DTV denotes the traditional total variation denoising method and the total variation denoising method based on data tilt angle, respectively; C is structural coherence, and τ is the coherence threshold.

[0066] In some embodiments, the preprocessing of observation data for each period also includes: obtaining a clean and effective signal based on denoising preprocessing; selecting the multiple wave development region in the acoustic data through a pre-windowing strategy; projecting the multiple reflection waves and the primary reflection waves to the feature space using the generalized singular value decomposition method; and then reconstructing the primary reflection wave component of the acoustic data through low-rank component reconstruction to obtain the reconstructed data. By employing the local correlation method, the local similarity between the reconstructed data and the effective signal is calculated. A similarity threshold is set, and a binary mask for the primary reflection wave data is constructed and applied to the reconstructed data. This effectively suppresses reconstruction noise and residual multiple reflection waves, thereby achieving effective separation of the primary reflection wave and multiple reflection waves.

[0067] Specifically, a pre-windowing strategy is used to select the regions where multiple reflections develop in the acoustic data. Then, the generalized singular value decomposition method is used to project the multiple and primary reflections onto the feature space. Finally, low-rank component reconstruction is used to reconstruct the primary reflection component of the acoustic data, resulting in reconstructed data, including: First, the multiple wave data in the effective signal are selected by box selection, and the effective signal is divided into a region where the primary and multiple reflection waves overlap and a region where multiple reflection waves overlap, and data matrices for the two regions are constructed. The two data matrices are decomposed using the generalized singular value decomposition method to obtain the eigenvector matrix and the diagonal matrix, respectively. By comparing the magnitudes of the singular values ​​in the diagonal matrix, the dominant component of the primary reflected wave is obtained based on the comparison results, and the reconstructed data is obtained.

[0068] The characteristics of single-reflection wave data and multiple-reflection wave data are as follows: the reconstructed data of single-reflection wave data has a stable shape and a high similarity with the original effective signal; the reconstruction noise is newly introduced noise during the reconstruction process and has no correlation; the residual multiple-reflection wave signal has a low similarity with the original effective signal due to its incomplete energy or structure.

[0069] The local similarity index is set within the range of [0,1], with a higher similarity between the two matrices at that point. Utilizing the high correlation of primary reflections and the low coherence of reconstruction noise and residual multiple reflections, a mask matrix is ​​constructed by setting a similarity threshold. ; in, α Let Γ be the attenuation coefficient, Γ be the local similarity index, and the mask matrix represent the spatial location information of the first reflected wave, i.e., 1 indicates the presence of a first reflected wave, and 0 indicates the absence of a first reflected wave. By using a similarity threshold, reconstruction noise and the spatial location information of residual multiple reflected waves are removed. Therefore, the mask matrix only contains the spatial location information of the first reflected wave. This information is then used to extract the first reflected wave and remove the multiple reflected waves.

[0070] The local correlation method is used to calculate the local similarity between the reconstructed data and the effective signal, including: Construct the global correlation coefficient of the signal; By using the vector matrix of global correlation coefficients, an approximate solution is obtained through an integer regularization method. Finally, a local similarity index is constructed using the approximate solution values.

[0071] In some embodiments, the process of picking up the arrival time of observation data in each period includes: using the kurtosis feature of the transmitted signal as the arrival time indicator, and constructing a path cost function based on the kurtosis feature of the observation data; Apply hard constraints to the arrival candidate points between adjacent channels; Under the hard constraints, based on the path cost function, the problem of picking up the arrival time of the reflected wave is regarded as a decision problem of picking up the wave one by one along the phase axis of the reflected wave. The result of picking up the arrival time of the reflected wave is regarded as an arrival curve. The decision problem is solved recursively by dynamic programming to obtain the preliminary automatic picking result. The initial automatic picking results are fine-tuned based on the reliability to obtain the final automatic picking results.

[0072] Furthermore, the process of using the kurtosis characteristics of the transmitted signal as a time-of-arrival indicator includes: the kurtosis characteristics of the signal are: ; in, μ Indicates signal The mean, m 2 and m 4 represents the second-order central moment and the fourth-order central moment. Indicates signal Total quantity Indicates time.

[0073] Furthermore, by performing a first-order difference on the signal kurtosis and removing the negative variation, the non-Gaussian abrupt change from noise to the event is highlighted, as follows: ; Here, max(·) represents taking the maximum of the two.

[0074] Furthermore, the process of constructing the path cost function based on the kurtosis characteristics of the observed data includes: Path cost function ; Where, d K + Indicates the kurtosis characteristics of the observed data. λ As a scaling factor, τ i This indicates the pickup result when the reflected wave corresponding to channel i arrives. Let x represent time, x represent the signal, and the first term be the confidence cost function, calculated by observing the kurtosis characteristic d of the data. K + The first term characterizes the reliability of the candidate point as the arrival time of the reflected wave. The higher the energy, the higher the probability of it being the correct arrival time of the reflected wave, and the smaller the loss. The second term is a smoothing function, which describes whether the travel time changes between two adjacent candidate points are smooth.

[0075] Furthermore, the process of applying hard constraints to the candidate points for arrival times between adjacent channels includes: hard constraints ; Where, Δ t max For the maximum time interval, τ i This represents the arrival time picking result of the reflected wave corresponding to the i-th channel. This hard constraint is used to achieve the selection of the arrival time of the reflected wave within a finite neighborhood.

[0076] Furthermore, the problem of picking up the arrival time of reflected waves can be viewed as a decision-making problem involving picking up each channel along the phase axis of the reflected waves, and the result of picking up the arrival time of reflected waves can be viewed as a process of an arrival curve: The problem of picking up the arrival time of reflected waves can be viewed as a decision problem involving picking up each wave one by one along the phase axis of the reflected wave. The result of picking up the arrival time of reflected waves can be considered as an arrival curve, which satisfies the following: ; in, Indicates the curve / path of the reflected wave arrival time. τ i This indicates the pickup result when each corresponding reflected wave arrives. N rec Given the number of detectors, based on the subproblem optimization idea of ​​dynamic programming, the picking result when each reflected wave arrives should satisfy the optimal path cost function.

[0077] Furthermore, the optimal path cost function is: ; Cumulative path cost function Minimum.

[0078] Furthermore, the process of solving the decision problem using dynamic programming recursion includes: using dynamic programming recursion, we obtain: ; in, Let Δ be the cost function of the optimal path. t max For the maximum time interval, τ i d represents the pickup result when the reflected wave corresponding to channel i arrives. K + Indicates the kurtosis characteristics of the observed data. λ This is a scaling factor; it records the optimal index for each path and obtains the entire path by backtracking.

[0079] Furthermore, the process of fine-tuning the preliminary automatic picking results based on credibility to obtain the final automatic picking results includes: introducing expert experience to correct or rectify the automatic picking results; manually picking multiple arrival points; automatically fine-tuning based on the credibility cost function to ensure that the manually selected points fall at the positions with the strongest kurtosis; and using the manually picked arrival points as anchor points and adding them to the path cost function to obtain: ; Where, d K + Indicates the kurtosis characteristics of the observed data. λ As a scaling factor, τ i This indicates the pickup result when the reflected wave corresponding to channel i arrives. This indicates manually selected anchor points. Nt Indicates the sampling point. γ As a scaling factor, when the complexity of the in-phase axis of the reflection data exceeds a predetermined level, γ >0, introduce soft constraints to make the path close to the prior, and realize accurate picking of the arrival time of complex reflected wave phase axis.

[0080] Example 2 A TBM-based close-array acoustic dynamic detection system for rock mass instability in a fault fracture zone includes: The preprocessing module is configured to acquire observation data from multiple periods and preprocess the observation data from each period. The time-out pickup module is configured to pick up the preprocessed data at the time-out point to obtain multiple time-out periods; The inversion module is configured to perform baseline inversion and monitoring inversion based on the multi-period travel times. The baseline inversion obtains the baseline wave velocity inversion result by iteratively fitting the arrival time of the observation data, and uses it as the initial model input into the monitoring inversion. The monitoring inversion reconstructs the wave velocity distribution of the TBM disturbance area at different times by fitting the double difference of the arrival time data of the two periods before and after, and obtains the monitoring wave velocity inversion result. The reverse time migration imaging module is configured to use the reference wave velocity evolution results and the monitoring wave velocity inversion results as background models, and to perform dynamic imaging of the TBM excavation disturbance zone development process using the reverse time migration imaging method.

[0081] Example 3 A dynamic imaging and control device for the development process of TBM excavation disturbance zone includes: Detector array 100, industrial control computer 200 and TBM actuator group 300.

[0082] The detector array 100 is used to receive seismic wave signals induced by surrounding rock disturbance in real time during TBM tunneling. It includes: multiple detectors 110: arranged on the surface of the TBM cutterhead spokes and the outer periphery of the shield in a non-equidistant manner, used to convert mechanical vibration into electrical signals; a multi-channel data acquisition unit 120: electrically connected to the detectors 110, amplifying, filtering and converting the electrical signals into analog and digital signals; and a high-speed data transmission interface 130: transmitting the raw waveform data to the industrial control computer 200 in real time via fiber optic Ethernet, with a transmission delay of ≤5ms.

[0083] The detector array 100 is directly connected to the industrial control computer 200. The raw waveform data is directly input into the real-time memory area of ​​the industrial control computer 200 without being processed by external relay equipment.

[0084] The industrial control computer 200 is the core processing and decision-making unit of this device, and it has both dynamic imaging and control functions. (1) Dynamic imaging function The industrial control computer 200 receives the raw waveform data input from the detector array 100, executes the dynamic imaging method described in Embodiment 1 of this application, reconstructs a three-dimensional time-series image of the disturbance area in front of the cutter head, and extracts the following feature parameters: Maximum value of disturbance index ( ); Maximum displacement rate gradient ( ); Collapse cavity volume ( ); Real-time torque (T) and torque change rate (dT / dt) of the cutter head.

[0085] (2) Instability Critical Judgment Function The industrial control computer 200 has a built-in discrimination engine that determines the status based on the aforementioned characteristic parameters: Pre-instability state: and ; Instability state: ,or ; Machine jamming risk status: The sudden increase rate of dT / dt relative to the rated torque within 5 seconds is >150%.

[0086] (3) Control command generation function After determining the engine output status code, the control instruction generation unit automatically generates the corresponding instruction set: Normal state: Maintain current tunneling parameters; Pre-instability condition: The cutter head speed drops to 40%–60% of the rated value, the total thrust drops to 70% of the rated value, and the maximum allowable load of a single hob drops to 80% of the rated value; Instability state: The cutterhead stops rotating and remains idle, the total thrust drops to the sustaining thrust, the advanced grouting valve is opened and the grouting pressure is dynamically adjusted according to the volume of the collapsed body, and the support shoe pressure is adjusted so that the shield eccentric load coefficient is ≤1.15; Machine jamming risk condition: The cutter head reverses, the propulsion cylinder performs pulsating propulsion, and lubricating slurry is injected.

[0087] In this embodiment, control commands are directly sent to the TBM actuator group 300 via the EtherCAT bus.

[0088] The TBM actuator group 300 receives control commands from the industrial control computer 200 and executes corresponding physical actions, including: Cutter head frequency converter driver 310: Adjusts the rotation state of the cutter head according to the speed and direction commands; Propulsion cylinder proportional valve group 320: Adjusts the oil supply pressure of each propulsion cylinder group according to the thrust command; Hydraulic servo valve assembly 330 for support shoes: Independently adjusts the reaction force of each support shoe according to the pressure distribution command; Advanced grouting solenoid valve 340: Opens and closes the grouting passage according to the grouting trigger command, and adjusts the grouting pressure in conjunction with the grouting pump.

[0089] In this embodiment, the control process further includes: During the TBM tunneling process, the detector array 100 continuously inputs the raw waveform data into the industrial control computer 200; the industrial control computer 200 performs dynamic imaging in real time, acquires the characteristic parameters of the disturbance area, and generates a status code through the discrimination engine, and automatically sends control commands to the TBM actuator group 300, forming a millisecond-level response closed loop. When the engine determines that the surrounding rock has entered an unstable state, the industrial control computer 200 synchronously outputs a stop command to the cutterhead frequency converter driver 310, a pressure reduction command to the propulsion cylinder proportional valve group 320, an opening command to the advanced grouting solenoid valve 340, and an off-center load balancing command to the support shoe hydraulic servo valve group 330. It should be noted that each actuator completes its action within 50 ms, and the TBM switches from the tunneling state to the reinforcement-stepping state.

[0090] When the results of consecutive multi-frame imaging are displayed When the volume of the collapsed body is significantly reduced below the recovery threshold, the engine automatically de-stresses and the industrial control computer 200 switches the control commands to either low-disturbance passage mode or normal tunneling mode.

[0091] In this embodiment, the detector array 100 is directly connected to the industrial control computer 200 via fiber optic Ethernet, and the industrial control computer 200 is directly connected to the TBM actuator group 300 via an EtherCAT bus. The industrial control computer 200 serves as a unified hardware carrier for waveform data processing, dynamic imaging reconstruction, instability detection, and control command generation, forming an inseparable electromechanical-hydraulic integrated closed loop with the TBM's mechanical structure and hydraulic system. Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0092] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations and / or block diagrams. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0093] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0094] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0095] The above description is merely a preferred embodiment of the present invention and is not intended to limit the invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made by those skilled in the art without creative effort within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for dynamic detection of rock mass instability in a fault fracture zone using dense array acoustic waves via TBM, characterized in that, Includes the following steps: Acquire observation data from multiple periods, preprocess the observation data from each period, and then pick up the time of arrival of the preprocessed data to obtain the travel times for multiple periods; Based on the multi-period travel times, benchmark inversion and monitoring inversion are performed. The benchmark inversion obtains the benchmark wave velocity inversion result by iteratively fitting the arrival time of the observation data, and uses it as the initial model input into the monitoring inversion. The monitoring inversion reconstructs the wave velocity distribution of the TBM disturbance area at different times by fitting the double difference of the arrival time data of the two periods before and after, and obtains the monitoring wave velocity inversion result. Using the monitored wave velocity inversion results as a background model, the reverse time migration imaging method was used to perform dynamic imaging of the TBM excavation disturbance zone development process.

2. The method for dynamic detection of rock mass instability in a fault fracture zone using dense array acoustic waves as described in claim 1, characterized in that, The preprocessing of observation data for each period includes: estimating the tilt information and structural coherence of acoustic signals in the observation data using the multi-scale structure tensor method, setting a coherence threshold using the correlation between local coherence and acoustic signals, and dividing the acoustic signals into signal regions and noise regions. A hybrid total variational denoising random noise suppression method based on structural tilt angle is adopted to adaptively denoise the signal region and the noise region, thereby suppressing the background noise. The hybrid total variational denoising random noise suppression method based on structural tilt angle guidance adaptively denoises the signal region and the noise region, including: For the signal region, where the tilt information is stable, a total variation denoising method guided by the data tilt angle is used for noise suppression. For the noise region, where there is no obvious tilt information, a total variation denoising method is directly used for noise suppression, thereby effectively balancing signal leakage and noise suppression.

3. The method for dynamic detection of rock mass instability in a dense array acoustic wave traversing a fault fracture zone using a TBM as described in claim 2, characterized in that... The tilt angle information and structural coherence of acoustic signals are estimated using the multi-scale structure tensor method. This method quantifies the local structural features of a pixel by analyzing the gradient in its neighborhood, including: First, the acoustic data is filtered using a derivative Gaussian window function. Then, a ring filter is used to filter the structure tensor. This filter can align with the feature edges and accurately normalize the edge energy to the feature center, thus obtaining the filtered set of structure tensors. The filtered structure tensors at different scales are normalized, and the trace of the structure tensor at each scale is solved at each acoustic data location. The optimal filter scale at each pixel is obtained through one-dimensional maximum search. By extracting the structure tensor corresponding to the optimal filter scale pixel by pixel, the optimal structure tensor is obtained. The eigenvalues ​​of the optimal structure tensor are then solved to obtain the signal structure tilt angle information. The coherence of the local structure is calculated using the eigenvalues ​​of the structure tensor, yielding the results of the local structure coherence.

4. The method for dynamic detection of rock mass instability in a fault fracture zone using dense array acoustic waves as described in claim 2, characterized in that, The tilt angle information and structural coherence of acoustic signals are estimated using the multi-scale structure tensor method. The multi-scale structure tensor method quantifies the local structural features of a pixel by analyzing the gradient in its neighborhood. First, the acoustic data is filtered using a derivative Gaussian window function. Then, a ring filter is used to filter the structure tensor. This filter can align with the feature edges and accurately normalize the edge energy to the feature center, thus obtaining the filtered set of structure tensors. The filtered structure tensors at different scales are normalized, and the trace of the structure tensor at each scale is solved at each acoustic data location. The optimal filter scale at each pixel is obtained through one-dimensional maximum search. By extracting the structure tensor corresponding to the optimal filter scale pixel by pixel, the optimal structure tensor is obtained, and the eigenvalues ​​of the optimal structure tensor are solved to obtain the signal structure tilt angle information. The coherence of the local structure is calculated using the eigenvalues ​​of the structure tensor, yielding the results of the local structure coherence.

5. The method for dynamic detection of rock mass instability in a fault fracture zone using dense array acoustic waves as described in claim 2, characterized in that, The preprocessing of observation data for each period also includes: obtaining a clean effective signal based on denoising preprocessing; selecting multiple wave data in the effective signal by pre-windowing strategy; extracting the primary reflection wave component by using generalized singular value decomposition and low-rank approximation; reconstructing the primary reflection wave component in the effective signal to obtain reconstructed data. The local correlation method is used to calculate the local similarity between the reconstructed data and the effective signal. A similarity threshold is set, a binary mask of the primary reflection wave data is constructed and applied to the original clean signal. At the same time, the reconstructed noise and the spatial location information of the residual multiple reflection waves are removed, so as to achieve effective separation of the primary reflection wave and the multiple waves.

6. The method for dynamic detection of rock mass instability in a dense array acoustic wave traversing a fault fracture zone using a TBM as described in claim 5, characterized in that... Multiple wave data in the effective signal are selected by pre-windowing strategy. The primary reflection wave component is extracted using generalized singular value decomposition and low-rank approximation. The primary reflection wave component in the effective signal is then reconstructed to obtain the reconstructed data, including: First, the multiple wave data in the effective signal are selected by box selection, and the effective signal is divided into a region where the primary and multiple reflection waves overlap and a region where multiple reflection waves overlap, and data matrices for the two regions are constructed. The two data matrices are decomposed using the generalized singular value decomposition method to obtain the eigenvector matrix and the diagonal matrix, respectively. By comparing the magnitudes of the singular values ​​in the diagonal matrix, the dominant component of the primary reflected wave is obtained based on the comparison results, and the reconstructed data is obtained.

7. The method for dynamic detection of rock mass instability in a fault fracture zone using dense array acoustic waves as described in claim 1, characterized in that, The process of acquiring observation data from multiple periods and picking up arrival times for each period of observation data includes: using the kurtosis characteristics of the transmitted signal of the observation data as the arrival time indicator, and constructing a path cost function based on the kurtosis characteristics of the observation data; Apply hard constraints to the arrival candidate points between adjacent channels; Under the hard constraints, based on the path cost function, the problem of picking up the arrival time of the reflected wave is regarded as a decision problem of picking up the wave one by one along the phase axis of the reflected wave. The result of picking up the arrival time of the reflected wave is regarded as an arrival curve. The decision problem is solved recursively by dynamic programming to obtain the preliminary automatic picking result. The initial automatic picking results are fine-tuned based on the reliability to obtain the final automatic picking results.

8. The method for dynamic detection of rock mass instability in a fault fracture zone using dense array acoustic waves as described in claim 7, characterized in that, The process of constructing the path cost function includes: Path cost function ; Where, d K + Indicates the kurtosis characteristics of the observed data. λ As a scaling factor, τ i This indicates the pickup result when the reflected wave corresponding to channel i arrives. Let x represent time, x represent the signal, and the first term be the confidence cost function, calculated by observing the kurtosis characteristic d of the data. K + The first term characterizes the reliability of the candidate point as the arrival time of the reflected wave. The higher the energy, the higher the probability of it being the correct arrival time of the reflected wave, and the smaller the loss. The second term is a smoothing function, which describes whether the travel time changes between two adjacent candidate points are smooth.

9. The method for dynamic detection of rock mass instability in a fault fracture zone using dense array acoustic waves as described in claim 7, characterized in that, The process of solving the decision problem using dynamic programming recursion includes: using dynamic programming recursion, we obtain: ； in, Let Δ be the cost function of the optimal path. t max For the maximum time interval, τ i d represents the pickup result when the reflected wave corresponding to channel i arrives. K + Indicates the kurtosis characteristics of the observed data. λ This is a scaling factor; it records the optimal index for each path and obtains the entire path by backtracking.

10. The method for dynamic detection of rock mass instability in a fault fracture zone using dense array acoustic waves as described in claim 7, characterized in that, The process of fine-tuning the initial automatic picking results based on credibility to obtain the final automatic picking results includes: introducing expert experience to correct or rectify the automatic picking results; manually picking multiple arrival points; automatically fine-tuning based on the credibility cost function to ensure that the manually selected points fall at the positions with the strongest kurtosis; and using the manually picked arrival points as anchor points and adding them to the path cost function to obtain: ； Where, d K + Indicates the kurtosis characteristics of the observed data. λ As a scaling factor, τ i This indicates the pickup result when the reflected wave corresponding to channel i arrives. This indicates manually selected anchor points. N t Indicates the sampling point. γ As a scaling factor, when the complexity of the in-phase axis of the reflection data exceeds a predetermined level, γ >0, introduce soft constraints to make the path close to the prior, and realize accurate picking of the arrival time of complex reflected wave phase axis.

11. The method for dynamic detection of rock mass instability in a fault fracture zone using dense array acoustic waves as described in claim 1, characterized in that, Based on the aforementioned multi-period travel times, the process of benchmark inversion and monitoring inversion includes: referring to the data obtained from a certain detection as the benchmark data arrival time, referring to the subsequent detection data as the monitoring data arrival time, using the benchmark data arrival time as the benchmark state, and using the wave velocity and interface results obtained from its inversion as the background model; extracting the travel time change from the observation data of the subsequent monitoring period to minimize the impact of static structural errors and systematic time delays on the inversion; and solving for the time shift change of the disturbance region based on the difference between the inversion results of the two periods of detection.

12. The method for dynamic detection of rock mass instability in a dense array acoustic wave traversing a fault fracture zone using a TBM as described in claim 1, characterized in that, Based on the multi-period travel times, the process of benchmark inversion includes: obtaining the basic inversion parameters, which include the deployment parameters of the TBM tunnel array acoustic observation system, the prior information of the surrounding rock medium layering, the original observation time of the incident acoustic wave and the original observation time of the reflected acoustic wave. Based on the inversion basic parameters, and combined with the equation, forward modeling models of incident wave travel time field and reflected wave travel time field are constructed respectively, and the simulated travel time of incident wave and simulated travel time of reflected wave are obtained by solving them. A hierarchical joint loss function is constructed that includes both absolute travel time loss term and differential travel time loss term. The average wave velocity update gradient and the B-spline reflection interface control coefficient update gradient are solved according to the hierarchical joint loss function. The average wave velocity of the medium and the control coefficient of the B-spline reflection interface are updated alternately by adopting a layer-by-layer alternating iterative strategy to complete the joint inversion of the wave velocity distribution of the surrounding rock and the geometric morphology of the reflection interface in the TBM excavation disturbance zone.

13. The method for dynamic detection of rock mass instability in a fault fracture zone using dense array acoustic waves as described in claim 1, characterized in that, The monitoring and inversion process, which reconstructs the wave velocity distribution of the TBM disturbance zone at different times by fitting the double difference of arrival data from two consecutive periods, includes: using the tunnel reflection wave double difference travel-time inversion loss function, specifically: ； in, J f and J r Let these represent the event arrival loss and the receiver double difference loss function, respectively. For the receiving point residual, For the residual at that time, For coefficients; Based on the aforementioned inversion fundamental parameters, and in conjunction with the equation, forward modeling models of the incident wave travel time field and the reflected wave travel time field are constructed respectively, including: The equation is solved by taking the location of the seismic source as the global initial zero-to-time boundary condition, and the global solution yields the incident simulation travel time. Using the incident simulated travel time corresponding to the geological reflection interface as the global initial boundary condition for the reflected wave, the equation is solved again to obtain the global simulated travel time of the reflection. Based on the incident simulation travel time and the reflection simulation travel time, the theoretical arrival times of the corresponding direct waves and layered reflected waves of each detector are extracted, providing theoretical numerical support for the calculation of the loss function.

14. The method for dynamic detection of rock mass instability in a fault fracture zone using dense array acoustic waves as described in claim 13, characterized in that, The step of solving the medium average wave velocity update gradient based on the hierarchical joint loss function includes: Based on the adjoint state method, the incident wave specific adjoint field equation and the reflected wave specific adjoint field equation are constructed respectively. Substituting them into the incident simulation travel time and the reflection simulation travel time, the incident adjoint field and the reflection adjoint field are obtained by solving. Combining the mathematical relationship that the slowness of sound waves and the velocity of medium waves are reciprocals, the original gradient of global wave velocity is generated by the global superposition of the incident accompanying field and the reflected accompanying field. The spatial medium partitioning is completed by retrieving the prior information of the surrounding rock medium layering. The original global wave velocity gradients corresponding to all grids within the same medium partition are accumulated and summed to generate the medium average wave velocity update gradient that is uniform across the partition.

15. The method for dynamic detection of rock mass instability in a fault fracture zone using dense array acoustic waves as described in claim 13, characterized in that, The step of updating the gradient by solving the B-spline reflection interface control coefficients according to the hierarchical joint loss function includes: Gradient transformation calculation is performed based on the wave velocity difference between two adjacent media layers, and the original gradient of the full-space interface geometry is generated from the original global wave velocity gradient. Multi-order, multi-segment B-spline basis functions are used to complete the parametric characterization of the geometry of all reflective interfaces, and a one-to-one mapping mathematical relationship between interface spatial coordinates and B-spline control coefficients is established. The original geometric gradient of the full-space interface is projected along the dimension of the B-spline basis function, and a second-order curvature regularization constraint term to suppress interface overfitting is superimposed to finally generate the updated gradient of the B-spline reflection interface control coefficient.

16. The method for dynamic detection of rock mass instability in a fault fracture zone using dense array acoustic waves as described in claim 13, characterized in that, The construction of the hierarchical joint loss function, which simultaneously includes absolute travel loss and differential travel loss, includes: For the observation data corresponding to each layer of reflected wave, a single-layer absolute travel time loss is generated. The absolute travel time loss is the weighted average of the simulated travel time of the reflection and the square of the original observed time difference of the reflected sound wave under the same source and the same detector combination. The differential travel loss is generated in layers. The differential travel loss of the first layer is the weighted average of the simulated travel time of direct incident and the square of the original direct arrival time difference. The differential travel loss of the remaining layers is the weighted average of the simulated travel time of reflection in this layer and the square of the original observed time difference of the reflected sound wave in the previous layer. Set a differential loss growth rate that increases linearly with the iteration rounds. In each iteration, add the corresponding weight to the differential time loss according to the growth rate. The absolute time loss and the weighted differential time loss of the same layer are weighted and fused to generate the total loss of a single layer. The total loss of each single layer corresponding to all reflection layers is combined and summarized to generate the complete joint loss function of the layers.

17. The method for dynamic detection of rock mass instability in a fault fracture zone using dense array acoustic waves as described in claim 13, characterized in that, The method of alternately updating the average wave velocity of the medium and the B-spline reflection interface control coefficient using a layer-by-layer iterative strategy includes: According to the spatial layering order of the reflective interface from near to far and from shallow to deep, complete iterative calculations are performed on each reflective interface layer separately. Within a single layer, continuous iteration rounds are divided, and differential parameter update operations are performed between odd-numbered and even-numbered iteration rounds. In a single-layer even-numbered iteration cycle, the gradient of the B-spline reflection interface control coefficient is retrieved to correct and update the currently stored B-spline interface control parameters. In a single-layer odd-numbered iteration round, the average wave velocity update gradient of the medium is retrieved, and the currently stored average wave velocity of the layered medium is corrected and updated. After a single layer completes all preset iteration rounds of calculation, the average wave velocity of the medium and the control coefficients of the B-spline reflection interface obtained from the final update of this iteration are all saved as the initial input model for the next layer of reflection interface iteration calculation.

18. A TBM-based system for detecting rock mass instability in a fault fracture zone using a dense array of acoustic waves, characterized in that... include: The preprocessing module is configured to acquire observation data from multiple periods and preprocess the observation data from each period. The time-out pickup module is configured to pick up the preprocessed data at the time-out point to obtain multiple time-out periods; The inversion module is configured to perform baseline inversion and monitoring inversion based on the multi-period travel times. The baseline inversion obtains the baseline wave velocity inversion result by iteratively fitting the arrival time of the observation data, and uses it as the initial model input into the monitoring inversion. The monitoring inversion reconstructs the wave velocity distribution of the TBM disturbance area at different times by fitting the double difference of the arrival time data of the two periods before and after, and obtains the monitoring wave velocity inversion result. The reverse time migration imaging module is configured to use the monitored wave velocity inversion result as a background model and perform dynamic imaging of the TBM excavation disturbance zone development process using the reverse time migration imaging method.

Images