A coal and gas outburst early warning method based on multi-dimensional time-frequency feature fusion

Abstract

The present application relates to the technical field of intelligent early warning control, in particular to a coal and gas outburst early warning method based on multi-dimensional time-frequency feature fusion, comprising: firstly collecting time sequence alignment data of cutting working conditions, acoustic emission and gas concentration; then extracting nonlinear damage features, combining transfer entropy causal measure and frequency domain span analysis to capture disaster precursors; further constructing a physical information driven rheological instability model to invert the residual activation energy of the coal rock system evolution to the critical point; finally using a multi-agent reinforcement learning network to perform spatial gridding risk assessment, outputting early warning levels and accurate coordinates of risk sources. The present application realizes accurate tracing of disaster risks through the fusion of physical laws and intelligent algorithms.

Description

Technical Field

[0001] This invention relates to the field of intelligent early warning and control technology, specifically to a coal and gas outburst early warning method based on multi-dimensional time-frequency feature fusion. Background Technology

[0002] As coal mining in China extends into deeper strata, the mechanical environment of deep tunneling faces is becoming increasingly complex, exhibiting significant characteristics such as high ground stress, high gas pressure, and strong mining disturbance. Coal and gas outbursts, as a complex mine dynamic disaster, are characterized by frequency domain drift of coal seam fracture acoustic signals, fluctuations in dynamic gas outbursts, and nonlinear reshaping of ground stress during their incubation process. To ensure inherent mine safety, various sensors, including acoustic emission, microseismic, gas concentration, and tunneling machine condition monitoring sensors, are now widely deployed at tunneling faces.

[0003] In the field of data processing, extracting effective feature information from these multi-source, heterogeneous monitoring sequences with extremely wide sampling frequencies is a key focus of current disaster prevention and mitigation research. Existing pattern recognition methods typically utilize physical quantity sequences acquired by various sensors to perform time-domain or frequency-domain index calculations, and identify disaster risks by setting critical thresholds or constructing statistical discrimination models. However, in deep mining environments, the frequent operation of tunneling machinery generates strong mechanical vibration noise, and the failure process of coal and rock masses in tectonic stress zones exhibits significant time dependence and spatial nonlinear characteristics.

[0004] From the perspective of pattern recognition and digital signal processing, these monitoring data not only contain single-dimensional changes in physical quantities, but also imply dynamic coupling relationships between multiple physical fields. How to effectively mine these heterogeneous data in depth to capture key precursor information that can objectively reflect the evolution of the coal and rock system from steady state to unsteady state is a prerequisite for achieving accurate disaster early warning.

[0005] Therefore, the urgent technical problem to be solved in this field is how to extract high-dimensional features that can accurately map the catastrophic process of the coal and rock system from multi-source heterogeneous monitoring data under the environment of high-intensity mechanical disturbance and complex stress field coupling in deep tunneling face, and to construct an early warning and judgment mechanism that can characterize the instability logic of the system.

[0006] To address this, a coal and gas outburst early warning method based on multi-dimensional time-frequency feature fusion is proposed. Summary of the Invention

[0007] The purpose of this invention is to provide a coal and gas outburst early warning method based on multi-dimensional time-frequency feature fusion. By integrating physical laws and intelligent algorithms, it achieves accurate source tracing of disaster risks. The method includes: firstly, collecting time-series aligned data on cutting conditions, acoustic emission, and gas concentration; then, extracting nonlinear damage features and combining transfer entropy causal measures and frequency domain span analysis to capture early warning signs; next, constructing a physical information-driven rheological instability model to invert the remaining activation energy of the coal-rock system evolving to the critical point; and finally, using a multi-agent reinforcement learning network to perform spatial gridded risk assessment, outputting the early warning level and precise coordinates of the risk source.

[0008] To achieve the above objectives, the present invention provides the following technical solution: A coal and gas outburst early warning method based on multi-dimensional time-frequency feature fusion includes: Acquire the cutting head rotation frequency, the instantaneous power sequence of the cutting motor, the cutting head advance coordinates, the broadband acoustic emission signal, and the gas concentration sequence; Using the rotation frequency of the cutting head as the fundamental frequency, the ratio of the sideband modulation component to the fundamental frequency energy is extracted from the broadband acoustic emission signal to obtain the nonlinear modulation coefficient sequence; Sliding window sampling is performed on the nonlinear modulation coefficient sequence, the autocorrelation coefficient and variance evolution value of adjacent window sequences are calculated, the monotonically increasing trend of the autocorrelation coefficient and the surge point of the variance evolution value are identified, and a critical instability physical triggering operator is generated. The transfer entropy algorithm is used to perform directed correlation analysis on the nonlinear modulation coefficient sequence and the gas concentration sequence, and outputs the causal correlation weight vector; the instantaneous frequency domain span feature of the spectral centroid trajectory of the broadband acoustic emission signal drops sharply to the low frequency band is extracted; By inputting the instantaneous power sequence, critical instability physical triggering operator, causal correlation weight vector, and instantaneous frequency domain span characteristics into the rheological instability model, the remaining activation energy of the coal-rock system evolving to the instability critical point is obtained. A multi-agent reinforcement learning network is constructed, using the decay rate of the remaining activation energy as the reward function, and combined with the cutting head advance coordinates to output the early warning level and location of the early warning area for coal and gas outbursts.

[0009] Preferably, the process of acquiring the cutting head rotation frequency, the instantaneous power sequence of the cutting motor, the cutting head advance coordinates, the broadband acoustic emission signal, and the gas concentration sequence includes: real-time retrieval of the inverter output parameters through the communication interface of the tunneling machine's electrical control system to acquire the cutting head rotation frequency and the instantaneous power sequence of the cutting motor; acquisition of the cutting head advance coordinates dynamically changing with the tunneling machine's advance through a displacement monitoring device arranged on the tunneling machine body; acquisition of broadband acoustic emission signals through piezoelectric acoustic emission sensors fixed to the inner wall of the roadway or in the advance borehole, and analog-to-digital conversion processing through a high-speed data acquisition card; acquisition of the gas concentration sequence through a gas sensor arranged in the return air area of ​​the working face; and connection of all acquired data to a unified clock source server, adding a global timestamp to each data stream to achieve time-series alignment acquisition of multi-source heterogeneous data.

[0010] Preferably, the process of obtaining the nonlinear modulation coefficient sequence includes: synchronously framing the broadband acoustic emission signal using the rotation frequency of the cutting head, and extracting acoustic emission time-domain sample segments corresponding to the cutting period; performing spectral transformation processing on the acoustic emission time-domain sample segments to identify the fundamental frequency energy distribution with the rotation frequency of the cutting head as the center frequency; using narrowband bandpass filtering to obtain symmetrical sideband components distributed on both sides of the center frequency, wherein the symmetrical sideband components are generated by the nonlinear modulation effect of the internal fissures of the coal and rock mass on the cutting vibration; integrating the energy at the center frequency to obtain the total fundamental frequency energy, and simultaneously integrating the energy at the symmetrical sideband components to obtain the total sideband modulation energy; calculating the ratio of the total sideband modulation energy to the total fundamental frequency energy, and mapping the calculation result with the time axis to generate the nonlinear modulation coefficient sequence.

[0011] Preferably, the process of generating the critical instability physical trigger operator includes: setting a fixed window width and sliding step size, truncating the nonlinear modulation coefficient sequence to obtain overlapping adjacent window sub-sequences; calculating the autocorrelation coefficient of each sub-sequence under first-order lag, and calculating the variance of each sub-sequence relative to its own mean, constructing a sequence of autocorrelation coefficient evolution over time and a sequence of variance evolution, respectively; performing linear regression fitting on the sequence of autocorrelation coefficient evolution over time, and calculating the fitting slope; when the fitting slope is continuously positive within a preset number of consecutive windows, and the autocorrelation coefficient value of the current window enters a preset high-value threshold range, it is determined that there is a monotonically increasing trend; calculating the rolling mean and rolling standard deviation of the variance evolution sequence before the current window; when the variance value of the current window exceeds the sum of the rolling mean and the rolling standard deviation by a preset multiple, it is determined to be the surge point; when the monotonically increasing trend and the surge point are identified simultaneously, the autocorrelation coefficient and variance value of the current window are multiplied to generate the critical instability physical trigger operator.

[0012] Preferably, the process of outputting the causal correlation weight vector includes: performing symbolization processing on the nonlinear modulation coefficient sequence and the gas concentration sequence respectively, transforming the continuous numerical sequence into a discrete state sequence reflecting the fluctuation trend; calculating the forward transfer entropy from the nonlinear modulation coefficient state sequence to the gas concentration state sequence and the backward transfer entropy from the gas concentration state sequence to the nonlinear modulation coefficient state sequence based on a preset time lag order; calculating the difference between the forward transfer entropy and the backward transfer entropy to obtain a net transfer entropy value used to characterize the asymmetry of information transmission; identifying the current dynamic dominant factor according to the positive or negative polarity of the net transfer entropy value; if the net transfer entropy value is positive, it is determined to be a stress-driven risk; if the net transfer entropy value is negative, it is determined to be a gas-driven risk; performing normalization processing on the forward transfer entropy and the backward transfer entropy, and combining the normalized value with the polarity identifier of the dynamic dominant factor to generate a causal correlation weight vector.

[0013] Preferably, the process of extracting the instantaneous frequency domain span feature of the spectral centroid trajectory of the broadband acoustic emission signal plunging to the low-frequency band includes: converting the broadband acoustic emission signal into a two-dimensional power spectral density matrix that evolves over time using a short-time Fourier transform; performing a weighted average calculation of the frequency components and corresponding energy amplitudes for each time slice in the two-dimensional power spectral density matrix to extract a spectral centroid trajectory sequence reflecting the drift characteristics of the energy concentration region; performing a first-order gradient operation on the spectral centroid trajectory sequence to identify the instantaneous moment when the negative first-order gradient exceeds a preset steep descent threshold, as the trigger time point for the steep descent in the low-frequency band; extracting the average spectral centroid within a preset time period preceding the trigger time point as a high-frequency steady-state value, and extracting the spectral centroid under the instantaneous slice following the trigger time point as a low-frequency transient value; calculating the difference between the high-frequency steady-state value and the low-frequency transient value to generate the instantaneous frequency domain span feature characterizing the abrupt change in the scale of coal body fracture.

[0014] Preferably, the rheological instability model includes: Feature Space Alignment Layer: Receives instantaneous power sequence, critical instability physical triggering operator, causal correlation weight vector and instantaneous frequency domain span feature, and performs dimensionality upscaling on each heterogeneous feature, mapping it to a unified feature hidden space; Physical causal interaction layer: Using the causal correlation weight vector as the query matrix, attention weighting is performed on the instantaneous frequency domain span features, and the features are concatenated with the critical instability physical triggering operator to generate a physical feature vector representing the damage accumulation state; Time-varying rheological evolution layer: Using the physical feature vector and the instantaneous power sequence as time-series input, the strain rate evolution and damage accumulation process of coal and rock mass under truncation disturbance are simulated through a cyclic calculation unit containing coal and rock rheological constitutive constraints, and the time-varying rheological state tensor is output. Energy Residual Decoding Layer: Performs nonlinear dimensionality reduction mapping on the time-varying rheological state tensor, calculates the physical difference between the current energy accumulation state and the preset critical destruction energy threshold, and outputs the remaining activation energy of the coal-rock system evolving to the instability critical point.

[0015] Preferably, the multi-agent reinforcement learning network includes: Spatial intelligent agent mapping layer: The current tunneling area is divided into multiple spatial monitoring grids according to the cutting head advance coordinates, and a virtual monitoring intelligent agent is configured for each spatial monitoring grid to perform risk assessment actions; Physical reward calculation layer: Obtain the time series of the remaining activation energy, calculate the negative rate of change of the remaining activation energy over time, and generate a globally shared reward value characterizing the severity of instability of the coal-rock system; Attention allocation layer: The attention mechanism is used to calculate the spatiotemporal correlation between the local features extracted by each virtual monitoring agent and the global shared reward value, and the global shared reward value is decomposed into the local contribution weights corresponding to each virtual monitoring agent; Collaborative Risk Decision Layer: By integrating the prediction outputs of each virtual monitoring agent under the constraints of the local contribution weights, joint probability voting is performed to generate the final coal and gas outburst early warning level and the risk source coordinates of the corresponding spatial monitoring grid.

[0016] Compared with the prior art, the beneficial effects of the present invention are as follows: 1. By using the rotation frequency of the cutting head as the reference fundamental frequency and extracting the ratio of the sideband modulation component to the fundamental frequency energy from the broadband acoustic emission signal, a feature extraction mechanism under active stimulated detection mode is established. This mechanism transforms the mechanical operation process of the tunneling machine into a dynamic physical detection process of the coal and rock mass. The periodic disturbance of the coal mass by the cutting head is used as a "carrier" to capture the nonlinear modulation effect of the internal cracks of the coal and rock mass on the carrier. It can decouple the feature sequence that directly reflects the nonlinear dynamic characteristics of the coal and rock mass from the background of chaotic mechanical vibration, and realize the physical correlation between the monitoring signal and the mechanical operation condition.

[0017] 2. By employing the transfer entropy algorithm to perform directed correlation analysis on the nonlinear modulation coefficient and gas concentration sequence, quantitative identification of the dominant dynamic factor of the disaster was achieved. Utilizing asymmetric measurement methods from information theory, the dynamic information flow between stress evolution signals and gas outburst signals was quantified. By outputting a causal correlation weight vector, it is possible to identify whether the current instability symptoms are dominated by ground stress compression or high-pressure gas. Thus, in the complex environment of multi-physics coupling, the true driving mechanism of the system's internal evolution is revealed, providing a causal basis for the logical judgment of the nature of the disaster.

[0018] 3. By constructing a hierarchical energy release rheological instability model, multi-source feature variables are mapped to a physical feature space and rheological state inversion is performed. Utilizing a feature space alignment layer and a time-varying rheological evolution layer, the rheological mechanical process of internal damage accumulation and energy dissipation in coal-rock masses under external cutting disturbances is simulated. The final calculated "residual activation energy" transforms abstract signal fluctuations into an energy index describing the physical distance of the coal-rock system from the instability critical point. This ensures that the early warning judgment logic is based on the constitutive relationship of solid mechanics, guaranteeing that the pattern recognition results have a clear physical evolutionary direction. Attached Figure Description

[0019] Figure 1 This is a schematic diagram of a coal and gas outburst early warning method based on multi-dimensional time-frequency feature fusion according to the present invention. Figure 2 This is a schematic diagram of the process for obtaining the nonlinear modulation coefficient sequence according to the present invention; Figure 3 This is a schematic diagram of the rheological instability model of the present invention. Detailed Implementation

[0020] 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 skilled in the art without creative effort are within the scope of protection of the present invention.

[0021] Please see Figures 1 to 3 This invention provides a coal and gas outburst early warning method based on multi-dimensional time-frequency feature fusion, the technical solution of which is as follows: Example

[0022] A method for early warning of coal and gas outbursts based on multi-dimensional time-frequency feature fusion, the specific process of which is as follows: Figure 1 As shown, it includes: Acquire the cutting head rotation frequency, the instantaneous power sequence of the cutting motor, the cutting head advance coordinates, the broadband acoustic emission signal, and the gas concentration sequence; Using the rotation frequency of the cutting head as the fundamental frequency, the ratio of the sideband modulation component to the fundamental frequency energy is extracted from the broadband acoustic emission signal to obtain the nonlinear modulation coefficient sequence; Sliding window sampling is performed on the nonlinear modulation coefficient sequence, the autocorrelation coefficient and variance evolution value of adjacent window sequences are calculated, the monotonically increasing trend of the autocorrelation coefficient and the surge point of the variance evolution value are identified, and a critical instability physical triggering operator is generated. The transfer entropy algorithm is used to perform directed correlation analysis on the nonlinear modulation coefficient sequence and the gas concentration sequence, and outputs the causal correlation weight vector; the instantaneous frequency domain span feature of the spectral centroid trajectory of the broadband acoustic emission signal drops sharply to the low frequency band is extracted; By inputting the instantaneous power sequence, critical instability physical triggering operator, causal correlation weight vector, and instantaneous frequency domain span characteristics into the rheological instability model, the remaining activation energy of the coal-rock system evolving to the instability critical point is obtained. A multi-agent reinforcement learning network is constructed, using the decay rate of the remaining activation energy as the reward function, and combined with the cutting head advance coordinates to output the early warning level and location of the early warning area for coal and gas outbursts.

[0023] Furthermore, the process of acquiring the cutting head rotation frequency, the instantaneous power sequence of the cutting motor, the cutting head advance coordinates, the broadband acoustic emission signal, and the gas concentration sequence includes: real-time retrieval of the inverter's output parameters through the communication interface of the tunneling machine's electrical control system to acquire the cutting head rotation frequency and the instantaneous power sequence of the cutting motor; acquisition of the cutting head advance coordinates dynamically changing with the tunneling machine's advance through a displacement monitoring device arranged on the tunneling machine body; acquisition of broadband acoustic emission signals through piezoelectric acoustic emission sensors fixed to the inner wall of the roadway or in the advance borehole, and analog-to-digital conversion processing through a high-speed data acquisition card; acquisition of the gas concentration sequence through a gas sensor arranged in the return air area of ​​the working face; and connection of all acquired data to a unified clock source server, adding a global timestamp to each data stream to achieve time-series alignment acquisition of multi-source heterogeneous data.

[0024] The programmable logic controller (PLC) of the tunneling machine's electrical control system is connected to the communication port of the industrial control computer via shielded twisted-pair cable. An industrial Ethernet communication protocol based on Modbus-TCP or Profinet is used. The industrial control computer, acting as a client, periodically sends read commands to the inverter's register addresses. It polls the inverter's output current and frequency mapping addresses in real time to obtain the instantaneous power value of the cutting motor and the current rotation frequency of the cutting head. The acquired raw data enters the industrial control computer's memory buffer in byte stream form and is converted into a floating-point sequence by a parsing program. Specifically, the parsing program pre-configures a mapping table between Modbus register addresses and physical quantities, combines two consecutive 16-bit registers into a 32-bit single-precision floating-point number conforming to the IEEE 754 standard, and sets a linear conversion scaling factor according to the inverter's rated parameters. The parsed value is then multiplied by the corresponding scaling factor to restore the frequency value in Hertz and the power value in kilowatts.

[0025] A rope-type displacement encoder or a high-precision laser rangefinder is installed on the side of the hydraulic propulsion cylinder of the tunneling machine as a displacement monitoring device. When the tunneling machine advances forward, the encoder outputs a pulse signal or analog voltage signal that is linearly proportional to the advance distance. This signal is connected to the tunneling machine's built-in auxiliary control unit. After Gray code conversion or A / D conversion, combined with the preset geometric parameters of the tunneling machine body, the real-time three-dimensional position coordinates of the cutting head in the three-dimensional coordinate system of the tunneling roadway are calculated and uploaded to the data processing center in real time via the industrial bus. Specifically, taking the starting position of the tunneling machine's advance as the global coordinate origin, the design azimuth angle and design slope angle of the roadway given by the mine survey are obtained; the linear advance distance measured by the displacement encoder is used as the vector modulus, which is multiplied by the cosine value of the azimuth angle and the cosine value of the slope angle to obtain the components on the horizontal and vertical axes of the roadway; combined with the swing offset of the cutting head relative to the centerline of the machine body, the absolute coordinate value of the cutting head in the three-dimensional space of the roadway is calculated through coordinate rotation matrix transformation.

[0026] Piezoelectric acoustic emission sensors are installed in the pre-drilled borehole ahead of the tunnel face. The sensors are attached to the borehole wall with a coupling agent to receive elastic waves generated by coal fracturing. The weak charge signal output by the sensors first enters the preamplifier for gain compensation, and then is transmitted to a multi-channel high-speed data acquisition card through a coaxial cable. Driven by an external trigger clock, the high-speed data acquisition card performs analog-to-digital conversion according to the preset quantization bit depth, converting the continuous acoustic pressure analog waveform into a discrete digital voltage sequence, which is stored in a high-speed solid-state drive in binary file format. Specifically, the gain of the preamplifier is set to 40 dB, and the bandwidth covers 20 kHz to 500 kHz to filter out low-frequency mechanical noise; the quantization bit depth of the high-speed data acquisition card is set to 16 bits, and the sampling mode is set to synchronous sampling to ensure that there is no phase deviation between the multi-channel acoustic emission signals; at the same time, the trigger mode of the acquisition card is configured to be triggered by hardware TTL level issued by the clock server to achieve strict control of the acquisition start time.

[0027] Intrinsically safe infrared gas sensors are suspended on the return air side support or the top of the roadway in the tunneling face. The sensors use the principle of non-dispersive infrared absorption to convert the sensed gas concentration into a standard industrial current signal (such as 4-20mA) or a digital signal (such as RS485 serial port signal). This signal is aggregated through the substations of the mine monitoring system and transmitted through the fiber optic network to the data receiving gateway of the ground or underground central station. Specifically, the A / D conversion module inside the substation samples the analog current signal at a frequency of not less than 10 Hz and encapsulates the collected gas concentration value into a data frame with a serial number. It is transmitted through fiber optic Ethernet using an industrial redundant network protocol to ensure that the data is not lost or out of sequence during transmission.

[0028] Deploy a unified clock source server based on PTP (Precision Time Protocol) or NTP (Network Time Protocol) within the local area network. All industrial control computers, PLC gateways, high-speed acquisition cards, and monitoring substations connected to the network are connected to this clock source to perform microsecond-level clock synchronization. When each piece of data (frequency, power, coordinates, acoustic emission, gas) enters the data processing unit, the underlying driver immediately calls the system global clock and appends a high-precision 64-bit nanosecond-level timestamp to the header of the data packet. According to the timestamp order of each data stream, the data processing unit uses linear interpolation or a zero-order hold to resample the sequences of different sampling frequencies to a unified time step, completing the time-series alignment and encapsulation of multi-source heterogeneous data. Specifically, the unified time step is set to 10 milliseconds, that is, the characteristic alignment frequency of the entire system is 100 Hz. For acoustic emission signals with sampling frequencies much higher than 100 Hz, the energy envelope value or root mean square value within each 10-millisecond interval is first calculated in the data processing unit, and then the characteristic value is matched with the timestamp of low-frequency data such as power and gas.

[0029] By clearly defining the hardware interaction paths of the frequency converter communication interface, displacement encoder, and high-frequency acquisition card, standardized acquisition of tunneling condition parameters and multi-dimensional physical monitoring signals was achieved. The introduction of a unified clock source server with an added global timestamp resolved the alignment challenge of high-frequency acoustic emission waveforms with low-frequency gas concentration and cutting power at different sampling frequencies, ensuring high consistency of heterogeneous data streams on the time axis. This provides a reliable synchronous data benchmark for subsequent deep fusion of multi-dimensional features and causal measurement.

[0030] Further, the process of obtaining the nonlinear modulation coefficient sequence includes: synchronously framing the broadband acoustic emission signal using the rotation frequency of the cutting head, and extracting acoustic emission time-domain sample segments corresponding to the cutting period; performing spectral transformation processing on the acoustic emission time-domain sample segments to identify the fundamental frequency energy distribution with the rotation frequency of the cutting head as the center frequency; using narrowband bandpass filtering to obtain symmetrical sideband components distributed on both sides of the center frequency, wherein the symmetrical sideband components are generated by the nonlinear modulation effect of the internal fissures of the coal and rock mass on the cutting vibration; integrating the energy at the center frequency to obtain the total fundamental frequency energy, and simultaneously integrating the energy at the symmetrical sideband components to obtain the total sideband modulation energy; calculating the ratio of the total sideband modulation energy to the total fundamental frequency energy, and mapping the calculation result with the time axis to generate the nonlinear modulation coefficient sequence, the specific process being as follows: Figure 2 As shown.

[0031] Specifically, the real-time rotation frequency of the cutting head obtained from the previous steps is read, and its corresponding rotation cycle time is calculated. This rotation cycle is used as the reference window length to slide and capture the broadband acoustic emission signal. During the framing process, a preset proportion of overlap interval is set between adjacent frames to ensure the continuity of signal features in the time dimension. Each frame's captured time-domain sample segment represents the complete acoustic emission response signal generated by the coal and rock mass during one rotation of the cutting head. Specifically, the preset proportion of overlap interval is set to 50% to 75% of the current reference window length to ensure that the time-domain signal energy remains stable during the synthesis process after subsequent windowing processing, avoiding signal feature loss due to window edge attenuation.

[0032] Windowing is pre-applied to each frame of acoustic emission time-domain sample to suppress spectral leakage. Specifically, a Hanning window is used to smooth the time-domain sample, utilizing its main lobe width and side lobe attenuation characteristics to reduce spectral leakage interference. In the amplitude-frequency characteristic curve, peak retrieval is performed within a search bandwidth of ±0.5 Hz to 2 Hz, centered on the theoretically calculated rotation frequency. The spectral frequency corresponding to the maximum amplitude is locked as the actual fundamental frequency, thereby eliminating the instantaneous speed fluctuation deviation caused by changes in motor load. Subsequently, a fast Fourier transform algorithm is used to convert the signal from the time domain to the frequency domain. In the obtained amplitude-frequency characteristic curve, the cutting head rotation frequency at the current step size is used as the center reference frequency. The frequency point with the maximum amplitude is searched within a preset extremely narrow frequency band near the center frequency and identified as the fundamental frequency component. This step ensures accurate locking of the actively detected carrier signal even under conditions of fluctuating tunneling machine speed by dynamically tracking changes in rotation frequency.

[0033] After determining the center position of the fundamental frequency component, a sideband offset is preset based on the nonlinear characteristics of the coal and rock mass. Specifically, based on the nonlinear resonance characteristics of the coal and rock mass under truncated compression, the sideband offset is set to a certain proportion or fixed step size of the fundamental frequency. Specifically, characteristic frequencies offset by 1 Hz to 5 Hz to the left and right of the fundamental frequency component are extracted. This offset range covers the low-frequency envelope modulation characteristics generated by micro-fractures within the coal and rock mass under active excitation. A digital narrowband bandpass filter is used to extract frequencies distributed around the center frequency of the fundamental frequency. The symmetrical frequency band signals on both sides, namely the left band component and the right band component, are specifically defined by the passband width of the digital narrowband bandpass filter being set to 1 Hz to 2 Hz, with its center frequency strictly aligned with the center positions of the left and right band components to filter out incoherent random noise in adjacent areas of the side bands and ensure the purity of the extracted components. These two side band components contain nonlinear modulation information caused by the development of fractures in the coal and rock mass. By extracting the spectral amplitude values ​​at the corresponding index positions in the frequency domain matrix, the decoupling and separation of the side band components and the fundamental frequency components are completed.

[0034] For the identified fundamental frequency component band, the amplitude of the spectral lines within its corresponding frequency range is summed to obtain the total fundamental frequency energy within that period. Specifically, the corresponding frequency range is defined as the frequency width extending left and right from the center spectral line of the corresponding component to a point where the amplitude decreases by 3 dB. By summing the amplitudes of all spectral lines within this effective frequency band, an energy index reflecting the true physical intensity of the component is obtained. Simultaneously, the same summing operation is performed on the amplitudes of the spectral lines within the left and right sideband component bands to obtain the total sideband modulation energy. Subsequently, the ratio of the total sideband modulation energy to the total fundamental frequency energy is calculated through division. This ratio reflects the degree of nonlinear damage to the coal and rock medium at the current moment. Specifically, based on the principle of nonlinear acoustic modulation, the breathing effect generated at the fracture interface during the damage evolution of the coal and rock mass is utilized. By quantifying the proportion of sideband component energy induced by truncating vibration in the acoustic emission signal, a quantitative characterization of the degree of nonlinear damage to the coal and rock medium is achieved.

[0035] As the cutting head continues to advance and rotate, the system cyclically executes the aforementioned framing, transformation, extraction, and calculation processes. The energy ratio results calculated for each framing window are arranged on the timeline according to their corresponding timestamps, thereby constructing a complete sequence of nonlinear modulation coefficients. This sequence serves as a feature input reflecting the evolution of the internal physical state of the coal and rock mass and is transmitted to subsequent logical judgment and decision-making stages.

[0036] By using the rotation frequency of the cutting head as an active excitation source to perform frame-segmentation and modulation analysis on the acoustic emission signal, the nonlinear characteristics generated by the breathing effect of micro-fractures inside the coal and rock mass were extracted. The amplitude interference caused by the power fluctuation of the tunneling machine was eliminated by using the energy ratio of the side band to the fundamental frequency. Nonlinear modulation coefficients that can quantitatively characterize the degree of damage to the medium were extracted, providing a data foundation with causal logic for subsequent identification of the physical process of the coal and rock system evolving from steady state to unsteady state.

[0037] Further, the process of generating the critical instability physical trigger operator includes: setting a fixed window width and sliding step size, truncating the nonlinear modulation coefficient sequence to obtain overlapping adjacent window sub-sequences, and determining the total number of data points in the sampling window width by multiplying the preset observation duration by the resampling frequency in the preceding step. The observation duration must satisfy a complete rheological cycle covering the acoustic emission characteristics of coal and rock fractures, typically set between ten and sixty seconds. By proportionally correlating the sampling window width with the resampling frequency, sufficient sample capacity is ensured within each sliding window, effectively mitigating numerical fluctuations caused by sensor random noise during statistical calculations and improving the stability of the autocorrelation coefficient and variance evolution sequence; calculating the autocorrelation coefficient of each sub-sequence under first-order lag and calculating the variance of each sub-sequence relative to its own mean, constructing the autocorrelation coefficient evolution sequence and variance evolution sequence over time respectively; performing linear regression fitting on the autocorrelation coefficient evolution sequence over time and calculating the fitting slope; when When the fitting slope remains positive for a consecutive preset number of windows, and the autocorrelation coefficient value of the current window enters a preset high-value threshold range, it is determined that there is a monotonically increasing trend; the rolling mean and rolling standard deviation of the variance evolution sequence before the current window are calculated; when the variance value of the current window exceeds the sum of the rolling mean and the rolling standard deviation by a preset multiple, it is determined to be the surge point; when the monotonically increasing trend and the surge point are identified simultaneously, the autocorrelation coefficient and variance value of the current window are multiplied to generate the critical instability physical trigger operator.

[0038] The nonlinear modulation coefficient sequence generated in the previous steps is stored in a continuous memory array. Based on the time scale of coal and rock mass fracturing evolution, a fixed sampling window width (e.g., containing 500 data points) and sliding step size (e.g., sliding once every 50 data points) are preset. Through a cyclic iterative algorithm, the starting index is moved in the memory array according to the step size. Each time, a data slice of the window width length is extracted to form a series of subsequences with overlapping regions. These subsequences are then stored sequentially in a temporary rolling buffer.

[0039] For each truncated subsequence, the autocorrelation coefficient under first lag is calculated. Specifically, the last element of the current subsequence is removed to form a preorder vector, and the first element is removed to form a postorder vector. The Pearson correlation coefficient between these two vectors is calculated to obtain the first-order autocorrelation value of the window. Simultaneously, the arithmetic mean of all elements in the subsequence is calculated, and the sum of squares of the differences between each element and the mean is obtained. The variance value is then divided by the window length. The calculated autocorrelation value and variance value are stored in the corresponding evolutionary sequence container.

[0040] Set an observation length (e.g., the most recent 10 sliding windows), extract the corresponding numerical set from the autocorrelation evolution sequence, and use the least squares method to perform linear regression fitting on this set of values. Calculate the slope parameter of the fitted line; specifically, while calculating the slope, calculate the coefficient of determination R during the fitting process. 2 Or the p-value test statistic; only when the slope is greater than zero and R0 2 Only when the autocorrelation coefficient exceeds a preset significance threshold is the increase considered to have a statistically significant trend, thus eliminating the interference of instantaneous numerical drift under non-stationary signal background. The system monitors the positive and negative polarities of the slope in real time. If the slope remains positive within a preset number of consecutive windows (e.g., 5 consecutive slides), and the autocorrelation coefficient value at the current moment exceeds the preset high value threshold, the logic judgment unit outputs a valid identifier bit for a monotonically increasing trend. Specifically, the high value threshold is obtained by probabilistic statistics on the distribution of autocorrelation coefficients of the tunneling face under normal tunneling conditions (without signs of disaster), and is usually taken as 1.2 to 1.5 times the upper quartile of the normal fluctuation range. Through this dynamic benchmark calibration, it is ensured that the warning threshold can adapt to the background signal characteristics under different geological structural zones.

[0041] Maintain a dynamic list recording the historical characteristics of the variance evolution sequence. Calculate the rolling mean and rolling standard deviation within a preset segment length before the current window. Specifically, the preset segment length (baseline window) should be significantly larger than the sampling window width, set to 5 to 10 times the sampling window width. Using this long-term historical variance level as a static background, compare the instantaneous variance jump within the current short-term analysis window to highlight the contrast effect of "surge". By setting a threshold judgment coefficient, the rolling mean plus a preset multiple (such as 3 times) of the rolling standard deviation is used as a dynamic threshold. The instantaneous variance value calculated in the current window is compared with this dynamic threshold in real time. Once the instantaneous value exceeds the threshold, it is determined to be a surge point in the variance evolution process, and a surge point trigger flag is output synchronously.

[0042] An internal synchronous logic gate is constructed to monitor the status of the monotonically increasing trend flag and the surge trigger flag in real time. When both flags are in a valid state at the same time or within a preset tolerance time window (specifically, the tolerance time window is set to a duration of 3 to 5 sliding steps) to compensate for the asymmetry in computational delay between the trend recognition algorithm and the surge recognition algorithm, it ensures that when the system physically enters the critical instability zone, the features of the two dimensions can achieve strong logical coupling output, indicating that the system has entered the collaborative stage of "critical slowing down" and "fluctuation enhancement". At this time, the autocorrelation coefficient value and variance value of the current window are read, a product operation is performed, the critical instability physical trigger operator is generated, and it is used as a high-dimensional feature input to the subsequent catastrophic dynamics analysis module.

[0043] By analyzing the evolution of autocorrelation coefficients and variances using a sliding window, the physical signs of "critical slowdown" before instability in coal-rock systems were captured. A collaborative logic combining linear regression slope determination and dynamic threshold surge identification was introduced to eliminate spurious triggering interference from single indicators under complex background noise. The triggering operator synthesized through product operations quantified the degree of decline in the system's ability to return to steady state and the surge in fluctuation intensity, providing a basis for accurately locating critical instability.

[0044] Further, the process of outputting the causal correlation weight vector includes: performing symbolization processing on the nonlinear modulation coefficient sequence and the gas concentration sequence respectively, transforming the continuous numerical sequence into a discrete state sequence reflecting the fluctuation trend; calculating the forward transfer entropy from the nonlinear modulation coefficient state sequence to the gas concentration state sequence and the backward transfer entropy from the gas concentration state sequence to the nonlinear modulation coefficient state sequence based on a preset time lag order; calculating the difference between the forward transfer entropy and the backward transfer entropy to obtain a net transfer entropy value used to characterize the asymmetry of information transmission; identifying the current dynamic dominant factor according to the positive or negative polarity of the net transfer entropy value; if the net transfer entropy value is positive, it is determined to be a stress-driven risk; if the net transfer entropy value is negative, it is determined to be a gas-driven risk; performing normalization processing on the forward transfer entropy and the backward transfer entropy, and combining the normalized value with the polarity identifier of the dynamic dominant factor to generate a causal correlation weight vector.

[0045] The nonlinear modulation coefficient sequence and the gas concentration sequence are obtained, and a discretization mapping rule is set for each sequence. Specifically, the mean and standard deviation of the sequence within a moving window are calculated. The length of the moving window is set to fifty to one hundred data points, and this length must be less than the sampling window width in the previous step. Calculating the local benchmark through a shorter window makes the symbolized sequence more sensitive to transient fluctuations induced by coal and rock damage, ensuring that states 0, 1, and 2 can accurately capture marginal changes in the signal. The current value is divided into a preset number of state intervals based on its deviation from the mean. In a preferred embodiment, a three-state symbolization method is used to discretize continuous values: values ​​falling within one standard deviation of the mean are mapped to symbol zero; values ​​exceeding one standard deviation are mapped to symbol one; and values ​​below one standard deviation are mapped to symbol two. The technical rationale for adopting the aforementioned three-state partitioning rule is that this method can effectively control the size of the state space while fully preserving the three core physical trend information of the sequence—rising, stable, and falling—thus avoiding data sparsity or overfitting in subsequent transition entropy probability statistics due to an excessive number of states. Through this mapping process, the continuous floating-point sequence is converted into a symbolic dynamic sequence composed of discrete integers, thereby filtering out high-frequency random fluctuations and retaining state information reflecting physical trends.

[0046] A time lag order is defined to indicate the time lead of causal effects. Specifically, the time lag order is preferably selected within a preset range of one to five uniform resampling step sizes. This range is defined based on the macroscopic physical time lag between acoustic emission characteristics and gas concentration response, as known in the art. In specific calculations, the lag order within this preset range is traversed to select the order that maximizes mutual information and achieves optimal statistical stability as the final calculation parameter. This optimization strategy is not simply data fitting, but rather aims to accurately pinpoint the true physical lag phase that exerts the strongest causal driving effect on the dynamic evolution of gas under actual working conditions where response delays differ across multiple physical fields. When calculating the forward transfer entropy, the system slides samples on the symbol sequence, statistically analyzing the frequency of the joint event consisting of "the state of the gas concentration sequence at the next moment," "the current state of the gas concentration sequence," and "the current state of the nonlinear modulation coefficient sequence." By statistically analyzing the frequency of each combination, a conditional probability distribution model is established. Following the logic of information entropy calculation, the product of the probability and its logarithm for all possible state combinations is accumulated to obtain a numerical value reflecting the ability of stress characteristics to explain gas dynamics, i.e., the forward transfer entropy. Specifically, the logarithmic operation uses a base-2 logarithm, and the calculation result is expressed in bits. The accumulation process traverses all possible ternary joint state spaces, calculating the gas state at the current moment. The difference between the entropy under known historical conditions and the joint entropy under known historical conditions and stress characteristics quantifies the reduction in uncertainty of the future state of the gas sequence due to stress characteristics. Similarly, by switching the roles of the two sequences and repeating the above probability statistics and accumulation process, the backward transfer entropy reflecting the ability of gas dynamics to explain stress characteristics is obtained. Specifically, to ensure the statistical significance of the probability estimation, the length of the symbol sequence involved in the statistics must be at least ten times the number of state combinations, that is, the sequence length must be no less than 270 symbol points. When encountering a state combination with a frequency of zero, Laplace smoothing is used to assign a very small probability correction value to the combination to avoid mathematical overflow or calculation interruption when performing logarithmic accumulation.

[0047] The forward and backward transfer entropies are obtained and subtracted to obtain the net transfer entropy value. The system determines the dominant direction of information transmission within the system by judging the sign of the net value. If the net value is positive, it indicates that the stress characteristics contribute more to the information of the gas sequence, and the decision logic determines that the current stage is a stress-driven risk stage. If the net value is negative, it indicates that the gas sequence has a stronger leading indication of the evolution of the power system, and it is determined to be a gas-driven risk. The polarity determination result is converted into the corresponding binary identifier for storage.

[0048] Normalization is performed on the forward and backward transfer entropy, specifically by calculating the proportion of each entropy value in the sum of the two, so that the processed values ​​are distributed between zero and one. Then, the normalized forward and backward values ​​are encapsulated with the polarity identifier of the dominant dynamic factor determined in the previous step. The resulting causal correlation weight vector is output in the form of a multidimensional array, which contains the weight proportion of the contribution of different physical features to the current instability and the polarity information of the dominant dynamic, providing a basis for the parameter weighting of the subsequent rheological model.

[0049] By quantifying the information transmission intensity between nonlinear modulation characteristics and gas dynamics using the transfer entropy algorithm, the asymmetric discrimination of the driving force of coal and gas outbursts is achieved. The introduction of symbolic discretization effectively filters out random interference from the original monitoring data. The positive and negative polarities of the net transfer entropy are used to accurately locate the current driving force of the system, providing quantifiable logical support for distinguishing between the two distinct disaster modes of "stress-induced" and "gas-dominated".

[0050] Further, the process of extracting the instantaneous frequency domain span feature of the spectral centroid trajectory of the broadband acoustic emission signal plunging to the low-frequency band includes: using short-time Fourier transform to convert the broadband acoustic emission signal into a two-dimensional power spectral density matrix that evolves over time; performing a weighted average calculation of the frequency components and corresponding energy amplitudes for each time slice in the two-dimensional power spectral density matrix to extract a spectral centroid trajectory sequence reflecting the drift characteristics of the energy concentration region; performing a first-order gradient operation on the spectral centroid trajectory sequence to identify the instantaneous moment when the negative first-order gradient exceeds a preset steep descent threshold, as the trigger time point for the steep descent in the low-frequency band; extracting the average spectral centroid within a preset time period before the trigger time point as a high-frequency steady-state value, and extracting the spectral centroid under the instantaneous slice after the trigger time point as a low-frequency transient value; calculating the difference between the high-frequency steady-state value and the low-frequency transient value to generate the instantaneous frequency domain span feature characterizing the abrupt change in the scale of coal body fracture.

[0051] The time-aligned digital broadband acoustic emission signal is acquired, with a sampling frequency of 500 kHz, and windowed and framed using a 1024-point Hamming window. The overlap rate between adjacent frames is set to 50%. The time-domain discrete signal is transformed into a frequency-domain complex array using a Fast Fourier Transform. By calculating the square of the amplitude at each frequency point, a two-dimensional power spectral density matrix reflecting the energy distribution with frequency and time is obtained, providing a high-dimensional data foundation for subsequent extraction of the centroid trajectory. Specifically, in generating the power spectral density matrix, the calculated power value needs to be divided by the product of the sampling frequency and the coherence gain of the window function to normalize the power spectral density, ensuring the comparability of energy characteristics under different sampling lengths.

[0052] The algorithm iterates through each time slice in the two-dimensional power spectral density matrix to obtain the frequency vector and corresponding power value vector at the current moment. In the specific calculation, each frequency value is multiplied by its corresponding power energy value. The product sequence is accumulated within the effective bandwidth, and the accumulated result is divided by the total power energy value in the current time slice. Specifically, the lower limit of the effective bandwidth is set to 20 kHz to filter out the mechanical vibration noise interference generated by the tunneling machine cutting head through high-pass filtering characteristics; the upper limit is set to 200 kHz, which covers the main energy range of the coal body fracture sound signal, ensuring that the centroid calculation focuses on the frequency domain evolution related to fracture. The weighted average calculation is performed iteratively, and the centroid frequencies corresponding to each time step are arranged in chronological order to construct a spectral centroid trajectory sequence that reflects the dynamic shift of the energy concentration area.

[0053] A first-order backward difference operation is performed on the generated spectral centroid trajectory sequence to calculate the frequency change rate between two adjacent sampling points. A negative steep drop threshold is preset, which is set based on the standard deviation of frequency fluctuations under normal tunneling conditions. Specifically, a five-minute period during which the tunneling machine operates normally without abnormal acoustic emission signals is selected as the calibration reference period, and the standard deviation of the spectral centroid gradient sequence during this period is calculated. The steep drop threshold is then set to three to five times this standard deviation. Through this statistical dynamic doubling method, normal random frequency fluctuations are strictly distinguished from physical steep drops caused by abrupt changes in fracture scale. The current frequency gradient value is compared in real time. When the gradient value is negative and its absolute value exceeds the steep drop threshold, the current sampling sequence number is immediately locked and recorded as the instantaneous trigger time point of the low-frequency steep drop, marking the abrupt change in the coal body fracture scale from micro to macro.

[0054] Using the identified trigger time point as a reference, a preset number of sampling points (e.g., fifty data points) are traced back along the time axis. Specifically, the physical duration corresponding to the preset number of sampling points should cover a stable monitoring period of one to two seconds. Considering a 50% frame overlap rate, the number of traced points needs to be dynamically calculated based on the current frame shift length to ensure that the extracted high-frequency steady-state value can truly represent the background frequency reference before the system enters the unstable critical state. The average value of the spectral centroid of this interval is extracted as the high-frequency steady-state value characterizing the stable state before coal body rupture. Simultaneously, the instantaneous value of the spectral centroid of the slice where the trigger time point is located and its immediately adjacent slice are obtained and used as the low-frequency transient value. Specifically, in order to eliminate the interference of single-frame instantaneous noise on the span calculation, the low-frequency transient value is obtained by performing an arithmetic mean of the spectral centroids of the trigger time point and the subsequent three to five consecutive frames, thereby obtaining the relatively stable low-frequency characteristics after the sudden change moment. Through this extraction strategy of "preceding mean" and "following transient", the influence of local frequency random disturbances on the span calculation is eliminated.

[0055] The high-frequency steady-state value and the low-frequency transient value are subtracted to calculate the absolute value of the frequency drop between them. This difference is the instantaneous frequency domain span characteristic, which is used to quantify the wavelength stretching and frequency redshift caused by the large-scale fracture penetration at the moment of instability of the coal-rock system. Finally, this span characteristic is used as a key physical indicator to input into the subsequent rheological instability model to assist in the calculation of residual activation energy and risk inversion.

[0056] To ensure the robustness of the frequency domain span feature extraction algorithm under complex and non-stationary field conditions, this embodiment further includes boundary handling logic for abnormal situations: First, when multiple consecutive sampling points simultaneously meet the steep drop threshold, the system extracts only the sampling point with the largest absolute gradient value as the unique low-frequency steep drop trigger time point to avoid a single rupture event being counted repeatedly.

[0057] Second, if no sampling point that meets the steep drop threshold is detected within the current observation window, the system will force the current instantaneous frequency domain span feature to be zero, or directly maintain the output value of the previous valid trigger moment to ensure the continuity of the feature sequence.

[0058] Third, if a trigger time point is identified, but the length of the preceding traceable data is insufficient to cover the preset high-frequency steady-state extraction duration, the system will actively discard the current trigger event and continue to wait for the next valid candidate point, thereby preventing misjudgment of the steady-state benchmark due to insufficient historical background data.

[0059] Fourth, if the physical difference between the calculated high-frequency steady-state value and the low-frequency transient value is less than zero under extreme disturbance conditions, the system strictly truncates this difference to zero, thereby avoiding inputting abnormal reverse span values ​​that violate the physical law of coal body fracture 'frequency redshift' into subsequent rheological instability models. By capturing the instantaneous gradient changes in the spectral centroid trajectory, a digital mapping of the transformation process of coal fracture from "micro-cracks" to "macro-major fractures" was achieved. Using the energy centroid difference between high-frequency steady-state and low-frequency transient states, the physical redshift phenomenon caused by the surge in coal and rock fracture scale was quantified. This frequency domain span characteristic eliminates spurious interference from amplitude fluctuations, providing physical parameters that directly reflect changes in fracture scale for energy release rheological instability models.

[0060] Furthermore, the rheological instability model includes: Feature Space Alignment Layer: Receives instantaneous power sequence, critical instability physical triggering operator, causal correlation weight vector and instantaneous frequency domain span feature, and performs dimensionality upscaling on each heterogeneous feature, mapping it to a unified feature hidden space; Physical causal interaction layer: Using the causal correlation weight vector as the query matrix, attention weighting is performed on the instantaneous frequency domain span features, and the features are concatenated with the critical instability physical triggering operator to generate a physical feature vector representing the damage accumulation state; Time-varying rheological evolution layer: Using the physical feature vector and the instantaneous power sequence as time-series input, the strain rate evolution and damage accumulation process of coal and rock mass under truncation disturbance are simulated through a cyclic calculation unit containing coal and rock rheological constitutive constraints, and the time-varying rheological state tensor is output. Energy Residual Decoding Layer: Performs nonlinear dimensionality reduction mapping on the time-varying rheological state tensor, calculates the physical difference between the current energy accumulation state and the preset critical destruction energy threshold, and outputs the remaining activation energy of the coal-rock system evolving to the instability critical point. The specific process is as follows: Figure 3 As shown.

[0061] A feature alignment module is constructed in the memory of the computing environment. This module contains four sets of parallel fully connected linear mapping units. The processor inputs the instantaneous power sequence, the critical instability physical trigger operator, the causal correlation weight vector, and the instantaneous frequency domain span feature into the corresponding mapping units. Each mapping unit unifies the dimension of each heterogeneous feature to a preset high-dimensional space (e.g., 64-dimensional or 128-dimensional) by multiplying the weight matrix and accumulating the bias term. During the dimensionality upscaling process, a layer normalization algorithm is used to adjust the numerical distribution after mapping to ensure that features from different physical dimensions have consistent dimensions and distribution characteristics in the feature hiding space, thereby eliminating the impact of magnitude differences on subsequent fusion. Specifically, before entering the fully connected mapping unit, the system determines the temporal sampling frequency of each feature sequence. For instantaneous power sequences with a large number of sampling points, one-dimensional max pooling is used for downsampling. For single-valued features such as causal correlation weights, the sequence is copied to expand them to a length consistent with the spatiotemporal step size to ensure that all tensors entering the feature hiding space have the same sequence length L and hidden dimension D.

[0062] The interaction layer performs logical weighting of physical information through the dot product attention operator. The processor uses the aligned causal correlation weight vector as the query vector and the instantaneous frequency domain span feature as the key vector and value vector. By calculating the similarity score between the query vector and the key vector, a weight distribution reflecting the sensitivity of the current causal dynamics to the frequency domain evolution is generated. Specifically, after performing a dot product operation on the query vector and the key vector, the result is divided by the square root of dimension D for scaling. Then, nonlinear normalization is performed using the Softmax function, so that the sum of the weights of all frequency components equals one. This processing ensures the sparsity of the attention distribution, which can accurately lock the frequency domain span feature that contributes the most under the current causal drive. The weight distribution is applied to the frequency domain span feature to obtain the frequency domain representation after causal correction. This representation is then concatenated with the aligned critical instability physical trigger operator in the channel dimension. Finally, a single-layer nonlinear mapping function compresses the concatenated data into a fixed-length physical feature vector to represent the current instantaneous state of damage accumulation in the coal-rock system.

[0063] A gated recurrent unit (GRU) is deployed as the core computational operator in the rheological evolution layer. This operator integrates the physical feature vector, instantaneous power sequence, and rheological state stored in the hidden layer at the previous moment through internal reset and update gate structures. To introduce coal-rock rheological constitutive constraints, the weight initialization parameters of the computational unit are preset according to the parameters of the coal-rock Nishihara model measured in the laboratory. This ensures that the iterative process conforms to the physical evolution law of decay creep and stable creep of coal-rock mass under high ground stress and truncation disturbance. The output of the hidden layer at each moment constitutes the time-varying rheological state tensor. This tensor dynamically records the evolution of strain rate and the accumulation process of damage variables within the coal-rock system. Specifically, the weight initialization is achieved by mapping the long-term elastic modulus and viscous constant to the initial values ​​of the bias terms of the GRU update and reset gates through a predefined proportional mapping function. At the same time, the analytical solution generates a set of standard rheological curves as a pre-training set. Through offline learning, the initial weight matrix of the computational unit approximates the numerical solution of the rheological differential equation, thereby realizing the constraint of physical constitutive structure on the network search space.

[0064] Specifically, in order to substantially embed the coal-rock rheological constitutive constraints into the cyclic computational unit, this embodiment uses the Nishihara rheological model as the physical basis for characterizing the time-varying evolution of the coal-rock mass. This model explicitly defines the total strain of the coal-rock mass under the coupled action of constant load and disturbance as the sum of three parts: instantaneous elastic strain, viscoelastic strain, and viscoplastic strain. In the network parameter initialization stage, the system extracts the elastic modulus parameter and viscosity coefficient from the Nishihara model and maps them to the initial weights and biases of the cyclic computational unit. The specific mapping rule is to maintain a monotonically corresponding relationship between the bias term of the network update gate and the physical viscosity coefficient, and to initialize the weight matrix of the candidate state proportionally to the elastic modulus.

[0065] During model training, the system constructs a total loss function that includes physical soft constraints. This total loss function consists of two superimposed parts: the first part is the conventional mean square error loss of network prediction; the second part is the physical consistency loss. The specific calculation method of the physical consistency loss is as follows: extract the strain rate currently predicted by the network in real time, and simultaneously calculate the average of the sum of squares of the deviations between the two based on the theoretical strain rate derived by the Nishihara model under the same stress level. The system multiplies this physical consistency loss term by a preset physical constraint weight coefficient and adds it to the conventional loss term. In this way, once the predicted trajectory of the neural network deviates from the analytical range of the coal and rock constitutive equation, a huge penalty gradient will be generated, thereby forcing the model to strictly follow the physical laws of solid mechanics while fitting the monitoring data.

[0066] The physical constraint weight coefficient is specifically used to balance the model's fit to the actual on-site monitoring data and its observability to the physical laws of Nishihara rheology. In the actual model training configuration, this weight coefficient is not an arbitrary fixed constant, but is set using a dynamic adjustment strategy: In the initial stage of model training, in order to ensure that the neural network can converge to the basic distribution characteristics of the actual data first, the coefficient is preset to a small initial empirical value (e.g., 0.1); as the number of training iterations increases, the system gradually increases the weight coefficient to the set upper limit value (e.g., 0.5) through a linear increment mechanism, thereby forcibly applying strict mechanical constitutive constraints in the later stage of model convergence. The initial empirical value and the upper limit value are obtained by technicians during the offline training stage by using a hyperparameter grid search method with the goal of minimizing the comprehensive prediction error of the model on the historical validation set.

[0067] The decoding layer employs a regression mapping network composed of multiple fully connected layers. This network receives a time-varying rheological state tensor as input and maps the high-dimensional state tensor to a scalar value through layer-by-layer dimensionality reduction and nonlinear activation. This value represents the current real-time energy accumulation state of the coal-rock system. Specifically, the regression mapping network contains at least three fully connected layers with decreasing numbers of neurons, with a ReLU activation function set between each layer. At the output of the last layer, a smoothing linear unit is used to constrain the output value to a non-negative number to conform to the physical fact that the energy accumulation value is always greater than or equal to zero. The network retrieves parameters from a preset parameter storage unit that match the geological conditions of the current tunneling face. The critical failure energy threshold is used as a benchmark. Subtracting the real-time energy accumulation value from this threshold, the calculated physical difference represents the remaining activation energy of the coal-rock system evolving to the critical point of instability. Specifically, the critical failure energy threshold is obtained by performing a uniaxial compression energy evolution test on coal samples from historical outburst accident sites in the mining area, representing the instability limit of the coal-rock system's stored energy release. For coal seams of different thicknesses and different geostress levels, the corresponding static threshold is retrieved by looking up a table and corrected according to the current cutting depth. The result is output in the form of a time series, reflecting the energy safety margin of the system from physical instability and collapse.

[0068] To further clarify the above calculation process, the residual activation energy referred to in this paper specifically refers to the energy capacity margin remaining between the coal-rock system and the final instability and failure critical point under the current operating state. The specific calculation logic is as follows: the difference between the preset critical failure energy threshold and the cumulative damage energy characterization value of the system output by the rheological instability model at the current moment is obtained. The preset critical failure energy threshold is not fixed but is a baseline value determined through offline laboratory calibration based on the current coal seam burial depth, geostress level, gas pressure, and specific coal mechanical parameters.

[0069] The training method for the rheological instability model includes the following steps: First, in a laboratory environment, triaxial cyclic loading and unloading rheological tests are performed on coal and rock samples collected from the tunnel face, simultaneously collecting axial pressure, acoustic emission characteristics, and gas emission during the experiment. The physical process of the coal and rock mass from steady state to accelerated creep and eventual failure is divided into different time slices, and the true value of the "residual activation energy" corresponding to each slice is calculated and used as a regression label. Second, numerical simulation software (such as FLAC3D or Discrete Element Method) is used to generate a large number of simulated evolution sequences under different geostress gradients based on the Nishihara rheological constitutive equation. Finally, the preprocessed instantaneous power, physical triggering operators, causal weights, and frequency domain span features are aligned by timestamps to form supervised learning sample pairs containing feature vectors and corresponding physical state labels.

[0070] The model training method employs a supervised learning algorithm with nested physical information to train the rheological instability model. During initialization, laboratory-measured physical parameters are transformed into a priori distribution for the neural network. During training, the model receives input aligned to the feature space and calculates the time-varying rheological state tensor using hierarchical GRU units. The loss function consists of two parts: the first is the mean squared error loss between the predicted residual activation energy and the label value; the second is the physical consistency loss, generated by calculating whether the currently predicted strain rate evolution deviates from the analytical range of the constitutive equation, thus producing a physical constraint gradient. The optimizer updates the weight matrix through backpropagation based on the total loss value, ensuring that while fitting the data, the model's internal parameter distribution conforms to the physical laws of coal and rock rheology.

[0071] By constructing a physical information-driven rheological instability model, deep coupling between monitoring data and constitutive relations was achieved. The feature space alignment layer and the attention interaction layer worked together to effectively identify the influence weight of causal dynamics on the damage state. The time-varying rheological evolution layer introduced Nishihara model constraints to ensure that the evolution process conformed to the real physical laws of coal and rock creep. The final output residual activation energy quantified the physical boundary of the system from catastrophe.

[0072] Furthermore, the multi-agent reinforcement learning network includes: Spatial Agent Mapping Layer: Based on the cutting head's advance coordinates, the current tunneling area is divided into multiple spatial monitoring grids, and a virtual monitoring agent is configured for each spatial monitoring grid to perform risk assessment actions. Specifically, the action space of the virtual monitoring agent is defined as a set of discrete evaluation values ​​for the risk probability of the current grid (five gradient levels between 0 and 1 in this embodiment). The agent outputs the predicted action in the current state through its internal policy network, and continuously optimizes its prediction strategy using the dominant executor-critic algorithm based on the feedback of the globally shared reward value to maximize the long-term cumulative reward. Physical reward calculation layer: Obtain the time series of the remaining activation energy, calculate the negative rate of change of the remaining activation energy over time, and generate a globally shared reward value characterizing the severity of instability of the coal-rock system; Attention allocation layer: The attention mechanism is used to calculate the spatiotemporal correlation between the local features extracted by each virtual monitoring agent and the global shared reward value, and the global shared reward value is decomposed into the local contribution weights corresponding to each virtual monitoring agent; Collaborative Risk Decision Layer: By integrating the prediction outputs of each virtual monitoring agent under the constraints of the local contribution weights, joint probability voting is performed to generate the final coal and gas outburst early warning level and the risk source coordinates of the corresponding spatial monitoring grid.

[0073] The system reads the cutting head advance coordinates obtained from previous steps and constructs a three-dimensional voxelized spatial grid model in memory based on the designed cross-sectional dimensions of the tunnel and the preset advance detection depth. The active area affected by the current tunneling is divided into multiple spatial monitoring grids with equal side lengths (e.g., one meter or two meters). For each spatial monitoring grid, the system instantiates an independent virtual monitoring agent object at the software level. Each agent is associated with the local monitoring data (such as the historical cutting power and acoustic emission energy density at that location) within its corresponding grid through an index. Specifically, the system establishes a dynamic spatial mask centered on the real-time coordinates of the cutting head and uses inverse distance weighted interpolation to map the physical features acquired by each fixed monitoring point and random sensor in the tunnel to the center point of each spatial monitoring grid. This generates a state vector representing the local physical state of the area under the responsibility of each virtual monitoring agent, thus transforming the macroscopic monitoring task into a multi-point parallel local evaluation task.

[0074] The physical reward calculation layer receives the remaining activation energy sequence output by the rheological instability model in real time, performs time difference calculation, calculates the difference between the remaining activation energy at the current sampling time and the previous sampling time, and extracts its negative rate of change. This rate of change reflects the "acceleration" of the coal-rock system approaching the instability critical point. The larger the absolute value of the rate of change, the more violent the energy release and the higher the risk of system collapse. This value is set as the globally shared reward value in the reinforcement learning environment to guide all virtual monitoring agents to find the feature combination that best reflects the energy mutation through policy iteration. Specifically, the calculated negative rate of change is normalized and scaled (in this embodiment, the Tanh function is used to constrain the value between -1 and 1), and a reward smoothing coefficient is set to ensure that the reward signal has sufficient distinctiveness when energy is released rapidly, while avoiding the instability of the agent's policy learning caused by instantaneous extreme values.

[0075] To address the "credit allocation" challenge in multi-agent systems, the attention allocation layer constructs a multi-head attention scoring mechanism. Specifically, firstly, a multilayer perceptron maps the scalar global shared reward value into a feature vector consistent with the local feature dimensions of the agents. This transformation process encodes global risk information into the latent space, which then serves as a query vector in subsequent attention scoring calculations. Using the global shared reward value as the query vector and the local feature tensors extracted by each virtual monitoring agent as the key vectors, the dot product similarity between the global reward and each local feature is calculated, followed by nonlinear normalization. This yields a set of weight coefficients reflecting the contribution of each point in the space to the overall instability risk. These coefficients are the local contribution weights, which can identify which grid regions are the "risk sources" causing a sharp drop in global energy, thus decoupling the single global signal into training feedback specific to each grid.

[0076] During the decision-making phase, each virtual monitoring agent, constrained by its local contribution weight, outputs the risk probability distribution of its corresponding grid. The collaborative risk decision layer collects the outputs of all agents, constructs a joint probability matrix, and, through a weighted voting algorithm, determines the final coal and gas outburst warning level based on the highest global probability value. Simultaneously, it retrieves the agent numbers whose probability values ​​exceed a preset danger threshold and queries their corresponding spatial monitoring grid indexes to calculate the centroid coordinates of the grid in the roadway coordinate system, thus outputting the specific spatial location of the risk source. This achieves a precise transition from "whether there is risk" to "where the risk is." Specifically, the weights are directly derived from the local contribution weights generated in the preceding steps. The system multiplies the risk probability output by each agent with its corresponding local contribution weight and performs a normalized summation operation across the entire spatial grid. By assigning higher voting weight to "grids with large contributions," the final warning level is ensured to be dominated by the physical characteristics closest to the risk source.

[0077] The multi-agent reinforcement learning network dataset is constructed by creating a dynamic game dataset with spatiotemporal correlation. First, historical mine monitoring data is spatially reconstructed according to the tunneling machine's advance coordinates, transforming single-point time series into spatial feature tensors based on a three-dimensional grid. Each entry in the dataset includes: the current global spatial grid state, the local observations of each virtual agent, and the rate of change of remaining activation energy output by the rheological instability model. By augmenting historical accident cases, an interaction trajectory library is constructed, encompassing the complete dynamic process from "safety" to "criticality" and then to "instability," providing rich offline exploration samples for reinforcement learning.

[0078] The model training method employs a "centralized evaluation, distributed execution" framework to optimize the policies of the multi-agent network. In the offline simulation environment, all virtual monitoring agents simultaneously access their respective grid data to perform risk assessment actions. During training, the central evaluation network receives global state information and uses an attention allocation layer to calculate each agent's contribution score to the global shared reward (negative rate of change of remaining activation energy). Each agent's policy network updates its gradient based on its assigned local contribution weight. Through tens of thousands of rounds of interactive iteration, the agent swarm learns how to reduce global prediction errors through spatial collaboration even when local features are noisy (e.g., a grid sensor malfunctions). The trained weights are solidified into an inference model and deployed to an underground edge computing center to support real-time spatial risk source location decisions.

[0079] To ensure that the multi-agent reinforcement learning network has a complete and reproducible training and inference loop, this embodiment strictly defines the core interaction elements of each virtual monitoring agent.

[0080] First, in terms of state space definition, the local state vector acquired by each virtual monitoring agent specifically includes the following feature dimensions: local acoustic emission energy density, local mean of nonlinear modulation coefficient, gas concentration change rate, local mapping value of remaining activation energy within the spatial monitoring grid it is responsible for, and spatial relative distance from the center point of the spatial monitoring grid to the current real-time position of the cutting head.

[0081] Secondly, in terms of action space definition, the action set of each agent is strictly defined as a set of discrete risk assessment levels, specifically divided into five discrete actions from level one to level five, corresponding to local risk probability outputs of 0%, 25%, 50%, 75%, and 100%, respectively.

[0082] Secondly, regarding the reward allocation mechanism, the system utilizes the globally shared reward value (i.e., the negative rate of change of remaining activation energy) calculated in the preceding steps, combined with the local contribution weights output by the attention allocation layer, to decompose the reward. The actual local reward value obtained by each agent after each decision is equal to the product of the globally shared reward value and the corresponding local contribution weight of that agent, thereby solving the credit allocation problem in multi-agent collaboration.

[0083] Finally, regarding the closed-loop early warning triggering mechanism, when the local risk probability output by any spatial monitoring grid in the system exceeds a preset first danger probability threshold, and this high-risk state is maintained continuously for more than a preset judgment duration in the time dimension, the system will mark that grid as a candidate area for risk sources. When the total number of marked candidate areas globally reaches or exceeds a preset triggering threshold, the collaborative risk decision-making layer will immediately output the corresponding coal and gas outburst early warning level and the precise coordinates of the risk source.

[0084] By introducing a multi-agent reinforcement learning network, distributed intelligent monitoring and precise credit allocation for dynamic tunneling areas were achieved. Utilizing spatial gridded agent mapping, the macroscopic critical early warning task was decomposed into microscopic spatiotemporal feature identification; combining reward feedback based on the rate of change of remaining activation energy with attention weight decomposition, key spatial nodes leading to system instability were effectively identified. Example

[0085] As the tunneling machine advances along the tunnel axis into the structural zone, the system retrieves the output current and frequency mapping address of the cutting motor in real time via an industrial Ethernet communication interface. A displacement monitoring device records the reciprocating cutting coordinates of the cutting head within the tunnel cross-section. Simultaneously, piezoelectric acoustic emission sensors deployed in the forward borehole capture stress release signals generated by the cutting disturbance in the coal and rock mass. All data streams are globally timestamped via a PTP clock server, and microsecond-level sampling step alignment is performed in the industrial control computer's memory, forming a foundational dataset reflecting the dynamic interaction between the machine and the rock.

[0086] When the cutting head contacts the hard coal body, the high-frequency acoustic emission signal is synchronized and framed using the cutting head's rotation frequency. In the frequency domain, the system locks the fundamental frequency energy centered on the rotation frequency and uses narrowband filtering to extract the symmetrical sideband components generated by the "breathing effect" of micro-fractures within the coal and rock. The calculated nonlinear modulation coefficient sequence enters the sliding window analysis stage. When approaching the tectonic stress concentration zone, the system identifies that the first-order hysteresis autocorrelation sequence of this coefficient exhibits a significant linear growth trend, and the variance value exceeds the rolling dynamic threshold set based on the previous five minutes of historical background.

[0087] The nonlinear modulation coefficients and the gas concentration sequence in the return airflow were symbolized. By calculating the bidirectional transfer entropy, it was found that the current forward transfer entropy is significantly greater than the backward transfer entropy, and the net transfer entropy value is positive. Based on this, the system determines that the current instability dynamics are mainly driven by the concentration of ground stress (stress-driven risk). At the same time, the short-time Fourier transform monitoring showed that the spectral centroid of the acoustic emission signal dropped sharply from the high-frequency steady-state region to the low-frequency band. The frequency domain span characteristics reflect that the internal fractures of the coal body have changed from microscopic development to macroscopic interconnection.

[0088] The aforementioned heterogeneous feature vectors are input into a pre-trained rheological instability model with pre-defined model parameters. The feature space alignment layer maps power, triggering operators, causal weights, and frequency domain span features to a high-dimensional hidden space. The time-varying rheological evolution layer, combined with the current truncation perturbation intensity, simulates the decay creep process of the coal-rock mass under stress concentration conditions. The residual energy decoding layer calculates through nonlinear mapping that the remaining activation energy of the current coal-rock system, which is now below the preset safety threshold, indicates that the system has entered a physically critical state.

[0089] The multi-agent reinforcement learning network divides the area within ten meters in front of the cutting head into a voxel-based spatial monitoring grid based on the cutting head's position. Each virtual monitoring agent performs policy evaluation based on local monitoring features. The physical reward calculation layer provides strong negative feedback based on the sharp decrease in remaining activation energy. The attention allocation layer identifies the agent located in the top corner of the left side of the roadway as having the highest contribution through weight allocation. The collaborative risk decision layer ultimately outputs a "Level 1 Red Alert" and locates the center of the spatial monitoring grid corresponding to the risk source approximately 1.5 meters inside the coal seam on the outer side of the left side of the roadway. Based on this coordinate, targeted directional drilling for pressure relief was immediately performed on-site, successfully eliminating the potential dynamic disaster risk.

[0090] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A coal and gas outburst early warning method based on multi-dimensional time-frequency feature fusion, characterized in that, include: Acquire the cutting head rotation frequency, the instantaneous power sequence of the cutting motor, the cutting head advance coordinates, the broadband acoustic emission signal, and the gas concentration sequence; Using the rotation frequency of the cutting head as the fundamental frequency, the ratio of the sideband modulation component to the fundamental frequency energy is extracted from the broadband acoustic emission signal to obtain the nonlinear modulation coefficient sequence; Sliding window sampling is performed on the nonlinear modulation coefficient sequence, the autocorrelation coefficient and variance evolution value of adjacent window sequences are calculated, the monotonically increasing trend of the autocorrelation coefficient and the surge point of the variance evolution value are identified, and a critical instability physical triggering operator is generated. The transfer entropy algorithm is used to perform directed correlation analysis on the nonlinear modulation coefficient sequence and the gas concentration sequence, and outputs the causal correlation weight vector; the instantaneous frequency domain span feature of the spectral centroid trajectory of the broadband acoustic emission signal drops sharply to the low frequency band is extracted; By inputting the instantaneous power sequence, critical instability physical triggering operator, causal correlation weight vector, and instantaneous frequency domain span characteristics into the rheological instability model, the remaining activation energy of the coal-rock system evolving to the instability critical point is obtained. A multi-agent reinforcement learning network is constructed, using the decay rate of the remaining activation energy as the reward function, and combined with the cutting head advance coordinates to output the early warning level and location of the early warning area for coal and gas outbursts.

2. The method for early warning of coal and gas outbursts based on multi-dimensional time-frequency feature fusion according to claim 1, characterized in that, The process of acquiring the cutting head rotation frequency, the instantaneous power sequence of the cutting motor, the cutting head advance coordinates, the broadband acoustic emission signal, and the gas concentration sequence includes: real-time retrieval of the inverter's output parameters through the communication interface of the tunneling machine's electrical control system to acquire the cutting head rotation frequency and the instantaneous power sequence of the cutting motor; acquisition of the cutting head advance coordinates dynamically changing with the tunneling machine's advance through a displacement monitoring device installed on the tunneling machine body; acquisition of broadband acoustic emission signals through piezoelectric acoustic emission sensors fixed to the inner wall of the roadway or in the advance borehole, and analog-to-digital conversion processing through a high-speed data acquisition card; acquisition of the gas concentration sequence through a gas sensor installed in the return air area of ​​the working face; and connection of all acquired data to a unified clock source server, adding a global timestamp to each data stream to achieve time-series alignment acquisition of multi-source heterogeneous data.

3. The method for early warning of coal and gas outbursts based on multi-dimensional time-frequency feature fusion according to claim 1, characterized in that, The process of obtaining the nonlinear modulation coefficient sequence includes: synchronously framing the broadband acoustic emission signal using the rotation frequency of the cutting head, and extracting acoustic emission time-domain sample segments corresponding to the cutting period; performing spectral transformation processing on the acoustic emission time-domain sample segments to identify the fundamental frequency energy distribution with the rotation frequency of the cutting head as the center frequency; using narrowband bandpass filtering to obtain symmetrical sideband components distributed on both sides of the center frequency, wherein the symmetrical sideband components are generated by the nonlinear modulation effect of the internal fissures of the coal and rock mass on the cutting vibration; integrating the energy at the center frequency to obtain the total fundamental frequency energy, and simultaneously integrating the energy at the symmetrical sideband components to obtain the total sideband modulation energy; calculating the ratio of the total sideband modulation energy to the total fundamental frequency energy, and mapping the calculation result with the time axis to generate the nonlinear modulation coefficient sequence.

4. The method for early warning of coal and gas outbursts based on multi-dimensional time-frequency feature fusion according to claim 1, characterized in that, The process of generating the critical instability physical trigger operator includes: setting a fixed window width and sliding step size, truncating the nonlinear modulation coefficient sequence to obtain overlapping adjacent window sub-sequences; calculating the autocorrelation coefficient of each sub-sequence under first-order lag, and calculating the variance of each sub-sequence relative to its own mean, constructing a sequence of autocorrelation coefficient evolution over time and a sequence of variance evolution, respectively; performing linear regression fitting on the sequence of autocorrelation coefficient evolution over time, and calculating the fitting slope; when the fitting slope is continuously positive within a preset number of consecutive windows, and the autocorrelation coefficient value of the current window enters a preset high-value threshold range, it is determined that there is a monotonically increasing trend; calculating the rolling mean and rolling standard deviation of the variance evolution sequence before the current window; when the variance value of the current window exceeds the sum of the rolling mean and the rolling standard deviation by a preset multiple, it is determined to be the surge point; when both the monotonically increasing trend and the surge point are identified simultaneously, the autocorrelation coefficient and variance value of the current window are multiplied to generate the critical instability physical trigger operator.

5. The method for early warning of coal and gas outbursts based on multi-dimensional time-frequency feature fusion according to claim 1, characterized in that, The process of outputting the causal correlation weight vector includes: performing symbolization processing on the nonlinear modulation coefficient sequence and the gas concentration sequence respectively, transforming the continuous numerical sequence into a discrete state sequence reflecting the fluctuation trend; calculating the forward transfer entropy from the nonlinear modulation coefficient state sequence to the gas concentration state sequence and the backward transfer entropy from the gas concentration state sequence to the nonlinear modulation coefficient state sequence based on a preset time lag order; calculating the difference between the forward transfer entropy and the backward transfer entropy to obtain a net transfer entropy value used to characterize the asymmetry of information transmission; identifying the current dynamic dominant factor according to the positive or negative polarity of the net transfer entropy value; if the net transfer entropy value is positive, it is determined to be a stress-driven risk; if the net transfer entropy value is negative, it is determined to be a gas-driven risk; performing normalization processing on the forward transfer entropy and the backward transfer entropy, and combining the normalized value with the polarity identifier of the dynamic dominant factor to generate a causal correlation weight vector.

6. The method for early warning of coal and gas outbursts based on multi-dimensional time-frequency feature fusion according to claim 1, characterized in that, The process of extracting the instantaneous frequency domain span feature of the spectral centroid trajectory of the broadband acoustic emission signal plunging to the low-frequency band includes: using short-time Fourier transform to convert the broadband acoustic emission signal into a two-dimensional power spectral density matrix that evolves over time; performing a weighted average calculation of the frequency components and corresponding energy amplitudes for each time slice in the two-dimensional power spectral density matrix to extract a spectral centroid trajectory sequence reflecting the drift characteristics of the energy concentration region; performing a first-order gradient operation on the spectral centroid trajectory sequence to identify the instantaneous moment when the negative first-order gradient exceeds a preset steep descent threshold, which is taken as the trigger time point of the low-frequency band steep descent; extracting the average spectral centroid within a preset time period before the trigger time point as the high-frequency steady-state value, and extracting the spectral centroid under the instantaneous slice after the trigger time point as the low-frequency transient value; calculating the difference between the high-frequency steady-state value and the low-frequency transient value to generate the instantaneous frequency domain span feature characterizing the abrupt change in the scale of coal body fracture.

7. The method for early warning of coal and gas outbursts based on multi-dimensional time-frequency feature fusion according to claim 1, characterized in that, The rheological instability model includes: Feature Space Alignment Layer: Receives instantaneous power sequence, critical instability physical triggering operator, causal correlation weight vector and instantaneous frequency domain span feature, and performs dimensionality upscaling on each heterogeneous feature, mapping it to a unified feature hidden space; Physical causal interaction layer: Using the causal correlation weight vector as the query matrix, attention weighting is performed on the instantaneous frequency domain span features, and the features are concatenated with the critical instability physical triggering operator to generate a physical feature vector representing the damage accumulation state; Time-varying rheological evolution layer: Using the physical feature vector and the instantaneous power sequence as time-series input, the strain rate evolution and damage accumulation process of coal and rock mass under truncation disturbance are simulated through a cyclic calculation unit containing coal and rock rheological constitutive constraints, and the time-varying rheological state tensor is output. Energy Residual Decoding Layer: Performs nonlinear dimensionality reduction mapping on the time-varying rheological state tensor, calculates the physical difference between the current energy accumulation state and the preset critical destruction energy threshold, and outputs the remaining activation energy of the coal-rock system evolving to the instability critical point.

8. The method for early warning of coal and gas outbursts based on multi-dimensional time-frequency feature fusion according to claim 1, characterized in that, The multi-agent reinforcement learning network includes: Spatial intelligent agent mapping layer: The current tunneling area is divided into multiple spatial monitoring grids according to the cutting head advance coordinates, and a virtual monitoring intelligent agent is configured for each spatial monitoring grid to perform risk assessment actions; Physical reward calculation layer: Obtain the time series of the remaining activation energy, calculate the negative rate of change of the remaining activation energy over time, and generate a globally shared reward value characterizing the severity of instability of the coal-rock system; Attention allocation layer: The attention mechanism is used to calculate the spatiotemporal correlation between the local features extracted by each virtual monitoring agent and the global shared reward value, and the global shared reward value is decomposed into the local contribution weights corresponding to each virtual monitoring agent; Collaborative Risk Decision Layer: Construct a multi-agent reinforcement learning network, using the decay rate of remaining activation energy as the reward function, and combine it with the cutting head advance coordinates to output the coal and gas outburst warning level and warning area location.

Images