Advanced geological exploration and ground stress inversion method for coal mine tunneling roadway

Abstract

The invention discloses a coal mine tunneling roadway advanced geological exploration and ground stress inversion method, particularly relates to the technical field of seismic signal processing in geophysical exploration, and is used for solving the problems of geological exploration signal distortion and insufficient stress inversion precision caused by strong noise interference of tunneling equipment in the prior art. Rock mass vibration signals are collected through an optical fiber sensor arranged on a roadway wall, equipment noise characteristics are identified through spectral analysis, a quantitative mapping relation between working condition parameters and noise characteristics is established, a noise reference time period is determined according to the quantitative mapping relation, a spatial coherence characteristic template is constructed, high-coherence noise components are inhibited through template comparison, and the noise reference time period is determined according to the spatial coherence characteristic template. Finally, the geological structure and stress field distribution in front of the working face are inverted based on the purified seismic wave signals, a stress concentration area and a disaster risk area are identified, the signal-to-noise ratio of the seismic signals in the strong noise environment and the reliability of the geological inversion result are effectively improved, and accurate technical support is provided for safe mining of a coal mine.

Description

Technical Field

[0001] This invention relates to the field of seismic signal processing technology in geophysical exploration, and more specifically, to a method for advanced geological exploration and geostress inversion in coal mine tunneling. Background Technology

[0002] In the process of underground roadway excavation in coal mines, advanced geological exploration and stress inversion are key technical links to ensure safe and efficient production. Existing technologies typically employ geophysical exploration methods, especially those based on seismic wave detection principles, using sensor arrays deployed within the roadway to collect rock mass vibration signals, and then inverting the geological structure and stress field distribution characteristics ahead of the working face. Fiber optic sensing technology, due to its advantages such as distributed measurement and resistance to electromagnetic interference, is applied to stress-strain and vibration signal monitoring in mines. Based on the collected data, stress field models are constructed using inversion algorithms to identify potential stress concentration areas, and further, the risk level of disasters such as rockbursts and roof falls at the working face is assessed based on mechanical parameters, forming a complete technical system from data acquisition to disaster early warning.

[0003] However, the aforementioned technical methods face significant limitations in practical applications. When tunnel excavation equipment operates, the mechanical vibrations it generates create broadband, high-intensity interference noise. This noise overlaps significantly with the effective seismic signal in both the frequency and time domains, and its intensity is significantly higher than the useful signal, resulting in a severe deterioration of the signal-to-noise ratio of the acquired raw data. Existing signal processing techniques struggle to effectively separate this strong interference within the same frequency band, leading to decreased accuracy in seismic wave acquisition during travel, blurred identification of reflected wave groups, and ultimately, distortion of the geological structural interface morphology obtained through inversion and increased deviations in stress field calculations. This deficiency directly reduces the reliability of advanced detection and the accuracy of stress concentration zone identification, thereby weakening the confidence of subsequent disaster risk assessment results for the working face and posing potential hazards to safe and efficient mine mining. Summary of the Invention

[0004] In order to overcome the above-mentioned defects of the prior art, the present invention provides a method for advanced geological exploration and geostress inversion of coal mine tunneling roadways to solve the problems mentioned in the background art.

[0005] To achieve the above objectives, the present invention provides the following technical solution:

[0006] A method for advanced geological exploration and inversion of geostress in coal mine tunneling roadways includes the following steps:

[0007] S1. Vibration signals generated by the rock mass in front of the tunnel face are collected by fiber optic sensors installed on the tunnel wall;

[0008] S2. Perform spectrum analysis on the collected vibration signals to identify specific frequency band interference noise generated by the operation of the tunneling equipment;

[0009] S3. Real-time monitoring of multi-source operating parameters of tunneling equipment, and establishment of a quantitative mapping relationship between multi-source operating parameters and characteristic modes of interference noise in a specific frequency band;

[0010] S4. Based on the quantitative mapping relationship, determine the period when the interference noise is dominant as the noise reference period, and calculate the coherence between the multi-channel vibration signals within the noise reference period to establish the spatial coherence feature template of the noise.

[0011] S5. The vibration signal is compared with the spatial coherence feature template, and the signal components with high spatial coherence are suppressed to obtain the purified effective seismic wave signal.

[0012] S6. Based on the purified effective seismic wave signal, invert the geological structure and stress field distribution in front of the working face, and identify stress concentration areas and disaster risk areas based on the stress field distribution.

[0013] Furthermore, vibration signals generated by the rock mass in front of the tunnel face are collected using fiber optic sensors deployed on the tunnel walls, including:

[0014] Fiber optic sensors are deployed at equal intervals in the space between the tunnel roof and the rock walls on both sides.

[0015] Vibration signals, including sound waves and seismic waves, are continuously acquired using fiber optic sensors in a distributed measurement manner.

[0016] Furthermore, spectral analysis was performed on the collected vibration signals to identify specific frequency band interference noise generated by the tunneling equipment, including:

[0017] The vibration signal is subjected to a fast Fourier transform to obtain a spectrum.

[0018] Identify continuous spectral peaks in the spectrum that match the operating frequency characteristics of the tunneling equipment;

[0019] The frequency band boundary of interference noise in a specific frequency band is determined based on the distribution range of continuous spectral peaks.

[0020] Furthermore, multi-source operating parameters of the tunneling equipment are monitored in real time, and a quantitative mapping relationship is established between multi-source operating parameters and characteristic modes of interference noise in specific frequency bands, including:

[0021] Real-time acquisition of spindle speed and hydraulic pressure parameters of tunneling equipment as multi-source operating condition parameters;

[0022] The amplitude and dominant frequency characteristics of interference noise in a specific frequency band are recorded synchronously as characteristic modes;

[0023] A quantitative mapping relationship is established by fitting the mathematical relationship between multi-source operating condition parameters and characteristic modes using the least squares method.

[0024] Furthermore, a quantitative mapping relationship is established by fitting the mathematical relationship between multi-source operating condition parameters and characteristic modes using the least squares method, including:

[0025] Multi-source operating condition parameters are used as the input matrix, and characteristic modes are used as the observation matrix;

[0026] Construct a system of linear equations and solve the coefficient matrix using the least squares criterion;

[0027] The coefficient matrix obtained by the solution is used to construct a quantitative mapping relationship from multi-source operating condition parameters to characteristic modes.

[0028] Furthermore, based on the quantitative mapping relationship, the period in which the interference noise dominates is determined as the noise reference period, and the coherence between multi-channel vibration signals within the noise reference period is calculated to establish a spatial coherence feature template for the noise, including:

[0029] Based on the quantitative mapping relationship, the time period in which the multi-source operating condition parameter values ​​reach the preset threshold is selected as the noise reference time period;

[0030] Extract multi-channel vibration signal data within the noise reference time period;

[0031] Calculate the cross-correlation matrix between the vibration signals of each channel;

[0032] A spatial coherence feature template for noise is constructed based on the cross-correlation matrix.

[0033] Further, the cross-correlation matrix between the vibration signals of each channel is calculated, including:

[0034] Select vibration signal data segments for each channel within the noise reference period;

[0035] Calculate the normalized cross-correlation coefficient between every two different vibration signal data segments from different channels;

[0036] Arrange the cross-correlation coefficients of all channel pairs in matrix form according to channel order to generate a cross-correlation coefficient matrix.

[0037] Furthermore, the vibration signal is compared with a spatial coherence feature template to suppress signal components with high spatial coherence, thereby obtaining a purified effective seismic wave signal, including:

[0038] Calculate the similarity coefficient between each data segment in the vibration signal and the spatial coherence feature template;

[0039] Data segments with similarity coefficients exceeding a preset threshold are marked as high spatial coherence signal components;

[0040] Data segments marked as having high spatial coherence signal components are attenuated.

[0041] The attenuated data segments and the retained data segments are reconstructed into a purified effective seismic wave signal.

[0042] Furthermore, based on the purified effective seismic wave signal inversion, the geological structure and stress field distribution ahead of the working face are retrieved, and stress concentration areas and disaster risk areas are identified based on the stress field distribution, including:

[0043] Extract the seismic wave travel time and amplitude information from the purified effective seismic wave signal;

[0044] The morphology of the geological structural interface and stress field distribution in front of the working face were inverted using tomographic imaging algorithms.

[0045] Calculate the rate of change of stress gradient at each point in the stress field distribution;

[0046] Regions where the rate of change of stress gradient exceeds a preset threshold are marked as stress concentration areas;

[0047] Disaster risk areas are determined by combining the geological structure interface morphology and the distribution of stress concentration zones.

[0048] Furthermore, tomographic imaging algorithms were used to invert the morphology of the geological structural interface and stress field distribution ahead of the working face, including:

[0049] The rock mass region in front of the working face is discretized into multiple grid cells;

[0050] Linear inversion equations were constructed using seismic wave travel time and amplitude information as observation data.

[0051] The physical property parameters of the mesh elements are solved by an iterative reconstruction algorithm;

[0052] The geological structure interface morphology and stress field distribution are determined based on the distribution of physical property parameters.

[0053] Compared with the prior art, the present invention has the following beneficial effects:

[0054] 1. Effectively improves the accuracy and reliability of geological exploration and geostress inversion during coal mine tunneling. By establishing a quantitative mapping relationship between multi-source operating parameters and specific frequency band interference noise characteristic modes, it is possible to accurately identify the strong interference noise characteristics generated by the operation of tunneling equipment, and accurately determine the noise-dominant period based on this mapping relationship. Furthermore, it constructs a spatial coherence characteristic template for the noise, enabling targeted separation of interference components with high spatial coherence to the noise template. This significantly improves the signal-to-noise ratio of the effective seismic wave signal, overcoming the technical challenge of traditional methods in handling strong interference noise in the same frequency band, and providing a high-quality data foundation for subsequent geological inversion.

[0055] 2. Inversion of geological structures and stress field distribution based on purified effective seismic wave signals can yield more accurate information on the rock mass structure and stress state distribution ahead of the working face. By combining the stress gradient change rate with the morphological characteristics of geological structural interfaces, the spatial distribution range and risk level of stress concentration zones can be accurately identified, effectively improving the early warning capability for mine disasters such as rockbursts and roof falls. This forms a complete technical system from data acquisition and noise suppression to geological inversion and risk identification, providing reliable technical support for safe and efficient coal mining. Attached Figure Description

[0056] Figure 1 This is a flowchart of a method for advanced geological exploration and geostress inversion in coal mine tunneling, according to the present invention. Detailed Implementation

[0057] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0058] Example: Figure 1 This invention presents a method for advanced geological exploration and inversion of geostress in coal mine tunneling roadways, comprising the following steps:

[0059] S1. Vibration signals generated by the rock mass in front of the tunnel face are collected by fiber optic sensors installed on the tunnel wall;

[0060] S2. Perform spectrum analysis on the collected vibration signals to identify specific frequency band interference noise generated by the operation of the tunneling equipment;

[0061] S3. Real-time monitoring of multi-source operating parameters of tunneling equipment, and establishment of a quantitative mapping relationship between multi-source operating parameters and characteristic modes of interference noise in a specific frequency band;

[0062] S4. Based on the quantitative mapping relationship, determine the period when the interference noise is dominant as the noise reference period, and calculate the coherence between the multi-channel vibration signals within the noise reference period to establish the spatial coherence feature template of the noise.

[0063] S5. The vibration signal is compared with the spatial coherence feature template, and the signal components with high spatial coherence are suppressed to obtain the purified effective seismic wave signal.

[0064] S6. Based on the purified effective seismic wave signal, invert the geological structure and stress field distribution in front of the working face, and identify stress concentration areas and disaster risk areas based on the stress field distribution.

[0065] The first step in implementation is the deployment of fiber optic sensors. Specifically, multiple single-mode fiber optic sensors are fixedly installed at equal spatial intervals inside the rock mass of the stable section of the roadway roof and sidewalls behind the coal mine roadway excavation face. The spacing is determined based on seismic wave wavelength theory, typically choosing an interval of 5 to 10 meters. For example, an 8-meter interval can be chosen as a specific implementation interval. This distance range ensures effective spatial sampling of the vibration wave field and avoids spatial aliasing. During installation, holes are first drilled in the rock wall, with a depth of approximately 0.5 meters. Then, the fiber optic sensors, encapsulated in armored sleeves, are inserted into the holes and coupled and fixed using epoxy resin anchoring agents to ensure good mechanical coupling between the sensors and the rock mass, enabling effective detection of minute vibrations propagating within the rock mass.

[0066] After sensor deployment, the vibration signal acquisition phase begins. A high-precision distributed fiber optic sensor demodulator, employing phase-sensitive optical time-domain reflectometry (OTDR) technology, continuously scans and measures the deployed fiber optic sensor network. The demodulator emits 1560 nm wavelength laser pulses into the optical fiber and receives the backscattered Rayleigh light signal. Vibration information at various points along the fiber is sensed by demodulating the phase changes of the backscattered light. During acquisition, the sampling frequency is set to 2000 Hz. This frequency value is determined based on the Nyquist sampling theorem and considers the noise of the tunneling equipment and the highest frequency components of the effective seismic wave signal. For example, considering that the vibration frequency generated by the impact of the tunneling machine's cutting teeth on the rock mass typically does not exceed 1000 Hz, a sampling frequency of 2000 Hz meets the sampling requirements and avoids frequency aliasing. The spatial sampling interval is set to 1 meter, meaning one vibration measurement point is provided per meter of fiber length. The continuously acquired vibration signals comprise two main types: one is sound waves generated by mechanical vibrations during tunneling equipment operation, propagating through the rock mass. These waves have a high frequency and concentrated energy, primarily distributed in the 300 Hz to 800 Hz range. The other type is seismic waves generated by stress changes in the rock mass ahead of the working face or reflections from geological structural interfaces. These waves have a lower frequency and wider bandwidth, primarily distributed in the 20 Hz to 200 Hz range. All acquired raw vibration data is transmitted in real time to the data processing unit for storage, providing fundamental data for subsequent analysis and processing.

[0067] When acquiring vibration signals, it is crucial to ensure the stability of the measurement environment. Temperature variations within the tunnel should be controlled within ±5 degrees Celsius to avoid impacting the accuracy of the fiber optic sensor measurements. Simultaneously, the laser output power of the fiber optic demodulator should be kept stable within a range of 10 milliwatts ±0.5 milliwatts to ensure the signal-to-noise ratio (SNR) of the measured signal meets requirements. During data acquisition, signal quality indicators such as SNR and dynamic range are monitored in real time. When the SNR falls below 30 dB, a data re-acquisition mechanism is automatically triggered. The acquired vibration data is stored in binary format, with each data file containing a timestamp, channel number, sampling frequency, and vibration amplitude to ensure data integrity and traceability.

[0068] When performing spectral analysis on the acquired vibration signals, the vibration signals continuously acquired via fiber optic sensors in a distributed measurement manner are first read from the data storage unit. The read vibration signals contain time-domain data from multiple channels, with each channel corresponding to a measurement point of a fiber optic sensor. Before performing the Fast Fourier Transform, the vibration signals are preprocessed. Preprocessing includes removing DC components and trend terms, and using linear interpolation to fill in a small number of missing data points caused by transmission interference. For example, interpolation is performed when there are fewer than 5 consecutive missing data points; when there are too many missing data points, the data for that time period is marked as invalid. A Hanning window function is applied to the preprocessed vibration signals, with a window length of 1024 sampling points and an overlap rate of 50% to reduce spectral leakage effects.

[0069] The windowed time-domain vibration signal was converted to the frequency domain using a Fast Fourier Transform (FFT) algorithm to obtain a spectrum. The FFT was set to 2048 points, with a frequency resolution of approximately 0.98 Hz, which is sufficient to identify specific frequency band interference noise generated by the tunneling equipment. The transformed spectrum is expressed as power spectral density in volts squared per Hz, with the horizontal axis representing frequency and the vertical axis representing power spectral density. The calculated spectrum contains both real and imaginary parts; the amplitude spectrum is obtained after taking the modulus value for subsequent analysis.

[0070] When identifying continuous spectral peaks in the frequency spectrum that match the operating frequency characteristics of the tunneling equipment, the background noise level is first determined. The background noise level is obtained by calculating the average power spectral density in the 50 Hz to 100 Hz band of the spectrum, which is typically less affected by equipment operation. Spectral peaks with power spectral density values ​​exceeding three times the background noise level are initially identified as potential interference peaks. These potential interference peaks are then assessed for continuity, requiring that continuous peaks maintain high power spectral density values ​​at at least five adjacent frequency points; for example, the power spectral density values ​​at adjacent frequency points should not fluctuate by more than 30%. Simultaneously, these peaks are checked for harmonic characteristics, i.e., the presence of multiple peaks with frequencies that are integer multiples of each other, a typical characteristic of tunneling equipment vibration.

[0071] When determining the frequency band boundary of interference noise in a specific frequency band based on the distribution range of continuous spectral peaks, first identify all spectral peaks that meet the continuity condition and record the center frequency and half-power bandwidth of each peak. Multiple spectral peaks with center frequencies differing by less than 20 Hz are grouped into the same interference frequency band. Calculate the lowest and highest frequencies of this interference frequency band. The lowest frequency is taken as the lowest center frequency in the group of peaks minus half its half-power bandwidth, and the highest frequency is taken as the highest center frequency in the group of peaks plus half its half-power bandwidth. For example, when a group of harmonic spectral peaks with center frequencies of 125 Hz, 250 Hz, and 375 Hz are identified, their half-power bandwidths are 15 Hz, 18 Hz, and 20 Hz, respectively. Therefore, the determined lower boundary of the interference frequency band is 125 - 15 / 2 = 117.5 Hz, and the upper boundary is 375 + 20 / 2 = 385 Hz. The finally determined frequency band boundary of the specific frequency band interference noise is used for subsequent quantitative mapping relationships and noise suppression processing.

[0072] When monitoring multi-source operating parameters of the tunneling equipment in real time, sensors installed on the equipment collect spindle speed and hydraulic pressure parameters. The spindle speed sensor uses a non-contact photoelectric encoder, installed on the drive shaft of the tunneling machine's cutting section, with a measurement range of 0-100 rpm, a sampling frequency of 10 Hz, and a measurement accuracy of ±0.5 rpm. The hydraulic pressure sensor uses a piezoresistive pressure transmitter, installed in the main oil circuit of the tunneling machine's hydraulic system, with a measurement range of 0-40 MPa, a sampling frequency of 10 Hz, and an accuracy class of 0.5. The collected multi-source operating parameters are transmitted to the data processing unit via an industrial fieldbus, maintaining time synchronization with the vibration signal acquisition system, with a time synchronization accuracy controlled within 10 milliseconds.

[0073] When synchronously recording the amplitude and dominant frequency characteristics of interference noise in a specific frequency band as feature modes, the real-time acquired vibration signal is first bandpass filtered based on the previously determined frequency band boundaries of the interference noise. The bandpass filter is a finite-length unit impulse response digital filter with an order of 100, a passband ripple of less than 0.1 dB, and a stopband attenuation greater than 60 dB. The filtered signal is used to calculate feature values ​​every 10 seconds. The amplitude feature is taken as the root mean square value of the vibration signal within that time period, calculated by averaging the squared signal and then taking the square root, with the unit being volts, consistent with the original vibration signal. The dominant frequency feature is obtained by finding the frequency value corresponding to the maximum power spectral density in the spectrum of that time period, with a frequency resolution of 0.98 Hz. These feature values ​​are time-aligned with the multi-source operating parameters of the corresponding time period to form a paired dataset.

[0074] When fitting the mathematical relationship between multi-source operating condition parameters and characteristic modes using the least squares method, the input matrix and observation matrix are first constructed. The input matrix consists of multi-source operating condition parameters, including spindle speed and hydraulic pressure, with 100 consecutive sampling points for each parameter, resulting in a matrix dimension of 100 rows and 2 columns. The observation matrix consists of characteristic modes, including amplitude and dominant frequency, also with a matrix dimension of 100 rows and 2 columns. To eliminate dimensional differences, all parameters are standardized by subtracting the mean and dividing by the standard deviation, ensuring that each parameter has a mean of 0 and a variance of 1.

[0075] A multiple linear regression model is established by constructing a system of linear equations and solving for the coefficient matrix using the least squares criterion. The model expression is Y = X × B + E, where Y is the observation matrix, X is the input matrix, B is the coefficient matrix to be solved, and E is the error matrix. The least squares criterion is used to find the coefficient matrix B that minimizes the sum of squared errors. The specific calculation process involves first calculating the transpose of the input matrix X multiplied by X, then finding the inverse matrix of this input matrix, multiplying it by the transpose of X, and finally multiplying it by the observation matrix Y. A regularization term is added during the calculation process, with a regularization parameter of 0.01, to avoid ill-conditioned problems during matrix inversion.

[0076] When constructing a quantitative mapping relationship from multi-source operating parameters to characteristic modes using the obtained coefficient matrix, the coefficient matrix is ​​stored as a mapping relationship model. This model contains a 2x2 coefficient matrix, as well as standardized parameters such as the mean and standard deviation of each parameter. When new multi-source operating parameters are input, the input parameters are first subjected to the same standardization process, and then the predicted characteristic mode values ​​are calculated through matrix multiplication. For example, when the spindle speed is 20 rpm and the hydraulic pressure is 25 MPa, the mapping relationship can predict that the amplitude of interference noise in a specific frequency band is approximately 0.5 volts, and the dominant frequency is approximately 150 Hz. The established quantitative mapping relationship is updated every 24 hours to adapt to changes in equipment operating conditions. A sliding window mechanism is used during updates, retaining the latest 100 sets of data and discarding the oldest data to ensure the timeliness of the mapping relationship. Simultaneously, a mapping relationship quality assessment mechanism is established. When the prediction error exceeds a threshold, the mapping relationship is re-established to ensure the accuracy and reliability of the quantitative mapping relationship.

[0077] When determining the period when interference noise dominates as the noise reference period based on the quantitative mapping relationship, the first step is to use the established quantitative mapping relationship between multi-source operating parameters and characteristic modes to predict the characteristic modes of the real-time acquired spindle speed and hydraulic pressure parameters. The preset thresholds are set based on historical data analysis. By statistically analyzing the distribution of characteristic mode values ​​over multiple normal operating cycles, the 95th percentile of the amplitude characteristic value is taken as the amplitude threshold (e.g., 0.6 volts). The 90th percentile of the dominant frequency characteristic value is taken as the dominant frequency threshold (e.g., 180 Hz). When both the predicted amplitude and dominant frequency characteristic values ​​exceed their respective thresholds, the period is determined to be the period when interference noise dominates. To ensure the reliability of the fault determination, this state must last at least 10 seconds to avoid misjudgments caused by instantaneous fluctuations. Simultaneously, the start and end timestamps of each qualified period are recorded to form a list of noise reference periods.

[0078] When extracting multi-channel vibration signal data within a noise reference period, the raw vibration data of all channels for the corresponding period are read from the data storage system according to the noise reference period list. The data extraction time range for each channel is extended by 2 seconds before and after the noise reference period to preserve the complete signal transition process. The extracted data segments are time-aligned to ensure that the vibration signals of all channels have the same time base. The extracted multi-channel vibration signal data undergoes quality checks, including checking data integrity, signal-to-noise ratio level, and outliers. When the data quality of a certain channel is substandard, interpolation compensation is performed using data from adjacent channels. For example, if data from a certain sensor is missing, the average data from the two sensors before and after it is used to supplement it. Finally, a three-dimensional data matrix is ​​obtained, containing time dimension, channel dimension, and signal amplitude dimension.

[0079] When calculating the cross-correlation matrix between vibration signals of each channel, the multi-channel vibration signal data within each noise reference time period are first preprocessed. Preprocessing includes removing linear trends and mean normalization, ensuring the mean of each channel signal is zero. The normalized cross-correlation coefficient is calculated using a time-domain method. For any two channel signal sequences, their correlation coefficient at zero time delay is calculated. Specifically, the corresponding points of the two signal sequences are multiplied and summed, then divided by the product of the standard deviations of the two signal sequences and the number of data points. Cross-correlation coefficients are calculated for all possible channel combinations. For a system with 20 channels, a total of 190 cross-correlation values ​​need to be calculated. Each cross-correlation value ranges from -1 to +1; the closer the value is to +1, the stronger the correlation between the two channel signals.

[0080] When constructing the spatial coherence feature template for noise based on the cross-correlation coefficient matrix, all calculated cross-correlation values ​​are arranged into a symmetric matrix according to channel number order. The rows and columns of this matrix correspond to channel numbers, and the diagonal elements represent the cross-correlation coefficient between each channel and itself, with a constant value of 1. The spatial coherence feature template is constructed by averaging the cross-correlation coefficient matrices calculated over multiple noise reference periods to improve template stability. For example, the cross-correlation coefficient matrices calculated over the most recent 10 noise reference periods are selected, and the average value of each matrix element is used as the final spatial coherence feature template value. The template data is stored in matrix form, simultaneously recording the noise reference period information, channel configuration information, and calculation timestamp used to construct the template. To adapt to changes in the tunnel environment, the spatial coherence feature template is updated every 8 hours using a weighted average method. The new template has a weight of 0.7, and the old template has a weight of 0.3, ensuring that the template can both track environmental changes and maintain sufficient stability. The final spatial coherence feature template is used for subsequent identification and suppression of noise components in vibration signals.

[0081] When comparing vibration signals with spatial coherence feature templates, the real-time acquired multi-channel vibration signals are first segmented into continuous data segments for processing. Each data segment is set to a length of 1024 sampling points, corresponding to approximately 0.5 seconds, with 512 sampling points overlapping between adjacent segments to ensure signal processing continuity. Each data segment undergoes preprocessing, including removing DC components and linear trends, and is windowed using a Hanning window function to reduce spectral leakage. The preprocessed data segments are then used to calculate similarity with the spatial coherence feature templates.

[0082] When calculating the similarity coefficients between each data segment of the vibration signal and the spatial coherence feature template, a similarity measurement method based on the cross-correlation coefficient matrix is ​​adopted. For each data segment, the cross-correlation coefficient matrix between its multi-channel signals is calculated, and the calculation method of this matrix is ​​consistent with the method used when constructing the spatial coherence feature template. Then, the calculated cross-correlation coefficient matrix is ​​compared element-by-element with the spatial coherence feature template, and the similarity coefficient is defined as the reciprocal of the root mean square value of the difference between corresponding elements of the two matrices. The similarity coefficient ranges from 0 to 1, and the closer the value is to 1, the higher the similarity between the data segment and the noise template. To improve computational efficiency, matrix operations are used to accelerate the calculation, and data from all channels are processed in batches.

[0083] When data segments with similarity coefficients exceeding a preset threshold are marked as high spatial coherence signal components, the preset threshold is determined based on extensive experimental data analysis. By analyzing the similarity coefficient distribution characteristics of clean seismic wave signals and noise signals under different operating conditions, the threshold is set to 0.8, which can effectively distinguish between noise-dominated signals and seismic wave-dominated signals. A sliding window mechanism is used in the marking process; only when the similarity coefficients of three consecutive data segments exceed the threshold is that time segment marked as a high spatial coherence signal component, thus avoiding mismarking caused by transient interference. The marking information is stored in binary mask form, with each data segment corresponding to a marking value for subsequent signal processing.

[0084] When attenuating data segments marked as having high spatial coherence, an adaptive filtering method is employed. The attenuation strength is determined based on the similarity coefficient; the higher the similarity coefficient, the greater the attenuation. The attenuation coefficient is calculated by subtracting the similarity coefficient from 1. For example, when the similarity coefficient is 0.9, the attenuation coefficient is 0.1, meaning the amplitude of the data segment is attenuated to 10% of its original value. Attenuation is performed in the frequency domain. After performing a Fast Fourier Transform on each data segment, the attenuation coefficient is multiplied in the frequency domain, and then an Inverse Fast Fourier Transform is performed back to the time domain. For particularly strong noise components, when the similarity coefficient exceeds 0.95, complete filtering is applied, and the data segment is set to zero.

[0085] When reconstructing the purified effective seismic wave signal from the attenuated data segments and the retained data segments, an overlap-preservation method is used. Since there is overlap between data segments during processing, redundant data in the overlapping areas needs to be removed during reconstruction. For each processed data segment, the middle 512 sampling points are taken as the effective data segment, and the overlapping area is removed. Then, the effective data segments are stitched together in chronological order to form a complete vibration signal. To ensure signal continuity, a cosine smoothing transition is applied at the stitching points, with a transition region length of 50 sampling points. The final purified effective seismic wave signal is saved as a new data file, while retaining all parameter records from the processing process, including the similarity coefficient, label status, and attenuation coefficient of each data segment, for subsequent quality assessment and source tracing analysis. The reconstructed signal has a significantly improved signal-to-noise ratio, providing a high-quality data foundation for subsequent geological inversion.

[0086] When retrieving the geological structure and stress field distribution ahead of the working face based on the purified effective seismic wave signal, the first step is to extract the travel time and amplitude information of the seismic waves from the purified effective seismic wave signal. The travel time information is obtained by identifying the arrival time of specific phases in each channel signal, such as the first arrival time of the P-wave. An automatic picking algorithm based on the energy ratio method is used, with an energy ratio threshold of 10 and a sliding window length of 50 sampling points to ensure the first arrival time picking accuracy is within 2 milliseconds. Amplitude information extraction mainly targets the maximum amplitude value of the direct wave, while also recording the amplitude attenuation characteristics. The amplitude value is taken as the peak value of the signal envelope, and the unit is volts. All travel time and amplitude data are organized into an observation dataset according to channel number and source-receiver distance for subsequent inversion calculations.

[0087] When using tomographic imaging algorithms to invert the morphology of geological structural interfaces and stress field distribution ahead of the working face, the rock mass region ahead of the working face is first discretized into multiple grid cells. The discretization range is determined according to the tunnel advancement direction and detection depth; for example, a lateral range of -50 meters to +50 meters and a longitudinal range of 0 meters to 100 meters are set, with the grid cell size set to a 2-meter × 2-meter × 2-meter cube. Each grid cell is assigned an initial wave velocity value and an attenuation coefficient value. The initial values ​​are set based on the experimental results of the physical and mechanical parameters of the rock mass in the mining area; for example, the initial wave velocity of the intact rock mass is set to 4000 meters per second, and the initial attenuation coefficient is set to 0.5 dB per meter.

[0088] When constructing linear inversion equations using seismic wave travel time and amplitude information as observation data, a travel time equation based on ray theory and an amplitude attenuation equation based on wave theory are established. The travel time equation states that the observed travel time equals the sum of the products of the slowness of each grid cell along the ray path and the path length, where slowness is the reciprocal of the wave velocity. The amplitude equation states that the logarithm of the ratio of the observed amplitude to the initial amplitude equals the negative of the sum of the products of the attenuation coefficients of each grid cell along the ray path and the path length. The two equations are combined into a matrix form, where the coefficient matrix consists of the ray path length, and the unknown vector represents the perturbation values ​​of the physical property parameters of each grid cell.

[0089] When solving for the physical property parameters of grid cells using an iterative reconstruction algorithm, a joint iterative reconstruction technique is employed. Each iteration consists of two steps: first, calculating the theoretical travel time and theoretical amplitude based on the current model; second, calculating the residuals between the observed and theoretical values; and finally, updating the model parameters using the inversion equations. The iteration termination condition is set to a residual descent rate of less than 5% or a maximum number of iterations reaching 50. To ensure inversion stability, model smoothing constraints and damping constraints are added, with a smoothing constraint weight of 0.1 and a damping constraint weight of 0.01. Finally, the wave velocity and attenuation coefficient values ​​for each grid cell are obtained, with the wave velocity value in meters per second and the attenuation coefficient in decibels per meter.

[0090] When determining the morphology of geological structural interfaces and the distribution of stress fields based on the distribution of physical property parameters, geological structural interfaces are identified by analyzing the spatial distribution characteristics of wave velocity and attenuation coefficient fields. Areas with abnormal wave velocity indicate lithological changes or structural development; for example, areas with wave velocities below 3500 m / s may correspond to fracture zones, while areas with wave velocities above 4500 m / s may correspond to intact rock masses. The stress field distribution is obtained through the correlation between wave velocity and stress, establishing a conversion relationship between changes in wave velocity and stress. The conversion coefficient is calibrated through laboratory core tests; for example, every 100 m / s change in wave velocity corresponds to a 10 MPa change in stress.

[0091] When calculating the rate of change of stress gradient at each point in the stress field distribution, the central difference method is used to calculate the stress gradient at each grid point. For each grid point, the stress difference between it and its six adjacent grid points is calculated, divided by the grid spacing to obtain the gradient value, and the maximum gradient value is taken as the rate of change of stress gradient at that point. A 3×3×3 window is used for smoothing during gradient calculation to reduce the impact of data noise. The unit of the rate of change of stress gradient is megapascals per meter, representing the amount of stress change per meter.

[0092] When marking areas where the rate of change of stress gradient exceeds a preset threshold as stress concentration zones, the preset threshold is determined based on the geological conditions of the mining area and historical microseismic monitoring data. By statistically analyzing the gradient distribution characteristics of normal areas and stress concentration zones, the threshold is set to 5 MPa per meter. When the rate of change of stress gradient at a certain grid point exceeds the threshold for three consecutive iteration steps, that grid point is marked as a stress concentration point. All adjacent stress concentration points are connected to form a region, creating a stress concentration zone distribution map. Simultaneously, the maximum gradient value, spatial extent, and geometric characteristics of each stress concentration zone are recorded.

[0093] When determining disaster risk areas by combining the morphology of geological structural interfaces with the distribution of stress concentration zones, a multi-factor overlay analysis method is employed. First, the stress concentration near geological structural interfaces is identified, especially in areas with drastic changes in interface morphology, such as fault inflection points or intersections. Then, the spatial configuration relationship between stress concentration zones and geological structures is analyzed. When a stress concentration zone is located within 50 meters of a structural interface and the stress gradient change rate exceeds 10 MPa per meter, the area is classified as a high-risk zone. Simultaneously, considering the degree of mining impact, the risk level of stress concentration zones within 100 meters directly in front of the working face is increased based on the roadway advance direction and mining history. Finally, a disaster risk zoning map is generated, dividing risk areas into high-risk, medium-risk, and low-risk levels, each marked with a different color. The spatial coordinates, risk level, and descriptions of major risk factors for each risk area are output, providing a basis for decision-making in coal mine safety production.

[0094] All calculations involved in the embodiments are dimensionless numerical calculations, and the preset parameters and thresholds in the calculations are set by those skilled in the art according to the actual situation.

[0095] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented using software, the above embodiments can be implemented, in whole or in part, in the form of a computer program product.

[0096] Those skilled in the art will recognize that the modules and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and inventive constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0097] In addition, the functional modules in the various embodiments of this application can be integrated into one processing module, or each module can exist physically separately, or two or more modules can be integrated into one module.

[0098] In the several embodiments provided in this application, it should be understood that the disclosed systems, apparatuses, and methods can be implemented in other ways. For example, the apparatus embodiments described above are merely illustrative; for instance, the division of modules is only a logical functional division, and in actual implementation, there may be other division methods. For example, multiple modules or components may be combined or integrated into another system, or some features may be ignored or not executed. Furthermore, the coupling or direct coupling or communication connection shown or discussed may be through some interfaces; the indirect coupling or communication connection between apparatuses or modules may be electrical, mechanical, or other forms.

[0099] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

[0100] In conclusion, the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for advanced geological exploration and ground stress inversion of a coal mine tunneling roadway, characterized in that, The method comprises the following steps: S1, collecting vibration signals generated by the rock mass in front of the tunneling face through optical fiber sensors arranged on the tunnel wall; S2, performing frequency spectrum analysis on the collected vibration signals to identify specific frequency band interference noise generated by the operation of the tunneling equipment; S3, real-time monitoring of multi-source working condition parameters of the tunneling equipment, and establishing a quantitative mapping relationship between the multi-source working condition parameters and the characteristic mode of the specific frequency band interference noise; S4, determining a period in which interference noise is dominant as a noise reference period based on the quantitative mapping relationship, and calculating the coherence between multi-channel vibration signals in the noise reference period to establish a spatial coherence characteristic template of the noise; S5, comparing the vibration signals with the spatial coherence characteristic template, suppressing the signal components with high spatial coherence, and thus obtaining the purified effective seismic wave signal; S6, inverting the geological structure and stress field distribution in front of the working face according to the purified effective seismic wave signal, and identifying stress concentration areas and disaster risk areas based on the stress field distribution.

2. The method according to claim 1, characterized in that, The method comprises the following steps: The optical fiber sensors are arranged in the tunnel roof and the rock walls on both sides in a spatial equidistant manner; The vibration signals are continuously collected by the optical fiber sensors in a distributed measurement manner, wherein the vibration signals contain acoustic waves and seismic waves.

3. The method according to claim 1, characterized in that, The collected vibration signals are subjected to frequency spectrum analysis to identify specific frequency band interference noise generated by the operation of the tunneling equipment, including: Performing fast Fourier transform on the vibration signals to obtain a frequency spectrum graph; Identifying a continuous spectrum peak in the frequency spectrum graph that matches the frequency characteristic of the tunneling equipment operation; Determining the frequency band boundary of the specific frequency band interference noise according to the distribution range of the continuous spectrum peak.

4. The method according to claim 1, characterized in that, Real-time monitoring of multi-source working condition parameters of the tunneling equipment, and establishing a quantitative mapping relationship between the multi-source working condition parameters and the characteristic mode of the specific frequency band interference noise, including: Real-time acquisition of the main shaft speed and hydraulic pressure parameters of the tunneling equipment as multi-source working condition parameters; Synchronously recording the amplitude and main frequency characteristics of the specific frequency band interference noise as characteristic modes; Fitting the mathematical relationship between the multi-source working condition parameters and the characteristic modes by the least square method to establish the quantitative mapping relationship.

5. The method according to claim 4, characterized in that, Fitting the mathematical relationship between the multi-source working condition parameters and the characteristic modes by the least square method to establish the quantitative mapping relationship, including: Taking the multi-source working condition parameters as an input matrix and the characteristic modes as an observation matrix; Constructing a linear equation system and solving the coefficient matrix by the least square method criterion; Using the solved coefficient matrix to construct the quantitative mapping relationship from the multi-source working condition parameters to the characteristic modes.

6. The method according to claim 1, characterized in that, Based on the quantitative mapping relationship, determining a period in which interference noise is dominant as a noise reference period, and calculating the coherence between multi-channel vibration signals in the noise reference period to establish a spatial coherence characteristic template of the noise, including: Selecting a period in which the multi-source working condition parameter value reaches a preset threshold as the noise reference period according to the quantitative mapping relationship; Extracting multi-channel vibration signal data in the noise reference period; Calculating the cross-correlation coefficient matrix between the vibration signals of each channel; Based on the cross-correlation coefficient matrix, constructing the spatial coherence characteristic template of the noise.

7. The method according to claim 6, characterized in that, Calculating the cross-correlation coefficient matrix between the vibration signals of each channel, including: Select the vibration signal data segment of each channel in the noise reference period; Calculate the normalized cross-correlation coefficient between each two different channel vibration signal data segments; Arrange the cross-correlation coefficients of all channel pairs in matrix form according to the channel order to generate the cross-correlation coefficient matrix.

8. The method according to claim 1, characterized in that, Compare the vibration signal with the spatial coherence feature template, suppress the signal components with high spatial coherence, and obtain the purified effective seismic wave signal, including: Calculate the similarity coefficient of each data segment in the vibration signal and the spatial coherence feature template; Mark the data segments with similarity coefficients exceeding the preset threshold as high spatial coherence signal components; Perform attenuation processing on the data segments marked as high spatial coherence signal components; Reconstruct the data segments after attenuation processing and the reserved data segments into the purified effective seismic wave signal.

9. The method according to claim 1, characterized in that, According to the purified effective seismic wave signal, the geological structure and stress field distribution in front of the working face are inverted, and the stress concentration area and disaster risk area are identified based on the stress field distribution, including: Extract the seismic wave travel time and amplitude information in the purified effective seismic wave signal; Use tomographic imaging algorithm to invert the geological structure interface morphology and stress field distribution in front of the working face; Calculate the stress gradient change rate of each point in the stress field distribution; Mark the area with stress gradient change rate exceeding the preset threshold as the stress concentration area; Determine the disaster risk area in combination with the geological structure interface morphology and the stress concentration area distribution.

10. The method according to claim 9, characterized in that, Use tomographic imaging algorithm to invert the geological structure interface morphology and stress field distribution in front of the working face, including: Discretize the rock mass area in front of the working face into multiple grid units; Use seismic wave travel time and amplitude information as observation data to construct linear inversion equation; Solve the physical property parameter value of the grid unit by iterative reconstruction algorithm; Determine the geological structure interface morphology and stress field distribution according to the distribution of physical property parameter value.

Images