Three-dimensional dynamic numerical simulation method of coal spontaneous combustion in goaf during working face advancing

Abstract

The application discloses a goaf coal spontaneous combustion three-dimensional dynamic numerical simulation method during working face advancing, relates to the technical field of coal spontaneous combustion three-dimensional dynamic numerical simulation, and comprises the following steps: based on the region of the coded identification, combining the three-dimensional disturbance offset track and the ventilation reverse distribution characteristics in the structure of the boundary disturbance, generating the boundary form uncertainty description body, and determining the boundary form uncertainty grade under the condition of local collapse of the goaf boundary space structure in the advancing process according to the space offset amplitude, disturbance duration and coupling strength in the boundary form uncertainty description body. The application introduces multi-source disturbance identification, boundary form uncertainty quantification, oxygen diffusion dynamic adjustment and oxygen consumption linkage feedback mechanism, realizes the synchronous update of the oxygen diffusion parameters and the coal oxidation reaction process under the boundary mutation condition, and thus improves the accuracy and reliability of the goaf coal spontaneous combustion three-dimensional dynamic numerical simulation.

Description

Technical Field

[0001] This invention relates to the field of three-dimensional dynamic numerical simulation technology for spontaneous combustion of coal, specifically to a method for three-dimensional dynamic numerical simulation of spontaneous combustion of coal in the goaf during the advancement of the working face. Background Technology

[0002] Three-dimensional dynamic numerical simulation of spontaneous combustion of coal in the goaf during the advance of the working face refers to the simulation of the spontaneous combustion of coal in the goaf in real time during coal mining, as the working face continues to advance and the goaf expands and evolves. This simulation involves establishing a three-dimensional mathematical model that includes the spatial geometry of the goaf, the oxidation reaction characteristics of the coal, oxygen diffusion, heat conduction and accumulation, and gas flow patterns, among other multi-physics coupling relationships. The simulation predicts the location, timing, and evolution trend of spontaneous combustion, thus achieving early warning and prevention. In existing technologies, this type of simulation is typically achieved through the following steps: first, a three-dimensional modeling software (such as Surfer, FLAC3D, etc.) is used to establish the... The process involves several steps: first, a geometric model of the goaf is created, and the distribution of coal and the boundary conditions of gas and air channels are determined based on geological survey data. Second, a kinetic model of the exothermic coal oxidation reaction is introduced to simulate the coupled evolution of oxygen concentration and temperature field. Then, the goaf model is dynamically updated in conjunction with the working face advance rate, and time is introduced as a coupling dimension to simulate changes over time. Next, numerical solutions are obtained using the finite element method or finite difference method, and the temperature, oxygen concentration distribution, and coal spontaneous combustion index changes at each time step are output. Finally, a 3D visualization and spontaneous combustion risk assessment are performed through a visualization platform to assist decision-making. The entire process typically includes key steps such as model initialization, parameter assignment, advance rule setting, physical field solution, result output, and visualization analysis.

[0003] The existing technology has the following shortcomings: During the three-dimensional dynamic numerical simulation of spontaneous combustion of coal in the goaf during face advancement, when advancing to areas with severe residual coal compaction, the local boundaries of these areas are prone to nonlinear collapse or spatial distortion due to the disturbance of the working face, forming boundary structures that are difficult to accurately reconstruct in real time during the modeling process. This boundary deformation causes the local space to change from an open state to a near-closed state, making it difficult for oxygen to actually enter the area. However, due to the suddenness and unpredictability of this type of boundary indentation, existing three-dimensional dynamic numerical simulation technology for spontaneous combustion of coal in the goaf during face advancement cannot dynamically adjust the oxygen diffusion conditions of the simulated boundary based on the uncertainty of the boundary morphology when the spatial structure of the goaf boundary indentation occurs during the advancement process. It still relies on the static boundary settings in the original geometric model to continue to perform gas diffusion and oxygen transport calculations. This results in an erroneous overestimation of the oxygen supply at the boundary in the simulation, leading to a serious misjudgment of the local oxygen concentration distribution and coal oxidation reaction rate. Furthermore, it causes the simulation system to misidentify areas with actual heat accumulation risks as well-ventilated low-risk areas, delaying the identification and early warning of spontaneous combustion risks, and thus reducing the effectiveness and reliability of the entire numerical simulation results in guiding actual prevention and control.

[0004] The information disclosed in the background section is only intended to enhance the understanding of the background of this disclosure, and therefore may include information that does not constitute prior art known to those skilled in the art. Summary of the Invention

[0005] The purpose of this invention is to provide a three-dimensional dynamic numerical simulation method for spontaneous combustion of coal in the goaf during the advancement of the working face, so as to solve the problems in the background art mentioned above.

[0006] To achieve the above objectives, the present invention provides the following technical solution: a three-dimensional dynamic numerical simulation method for spontaneous combustion of coal in the goaf during working face advancement, specifically including the following steps: S1. By collecting the micro-seismic response value, ventilation reversal rate and stress release fluctuation value of compacted residual coal in the propulsion unit, a boundary disturbance identification structure is constructed. Based on the abnormal co-occurrence threshold in the boundary disturbance identification structure, it is determined whether the local subsidence of the goaf boundary space structure occurs during the propulsion process, and the area where the goaf boundary space structure has local subsidence is coded and marked. S2. Based on the coded and identified areas, and combined with the three-dimensional disturbance offset trajectory and ventilation reverse distribution characteristics within the boundary disturbance identification structure, a boundary morphology uncertainty description body is generated. Based on the spatial offset amplitude, disturbance duration, and coupling strength in the boundary morphology uncertainty description body, the boundary morphology uncertainty level is determined when a local subsidence of the goaf boundary spatial structure occurs during the advancement process. S3. Establish an oxygen diffusion control trigger set based on the uncertainty level of the boundary morphology, and associate the boundary disturbance state in the oxygen diffusion control trigger set with the diffusion parameters to realize the dynamic adjustment of the simulated boundary oxygen diffusion conditions. S4. The dynamically adjusted simulated boundary oxygen diffusion conditions are introduced into the coal oxidation reaction calculation process. The oxygen diffusion changes and oxidation reaction parameters are linked by constructing an oxygen consumption linkage factor set. S5. Based on the changes in the oxygen consumption linkage factor set, a simulated control equilibrium sequence is constructed, and the calculation frequency, spatial resolution, and calculation step size of the propulsion unit are dynamically controlled according to the response change rate in the simulated control equilibrium sequence.

[0007] Preferably, S1 specifically includes the following steps: S101. Collect microseismic response values, ventilation reversal rate, and compacted residual coal stress release fluctuation values ​​within the propulsion unit, and establish microseismic response value sequences, ventilation reversal rate sequences, and compacted residual coal stress release fluctuation value sequences respectively. Synchronously organize the three types of sequences according to the spatial location of the propulsion unit to form the basic dataset for boundary disturbance identification structure. S102. The three types of sequences in the basic dataset are jointly arranged according to the advancement time sequence to form a boundary disturbance identification structure. The boundary disturbance identification structure uses the time overlap relationship between the sudden increase point of microseismic response value, the change point of ventilation reverse rate and the abnormal point of stress release fluctuation value of compacted residual coal as an abnormal co-occurrence reference index. By comparing the abnormal co-occurrence threshold, it is determined whether the local inward subsidence of the boundary space structure of the goaf occurs during the advancement process. S103. After determining that a local subsidence of the boundary space structure of the goaf has occurred, the main identifier is the advancing unit number, and the auxiliary identifier is the corresponding abnormal co-occurrence location in the boundary disturbance identification structure. The main identifier and the auxiliary identifier are combined to form a coded identifier, and the coded identifier is recorded as the area marking result of the local subsidence of the boundary space structure of the goaf during the advancing process.

[0008] Preferably, S102 is as follows: Extract the timestamp information of the microseismic response value sequence, ventilation reverse rate sequence and compacted residual coal stress release fluctuation value sequence from the basic dataset, establish a unified time axis according to the advancement sequence, and arrange all observation point data of the three types of sequences onto the unified time axis. By screening out abrupt increases in the microseismic response value sequence, changes in the ventilation reverse rate sequence, and anomalies in the stress release fluctuation value sequence of compacted residual coal, the overlap of the three types of characteristic points on a unified time axis is used as an abnormal co-occurrence reference index for boundary disturbance identification structures according to the spatial coordinates of each propulsion unit. Based on a pre-set abnormal co-occurrence threshold, when the three types of feature points in the propulsion unit coincide on a unified time axis, it is determined that the propulsion unit has a local indentation of the boundary space structure of the goaf during the propulsion process, and the judgment result is written into the indentation marker of the boundary disturbance identification structure.

[0009] Preferably, S2 specifically includes the following steps: S201. Extract the three-dimensional disturbance offset trajectory within the area corresponding to the coded identifier, construct a disturbance trajectory vector group through the three-axis displacement sequence, and perform time-by-time offset analysis on each vector group to calculate the spatial offset amplitude and disturbance duration; simultaneously extract ventilation reverse distribution features, spatially match the disturbance trajectory with the ventilation reverse section, and construct a fusion structure of three-dimensional disturbance offset trajectory and ventilation reverse distribution features. S202. Based on the disturbance amplitude, disturbance duration and ventilation reversal degree of each spatial point in the fusion structure, calculate the corresponding coupling strength according to the angle between the spatial overlap ratio and the changing trend. Combine the maximum displacement amplitude and minimum duration of the three-dimensional disturbance offset trajectory to generate a boundary morphological uncertainty description body with structural integrity and quantitative characteristics. S203. Based on the spatial offset amplitude, disturbance duration and coupling strength extracted from the boundary morphology uncertainty description, a scoring model is constructed. Feature weights are assigned to each and weighted calculations are performed to obtain the comprehensive score value of boundary morphology uncertainty. The comprehensive score value is compared with the pre-set uncertainty level range and classified into high-level, medium-level or low-level uncertainty level results, and bound to the coded identification area.

[0010] Preferably, S201 specifically refers to: Based on the spatial region corresponding to the coded identifier, the sequence of displacement observation points arranged along the three axes is extracted, and integrated into a three-axis displacement sequence in chronological order. Based on the displacement changes of each axis in adjacent time periods, a disturbance trajectory vector group is constructed. The three-dimensional displacement changes are calculated vector by vector to obtain the spatial offset amplitude of the disturbance trajectory. The time periods in which the disturbance exists continuously are recorded to obtain the disturbance duration. Based on the spatial coordinates of the disturbance trajectory, the ventilation reversal rate data corresponding to the disturbance point in the ventilation monitoring network is called to extract the segment boundary of the ventilation direction change, and the spatial path of the disturbance trajectory is matched with the ventilation reversal segment to identify the overlap range and matching degree between the trajectory point and the ventilation reversal segment. The matching results are normalized, and the processed matching degree is jointly mapped with the three-axis perturbation trajectory. Based on the spatial overlap relationship between the trajectory spatial path and the ventilation reverse distribution, a fusion structure of the three-dimensional perturbation offset trajectory and ventilation reverse distribution features is constructed for subsequent generation of boundary morphology uncertainty descriptor.

[0011] Preferably, S202 specifically refers to: Traverse each spatial point in the fusion structure of the three-dimensional disturbance offset trajectory and ventilation reverse distribution characteristics, extract the disturbance amplitude, disturbance duration and ventilation reverse degree data of each spatial point, construct a trend vector based on the disturbance direction and the ventilation reverse change direction, calculate the spatial change trend angle between the disturbance direction and the ventilation reverse change direction, and determine the spatial overlap ratio based on the overlap ratio of the angle and the spatial point in the ventilation change section. Based on the spatial overlap ratio and the angle of change trend of each spatial point, a spatial coupling strength calculation function is constructed according to the product of cosine similarity and overlap, outputting the coupling strength value of the corresponding spatial point, and summarizing the coupling strength set of all spatial points to quantitatively characterize the reverse coupling characteristics of disturbance and ventilation. The range of global maximum and mean values ​​is extracted from the coupling strength set, and combined with the maximum offset amplitude and minimum disturbance duration in the three-dimensional disturbance offset trajectory, and uniformly mapped to the boundary disturbance spatial domain. The boundary morphology uncertainty description body with closed structure and complete numerical data is output in the form of a spatial grid.

[0012] Preferably, S3 is as follows: Based on the uncertainty level of the boundary morphology, an oxygen diffusion regulation trigger set is constructed. The boundary disturbance states are set as low disturbance state, medium disturbance state and high disturbance state respectively. Predefined diffusion regulation parameters are assigned to each boundary disturbance state, including the adjustment amplitude of oxygen concentration diffusion coefficient and the correction value of local airflow velocity. A one-to-one correspondence is established between the boundary perturbation states of each oxygen diffusion regulation trigger and the corresponding diffusion regulation parameters. A parameter lookup table is established using the propulsion unit code as an index. The orderly association between the boundary perturbation states and the diffusion parameters is achieved through this parameter lookup table. Based on the parameter retrieval table, the simulated boundary configuration is dynamically interpolated and updated. When the propulsion unit enters the simulation iteration stage, the boundary disturbance state and diffusion parameters corresponding to the current code are automatically extracted, and the simulated boundary oxygen diffusion conditions are updated in real time to keep them synchronized with the disturbance state.

[0013] Preferably, S4 is as follows: The oxygen concentration diffusion coefficient and local airflow velocity parameters contained in the dynamically adjusted simulated boundary oxygen diffusion conditions are extracted. An oxygen diffusion input sequence is established according to the coding order of the propulsion unit and the simulation time step. The oxygen diffusion input sequence is then spatially mapped to the boundary nodes in the coal oxidation reaction calculation process. In the calculation of coal oxidation reaction, the oxygen concentration diffusion coefficient and local gas flow velocity parameters corresponding to each time step in the oxygen diffusion input sequence are input into the reaction region model. The rate of change of oxygen content, the rate of change of oxygen partial pressure gradient and the time derivative of oxygen concentration in a unit volume region are calculated to form a set of basic reaction input factors characterizing the changes in oxygen diffusion. Based on the basic reaction input factor set, an oxygen consumption linkage factor set is constructed. The oxygen content change rate, oxygen partial pressure gradient change rate, and oxygen concentration time derivative are dynamically mapped to the reaction rate, reaction order, and oxygen consumption function parameters in the coal oxidation reaction model, respectively. Through the oxygen consumption linkage factor set, the linkage processing of oxygen diffusion changes and oxidation reaction parameters is realized, and the oxygen consumption behavior in the reaction process is updated synchronously.

[0014] Preferably, S5 is as follows: Based on the changes in the oxygen consumption linkage factors, the numerical differences of the oxygen content change rate, oxygen partial pressure gradient change rate, and oxygen concentration time derivative in continuous time steps are extracted. The change rate of each factor is calculated, and a response change rate sequence is constructed in units of propulsion units. The response change rate sequence in each time step is normalized to form a simulated control equilibrium sequence for dynamic control analysis. Based on the response change rate in the simulated control equilibrium sequence, a control mapping relationship is established between the response change rate and the calculation frequency, spatial resolution and calculation step size of the propulsion unit. The propulsion units are classified according to the range of response change rate, and corresponding calculation frequency adjustment value, spatial resolution adjustment value and calculation step size adjustment value are assigned to different levels of propulsion units, thus forming a set of control parameters that match the simulated control equilibrium sequence. The set of control parameters is applied to the simulation system. Using the propulsion unit as an index, the corresponding control parameters are automatically selected based on the rate of change of the response in the simulation control equilibrium sequence corresponding to the current time step. The calculation frequency, spatial resolution and calculation step size of the propulsion unit are dynamically controlled, and adaptive allocation of computing resources and synchronous adjustment of simulation accuracy are realized during the simulation iteration process.

[0015] The technical effects and advantages provided by the present invention in the above technical solution are as follows: 1. This invention constructs a boundary disturbance identification structure based on the fusion of multi-source data from microseismic response, ventilation reversal, and stress fluctuations. This structure can accurately identify the local indentation state of the boundary spatial structure of the goaf during the advancement process. Through spatial coupling analysis of three-dimensional disturbance offset trajectories and ventilation reversal characteristics, a boundary morphology uncertainty description with spatial integrity and quantitative characteristics is generated. The uncertainty level is then classified by a comprehensive score, thereby dynamically guiding the real-time update of oxygen diffusion control conditions. This responsive modeling mechanism for changes in the gas transport environment caused by abrupt boundary morphology changes overcomes the problem of oxygen supply misjudgment caused by static boundary settings in existing methods, improving the model's accuracy in identifying spontaneous combustion risk areas.

[0016] 2. This invention constructs a set of oxygen consumption linkage factors, dynamically mapping boundary adjustment parameters such as oxygen concentration diffusion coefficient and local airflow velocity to key reaction parameters in the coal oxidation reaction model. This enables direct driving and feedback correction of reaction behavior by changes in oxygen diffusion. Furthermore, a simulation-controlled equilibrium sequence is introduced to dynamically adjust the computational frequency, spatial resolution, and calculation step size of the propulsion unit according to the rate of change in oxygen consumption behavior. This achieves adaptive allocation of simulation resources and high-precision analysis of key areas. The overall scheme constructs a self-sensing, self-regulating, and self-converging coupled simulation system. While ensuring computational efficiency, it enhances the adaptability and responsiveness of the simulation scenario to complex real-world disturbances, providing strong technical support for the dynamic monitoring and intelligent prevention of coal spontaneous combustion risks in goaf areas. Attached Figure Description

[0017] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0018] Figure 1 This is a flowchart illustrating the three-dimensional dynamic numerical simulation method for spontaneous combustion of coal in the goaf during the working face advancement of this invention. Detailed Implementation

[0019] Exemplary embodiments will now be described more fully with reference to the accompanying drawings. However, these exemplary embodiments can be implemented in many forms and should not be construed as limited to the examples set forth herein; rather, they are provided so that the description of this disclosure will be more complete and fully convey the concept of the exemplary embodiments to those skilled in the art.

[0020] This invention provides, for example Figure 1 The three-dimensional dynamic numerical simulation method for spontaneous combustion of coal in the goaf during the working face advancement shown herein includes the following steps: S1. By collecting the micro-seismic response value, ventilation reversal rate and stress release fluctuation value of compacted residual coal in the propulsion unit, a boundary disturbance identification structure is constructed. Based on the abnormal co-occurrence threshold in the boundary disturbance identification structure, it is determined whether the local subsidence of the goaf boundary space structure occurs during the propulsion process, and the area where the goaf boundary space structure has local subsidence is coded and marked. In this embodiment, S1 specifically includes the following steps: S101. Collect microseismic response values, ventilation reversal rate, and compacted residual coal stress release fluctuation values ​​within the propulsion unit, and establish microseismic response value sequences, ventilation reversal rate sequences, and compacted residual coal stress release fluctuation value sequences respectively. Synchronously organize the three types of sequences according to the spatial location of the propulsion unit to form the basic dataset for boundary disturbance identification structure. In the three-dimensional dynamic numerical simulation of spontaneous combustion of coal in the goaf, to achieve dynamic identification of local boundary disturbances during the advancement process, it is necessary to simultaneously collect three highly sensitive indicators around each advancement unit: microseismic response value, ventilation reversal rate, and compacted residual coal stress release fluctuation value. Microseismic response values ​​can be obtained through seismic detectors deployed around the goaf to monitor the intensity of high-frequency microseismic activity; the ventilation reversal rate can be measured by underground wind speed and direction sensors to measure the directional change of airflow at different locations; and the compacted residual coal stress release fluctuation value is monitored by deploying a stress sensor array to track the unstable fluctuation trend of stress release in local residual coal under disturbance. After collection, the three types of data need to be organized according to the one-to-one correspondence between collection time and spatial coordinates, and divided into three sequences based on the advancement unit as the smallest spatial block, forming microseismic response value sequences, ventilation reversal rate sequences, and compacted residual coal stress release fluctuation value sequences. Subsequently, based on the spatial positioning of the advancement unit, the three sequences are normalized and aligned along the time axis to complete the synchronous organization of the three sequences. The resulting 3D dataset will serve as the foundation for identifying boundary disturbances, ensuring that subsequent identification and analysis can be conducted within a unified temporal and spatial framework, thus avoiding misjudgments of the spatiotemporal correspondence of local disturbances. For example, if a propulsion unit experiences a sudden increase in microseismic values ​​during a certain period, accompanied by airflow reversal and stress release fluctuations occurring simultaneously, forming a co-occurrence of three types of sequence anomalies, the system can infer that the unit has experienced a potential indentation of its boundary structure.

[0021] An advancing unit refers to the smallest computational space unit divided during the 3D modeling of the goaf to adapt to the advancing rhythm of the working face. It typically presents a rectangular or cubic structure. The microseismic response value reflects the intensity of elastic waves generated by the coal body or surrounding rock under microscale stress disturbance. It can be collected by multi-point geophones deployed around the mining area, and its unit is energy level or amplitude value; its sequence reflects the trend of microseismic activity changes of the unit at different time points. The ventilation reversal rate is a ratiomatic index used to measure whether the local airflow direction deviates from the preset airflow direction. It is generated by detecting the actual wind direction using an array of wind speed and direction sensors and comparing the offset angle with the desired wind direction. Its sequence reflects the continuity and repetition of airflow disturbance at different times. The stress release fluctuation value of compacted residual coal measures the stress release change of the residual coal body under dynamic compaction and disturbance. It is usually obtained by stress gauges buried in the residual coal body or its adjacent surrounding rock. Its sequence is constructed by analyzing the peak value variation amplitude and periodic abrupt change points in the stress release curve. When constructing sequences for the three types of data, spatial classification must be performed based on the unique identifier of the propulsion unit. Combined with the synchronization standard of the time axis in the propulsion rhythm, different types of data should be precisely paired to ensure that the three types of response variables can express boundary perturbation phenomena in the same time frame.

[0022] S102. The three types of sequences in the basic dataset are jointly arranged according to the advancement time sequence to form a boundary disturbance identification structure. The boundary disturbance identification structure uses the time overlap relationship between the sudden increase point of microseismic response value, the change point of ventilation reverse rate and the abnormal point of stress release fluctuation value of compacted residual coal as an abnormal co-occurrence reference index. By comparing the abnormal co-occurrence threshold, it is determined whether the local inward subsidence of the boundary space structure of the goaf occurs during the advancement process. S103. After determining that a local subsidence of the boundary space structure of the goaf has occurred, the main identifier is the advancing unit number, and the auxiliary identifier is the corresponding abnormal co-occurrence location in the boundary disturbance identification structure. The main identifier and the auxiliary identifier are combined to form a coded identifier, and the coded identifier is recorded as the area marking result of the local subsidence of the boundary space structure of the goaf during the advancing process.

[0023] During the advancement process, to ensure accurate tracking and response to areas of localized indentation in the boundary space structure during subsequent simulations, the corresponding areas need to be clearly identified after localized indentation identification. First, the advancement unit number is used as the primary identifier. This number is typically assigned based on the segmented order of the working face along the advancement direction, possessing spatial uniqueness and reflecting the positional relationship of each advancement unit within the overall advancement path. Second, the abnormal co-occurrence location detected in the boundary disturbance identification structure is used as a secondary identifier. This location generally exists in the form of a timestamp or sampling point index, precisely indicating the time node and corresponding data row where the localized indentation occurred. Combining the primary and secondary identifiers forms a two-dimensional code, such as "Advancement Unit Number_Abnormal Time Number" or "Advancement Segment ID_Feature Point Index." This code not only reflects the location of the indentation but also encompasses temporal information, giving the identification result both spatial and temporal characteristics. This encoding method can be implemented through string concatenation, bitwise operations, or matrix mapping, offering advantages such as high computational efficiency, clear storage structure, and convenient indexing operations. After encoding, the result is recorded as a regional marker of local subsidence in the boundary spatial structure of the goaf, and embedded as a core label in the boundary disturbance identification structure to drive subsequent boundary morphology uncertainty assessment, oxygen diffusion regulation logic invocation, and risk level classification. By constructing a unique and resolvable coded identifier, high-precision labeling and source tracing of local boundary anomaly areas can be achieved, ensuring that the simulation system has stable and coherent information transmission and response control capabilities in dynamic environments.

[0024] In this embodiment, S102 specifically refers to: Extract the timestamp information of the microseismic response value sequence, ventilation reverse rate sequence and compacted residual coal stress release fluctuation value sequence from the basic dataset, establish a unified time axis according to the advancement sequence, and arrange all observation point data of the three types of sequences onto the unified time axis. When constructing the boundary disturbance identification structure, to ensure the temporal consistency of the microseismic response value sequence, the ventilation reversal rate sequence, and the compacted residual coal stress release fluctuation value sequence, it is necessary to first extract the timestamp information corresponding to each observation point in the three types of sequences. The timestamp information is usually in milliseconds or seconds, representing the actual acquisition time of each data item, and is uniformly generated by the clock system built into the data acquisition equipment. After acquisition, all propulsion units are numbered and sorted according to the propulsion sequence, and the timestamp information is added to the propulsion sequence database. Next, based on the propulsion rhythm and the progress of each unit, a unified time axis is constructed as the data alignment benchmark for the three types of sequences. In specific implementation, linear interpolation can be used to merge the three types of data at different sampling frequencies to standard sampling points on a unified time axis. For example, each minute can be a time node, and the values ​​of each type of sequence at that time node can be merged or interpolated. For the microseismic response value, its maximum or average energy value within each minute can be extracted; for the ventilation reversal rate, the frequency of wind direction changes within that time period can be statistically analyzed; for the compacted residual coal stress release fluctuation value, the fluctuation amplitude of stress changes can be recorded. This alignment operation ensures that all data are synchronized with a consistent time resolution, avoiding misalignment of different indicators due to differences in data acquisition frequency. Taking a propulsion unit as an example, if all three types of data show abnormal values ​​within a certain minute on a unified time axis, it can be determined that a boundary disturbance event may have occurred in the corresponding region at that time point, thus providing a time-series basis for subsequent anomaly co-occurrence analysis. The core of this operation lies in building a unified and continuous time analysis platform, making the subsequent coupled analysis between the three types of data comparable and logically rigorous.

[0025] By screening out abrupt increases in the microseismic response value sequence, changes in the ventilation reverse rate sequence, and anomalies in the stress release fluctuation value sequence of compacted residual coal, the overlap of the three types of characteristic points on a unified time axis is used as an abnormal co-occurrence reference index for boundary disturbance identification structures according to the spatial coordinates of each propulsion unit. In constructing the boundary disturbance identification structure, to identify the nonlinear variation characteristics of the local boundary space structure in the goaf, it is necessary to extract representative abrupt change feature points from the microseismic response value sequence, ventilation reversal rate sequence, and compacted residual coal stress release fluctuation value sequence. Sudden increase points in the microseismic response value sequence refer to observation points where the microseismic energy level amplifies abnormally within a short period. These are typically determined by calculating the slope of the amplitude change rate over a continuous time period and combining it with an energy threshold to judge the nature of the sudden increase. These points often indicate rapid displacement release behavior of the coal body structure caused by disturbance. Change points in the ventilation reversal rate sequence refer to observation points where the airflow direction or intensity abruptly changes from positive flow to reverse or turbulent state within a short period. These can be extracted by the wind speed vector angle change rate and identified by combining it with the number of wind direction reversals. Abnormal points in the compacted residual coal stress release fluctuation value sequence refer to time points where nonlinear and severe fluctuations occur in the stress monitoring results. These points are generally identified by the abnormal amplitude of the statistical fluctuation variance within a continuous time window, representing areas where the coal body compaction structure has loosened, collapsed, or experienced sudden stress release. After locating the three types of feature points, these feature points are mapped onto the same propulsion unit based on its spatial coordinates, and each point is mapped to a specific time point on a unified time axis. The overlap is then calculated by statistically analyzing the ratio of the number of the three types of feature points occurring at the same time point. A higher overlap indicates that the propulsion unit is more likely to have experienced a boundary disturbance event at that time point. For example, if a propulsion unit simultaneously detects a sudden increase in microseismic values, a reversal of airflow direction, and severe stress fluctuations within a minute, the overlap is close to 100%, and this time point is identified as a high-confidence anomalous co-occurrence event point. This method enables unified spatiotemporal identification of anomalous events across multiple physical dimensions and provides a quantitative reference for determining whether local subsidence of the boundary spatial structure of the goaf exists during the propulsion process.

[0026] Based on a pre-set abnormal co-occurrence threshold, when the three types of feature points in the propulsion unit coincide on a unified time axis, it is determined that the propulsion unit has a local indentation of the boundary space structure of the goaf during the propulsion process, and the judgment result is written into the indentation marker of the boundary disturbance identification structure.

[0027] In the process of identifying disturbances at the boundary of goaf areas, in order to effectively determine whether there is local subsidence of the boundary spatial structure during the advancement process, it is necessary to construct a judgment standard for quantifying the co-occurrence degree of multiple types of feature points, namely, a pre-set abnormal co-occurrence threshold. This threshold represents the minimum requirement for the synchronous occurrence of microseismic response value surge points, ventilation reversal rate change points, and stress release fluctuation value anomalies on a unified time axis. Generally, this can be achieved by extracting the overlap ratio of the three types of feature points before the actual collapse by backtracking analysis of historical collapse cases, and then setting it through statistical modeling in conjunction with the boundary disturbance sensitivity under different working conditions. For example, in a typical collapse area, if the overlap rate of microseismic surge points, ventilation reversal rate change points, and stress fluctuation anomalies within a certain time window reaches more than 60%, then 60% can be used as a threshold reference. In the actual judgment process, after aligning all feature points in each advancement unit according to the time axis, the occurrence of the three types of feature points at each time point is checked sequentially, and the frequency of two or more of the three types of points occurring simultaneously at the same time point is counted. If the overlap at a certain point in time reaches or exceeds a pre-set abnormal co-occurrence threshold, it is determined that the propulsion unit has a high probability of boundary disturbance at that point in time, and it is further inferred that a local indentation of the boundary spatial structure may have occurred in that region. Once this judgment is confirmed, the point in time and the corresponding spatial unit location are recorded as an anomaly in the boundary disturbance identification structure, and a dedicated indentation marker is created in the structural data to indicate that the unit has experienced a structural anomaly, providing a precise reference location and time basis for subsequent boundary morphology uncertainty assessment and oxygen diffusion adjustment. By introducing a judgment mechanism with a pre-set threshold, the stability and accuracy of the identification can be significantly improved, avoiding misjudgments and omissions, and it is especially suitable for the intelligent identification of sudden, nonlinear, and unpredictable boundary deformation events.

[0028] S2. Based on the coded and identified areas, and combined with the three-dimensional disturbance offset trajectory and ventilation reverse distribution characteristics within the boundary disturbance identification structure, a boundary morphology uncertainty description body is generated. Based on the spatial offset amplitude, disturbance duration, and coupling strength in the boundary morphology uncertainty description body, the boundary morphology uncertainty level is determined when a local subsidence of the goaf boundary spatial structure occurs during the advancement process. In this embodiment, S2 specifically includes the following steps: S201. Extract the three-dimensional disturbance offset trajectory within the area corresponding to the coded identifier, construct a disturbance trajectory vector group through the three-axis displacement sequence, and perform time-by-time offset analysis on each vector group to calculate the spatial offset amplitude and disturbance duration; simultaneously extract ventilation reverse distribution features, spatially match the disturbance trajectory with the ventilation reverse section, and construct a fusion structure of three-dimensional disturbance offset trajectory and ventilation reverse distribution features. S202. Based on the disturbance amplitude, disturbance duration and ventilation reversal degree of each spatial point in the fusion structure, calculate the corresponding coupling strength according to the angle between the spatial overlap ratio and the changing trend. Combine the maximum displacement amplitude and minimum duration of the three-dimensional disturbance offset trajectory to generate a boundary morphological uncertainty description body with structural integrity and quantitative characteristics. S203. Based on the spatial offset amplitude, disturbance duration and coupling strength extracted from the boundary morphology uncertainty description, a scoring model is constructed. Feature weights are assigned to each and weighted calculations are performed to obtain the comprehensive score value of boundary morphology uncertainty. The comprehensive score value is compared with the pre-set uncertainty level range and classified into high-level, medium-level or low-level uncertainty level results, and bound to the coded identification area.

[0029] To quantitatively classify boundary morphological uncertainty, three key features are first extracted from the boundary morphological uncertainty descriptor: spatial offset amplitude, disturbance duration, and coupling strength. These three features represent the spatial severity of boundary deformation, temporal duration, and the synergistic relationship between disturbance and ventilation effects, respectively. A scoring model is constructed, assigning differentiated feature weights to each of the three feature values. For example, in areas prone to heat accumulation, coupling strength may be more critical and therefore assigned a higher weight. After feature normalization, the three feature values ​​are weighted to obtain a comprehensive score for boundary morphological uncertainty. This comprehensive score is then compared with a predefined level interval, and the uncertainty level of the area is determined based on the numerical value. For instance, if a coded area has a spatial offset amplitude of 0.9 meters, a disturbance duration of 8 minutes, and a coupling strength of 0.75, the calculated score is 0.82. The system maps this score to a higher level within the preset interval, ultimately binding the higher-level result to the coded area for accurate identification and classification.

[0030] Spatial offset amplitude refers to the maximum cumulative displacement of boundary disturbances in three-dimensional space, usually obtained by synthesizing the incremental values ​​of trajectory vectors along each axis; disturbance duration represents the stable existence time of a disturbance event within a continuous time period, used to determine the dynamic extensibility of the disturbance; coupling strength comprehensively reflects the degree of spatial coordination between the disturbance direction and the ventilation reverse direction. The scoring model is the core structure that integrates these three parameters. The weighted calculation involves multiplying each normalized indicator by a preset weight according to its risk contribution ratio and then summing them to generate a unified comprehensive score. The pre-defined uncertainty level range is generally determined by combining historical data and expert experience. For example, a score of 0.0-0.3 can be classified as low level, 0.3-0.6 as medium level, and 0.6-1.0 as high level to ensure the scientific nature and adaptability of the classification. The binding operation directly maps the level results to the coded identification area, so that each potential anomaly area has a queryable and traceable level attribute for adaptive parameter updates in subsequent dynamic simulations.

[0031] In this embodiment, S201 specifically refers to: Based on the spatial region corresponding to the coded identifier, the sequence of displacement observation points arranged along the three axes is extracted, and integrated into a three-axis displacement sequence in chronological order. Based on the displacement changes of each axis in adjacent time periods, a disturbance trajectory vector group is constructed. The three-dimensional displacement changes are calculated vector by vector to obtain the spatial offset amplitude of the disturbance trajectory. The time periods in which the disturbance exists continuously are recorded to obtain the disturbance duration. This process primarily revolves around constructing a set of disturbance trajectory vectors. Its core lies in using the spatial region corresponding to the coded identifier as the target region, and first extracting historical data sequences from pre-defined displacement observation points along the X, Y, and Z axes within that region. Each observation point records its spatial displacement value at various time points, which can be acquired using 3D laser scanning devices, ground-penetrating radar interferometry systems, or embedded micro-displacement sensors. After integrating these time sequences along the X, Y, and Z axes, a continuous displacement sequence in the three axes is formed. Subsequently, these sequences are segmented chronologically, and the displacement difference between adjacent time points in each time segment along the three axes is calculated, thus forming a set of vectors representing the three-dimensional disturbance trend. Each vector reflects the direction and magnitude of the target point's displacement in three-dimensional space per unit time. Next, using the superposition of three-dimensional geometric vectors, the spatial offset amplitude corresponding to each vector can be obtained, used to describe the disturbance intensity. Simultaneously, the system continuously tracks which time segments exhibit non-zero vector changes, and by statistically analyzing the duration of these continuous disturbance segments, the disturbance persistence is determined. For example, if an observation point moves more than 2 millimeters in the X direction for 8 minutes starting from a certain time, and there are corresponding disturbances in the Y and Z directions, it can be identified as a continuous disturbance event and recorded in the disturbance trajectory vector set. This process not only captures the instantaneous amplitude of boundary morphological changes, but also fully characterizes the temporal evolution of the disturbance, which is the basis for subsequent spatial anomaly judgment and coupling modeling.

[0032] Based on the spatial coordinates of the disturbance trajectory, the ventilation reversal rate data corresponding to the disturbance point in the ventilation monitoring network is called to extract the segment boundary of the ventilation direction change, and the spatial path of the disturbance trajectory is matched with the ventilation reversal segment to identify the overlap range and matching degree between the trajectory point and the ventilation reversal segment. This process first requires acquiring airflow direction data related to the disturbance trajectory using a ventilation monitoring network. The ventilation monitoring network is a spatial sensing system composed of wind speed and direction sensors distributed at different spatial points in the goaf, capable of recording real-time changes in airflow direction. Based on the three-dimensional spatial coordinates of the disturbance trajectory, the geographical location of the disturbance point can be accurately located and spatially correlated with each sensor point in the ventilation monitoring network. By comparing the distance, relative angle, and deployment density between the disturbance point and the ventilation monitoring points, the ventilation reversal rate data closest to or most representative of each disturbance point can be selected. The time points when the ventilation direction changes from forward to reverse or vice versa within adjacent time periods can be extracted, thus defining the boundaries of ventilation direction change segments. Based on this, a spatial boundary model of the ventilation reversal segment is constructed, and the point set of the disturbance trajectory along the spatial path is compared one-to-one with the ventilation reversal segment. The comparison process can be calculated using a spatial overlap rate index, such as counting how many points in the disturbance trajectory vector fall into the ventilation reversal segment, and using this to assess the degree of matching between the disturbance trajectory and the ventilation change segment. If most disturbance points are located within the ventilation reversal zone, the disturbance trajectory is considered to be highly affected by ventilation disturbances; otherwise, the impact is relatively small. For example, in a propulsion unit, if 90% of the points on a certain disturbance trajectory overlap with the ventilation reversal zone, the trajectory can be considered significantly controlled by ventilation disturbances. This judgment provides a reliable spatial linkage basis for subsequent fusion structure construction. This method enables the spatial linkage capture between disturbance characteristics and airflow disturbances, effectively improving the ability to characterize local oxygen transport variability in simulations.

[0033] The matching results are normalized, and the processed matching degree is jointly mapped with the three-axis perturbation trajectory. Based on the spatial overlap relationship between the trajectory spatial path and the ventilation reverse distribution, a fusion structure of the three-dimensional perturbation offset trajectory and ventilation reverse distribution features is constructed for subsequent generation of boundary morphology uncertainty descriptor.

[0034] When processing the matching results between spatial paths and ventilation reversal distributions, the spatial overlap between each disturbance trajectory point and the ventilation reversal segment must first be quantified into a matching value. This matching value can be jointly evaluated using multiple spatial relationship indicators, such as whether the disturbance point is located within the ventilation reversal segment, the distance between the point and the segment boundary, and the angle between the trajectory and the direction of wind change. Since these indicators have different dimensions and distributions, to avoid computational bias, all matching values ​​need to be normalized so that their values ​​are uniformly distributed between 0 and 1, ensuring the balance of subsequent mappings. The normalized matching values ​​are then jointly mapped with the disturbance trajectory data along the three axes. That is, at each disturbance point, its spatial location, displacement direction and amplitude, and matching degree are used to construct a high-dimensional attribute vector. Through this joint mapping operation, a spatial fusion expression model that simultaneously includes displacement behavior characteristics and ventilation disturbance response characteristics can be established in a three-dimensional spatial coordinate system—a fusion structure of three-dimensional disturbance offset trajectory and ventilation reversal distribution characteristics. In the fused structure, each trajectory unit records both the time and direction of the disturbance, and also expresses the sensitivity of that disturbance to the response in the ventilation-opposite zone. For example, if the matching value of a disturbance vector is 0.9, and its spatial displacement amplitude is 0.6 meters, pointing upwind, then the disturbance can be presumed to have significant ventilation blockage sensitivity. This structure not only unifies spatial disturbances and airflow anomalies within the same computational framework, but also provides a precise spatial coupling basis for the quantitative characterization of boundary morphology uncertainties. This approach effectively avoids the causal separation caused by the separate modeling of disturbances and ventilation in traditional analyses, thereby improving the model's analytical capability for changes in oxygen transport states under nonlinear boundary variations.

[0035] In this embodiment, S202 specifically refers to: Traverse each spatial point in the fusion structure of the three-dimensional disturbance offset trajectory and ventilation reverse distribution characteristics, extract the disturbance amplitude, disturbance duration and ventilation reverse degree data of each spatial point, construct a trend vector based on the disturbance direction and the ventilation reverse change direction, calculate the spatial change trend angle between the disturbance direction and the ventilation reverse change direction, and determine the spatial overlap ratio based on the overlap ratio of the angle and the spatial point in the ventilation change section. When analyzing the fusion structure of the 3D disturbance offset trajectory and the ventilation reversal distribution characteristics, it is first necessary to extract the disturbance amplitude, disturbance duration, and ventilation reversal degree for each spatial point. The disturbance amplitude refers to the absolute value of the displacement of the point in the 3D coordinate system, which can be obtained by calculating the magnitude of the three-axis displacement vector; the disturbance duration represents the length of time during which the point continuously exhibits significant displacement changes, usually measured by the total duration of consecutive threshold segments in the time series; the ventilation reversal degree is calculated based on the depth and frequency of the region where the wind direction reverses in the ventilation data. Next, two trend vectors need to be constructed: one is the disturbance direction vector, representing the dominant movement direction of the disturbance trajectory; the other is the ventilation reversal change direction vector, representing the main spatial direction of the ventilation reversal trend. Calculating the angle between these two vectors in 3D space, i.e., the spatial trend angle, can reveal whether the disturbance significantly intersects or opposes the ventilation reversal direction. For example, if the angle is close to 180 degrees, it indicates that the disturbance and the ventilation reversal direction overlap in opposite directions, which may lead to the ventilation path being blocked by the disturbance. To further quantify this spatial consistency, it is necessary to calculate the spatial overlap ratio by combining the actual coverage proportion of spatial points within the ventilation change zone. This ratio represents the proportion of the length of the ventilation reversal zone where the trajectory point is located to the total length of the disturbance path. If 80% of the length of a disturbance path is within the ventilation reversal zone, it indicates a high correlation with ventilation anomalies. Using the angle and overlap ratio as indicators, we can not only determine whether the disturbance and ventilation change match in direction and location, but also quantify the strength of this match, providing a spatial coupling basis for subsequent calculations of coupling strength. This approach effectively avoids the inability of a single indicator to accurately reflect the true relationship between spatial disturbances and ventilation distribution, thereby improving the accuracy of understanding nonlinear boundary response behavior.

[0036] Based on the spatial overlap ratio and the angle of change trend of each spatial point, a spatial coupling strength calculation function is constructed according to the product of cosine similarity and overlap, outputting the coupling strength value of the corresponding spatial point, and summarizing the coupling strength set of all spatial points to quantitatively characterize the reverse coupling characteristics of disturbance and ventilation. When evaluating the coupling relationship between three-dimensional disturbances and ventilation reversal, a spatial coupling strength calculation function needs to be constructed for each spatial point. This function models the spatial overlap ratio and the angle of change trend as core variables. The spatial overlap ratio refers to the proportion of coverage between the disturbance trajectory and the ventilation reversal segment on the spatial path, reflecting positional consistency, while the angle of change trend is the angular difference between the disturbance direction and the ventilation reversal direction, reflecting directional consistency. To numerically integrate information from these two dimensions, the spatial coupling strength calculation function can be constructed by multiplying the overlap ratio by the cosine similarity of the trend angle. The cosine similarity maps the angle to between 0 and 1; the closer the value is to 1, the more consistent the directions. The overlap ratio itself is a normalized value, also between 0 and 1. The product operation ensures that the coupling strength value remains within the standardized range, facilitating subsequent analysis. Taking a spatial point as an example, if its angle is 30 degrees (cosine value approximately 0.87) and the overlap ratio is 0.9, then the spatial coupling strength of this point is 0.783, indicating that the disturbance at this point has a high coupling with the ventilation reversal path. Repeating this process allows for the calculation of coupling strength at each spatial point within the fused structure, ultimately forming a spatial coupling strength set. This set comprehensively reflects the potential interference levels of disturbances in different regions on the ventilation path. This structured and quantitative approach more accurately captures the dynamic changes in oxygen diffusion paths caused by local disturbances in the goaf, providing crucial support for subsequent boundary morphology uncertainty classification. The spatial coupling strength calculation function is a core method for mapping multidimensional spatial behavior into a single comparable quantity, preserving the physical meaning of the model while facilitating algorithmic processing with controllable accuracy.

[0037] The range of global maximum and mean values ​​is extracted from the coupling strength set, and combined with the maximum offset amplitude and minimum disturbance duration in the three-dimensional disturbance offset trajectory, and uniformly mapped to the boundary disturbance spatial domain. The boundary morphology uncertainty description body with closed structure and complete numerical data is output in the form of a spatial grid.

[0038] When constructing a boundary morphology uncertainty descriptor, two key statistics need to be extracted from the spatial coupling strength set: the global maximum and the mean range. The global maximum reflects the extreme point of the inverse coupling between disturbance and ventilation, representing the most significant uncertainty area; the mean range is used to represent the overall strength benchmark of normal disturbance behavior. These two statistics will serve as the normalized reference boundary for subsequent spatial mapping. Next, the maximum offset amplitude and minimum disturbance duration are obtained from the three-dimensional disturbance offset trajectory. The former is used to characterize the extreme scale of spatial offset, while the latter represents whether the disturbance change can continue to the lower limit of the time that affects the airflow structure. After normalizing these four indicators, they are uniformly mapped to the boundary disturbance spatial domain and projected according to the set spatial grid division method. The combined scores of each indicator are filled grid by grid to form a boundary morphology uncertainty descriptor with spatial continuity and quantitative integrity. Taking an actual coal mine advance area as an example, if the maximum coupling strength of a certain area is 0.85, the maximum offset amplitude is 1.2 meters, and the minimum duration is 5 minutes, its score after mapping is much higher than that of the surrounding areas, and the grid of that area will be marked as high uncertainty level. The use of spatial grid output facilitates structural integration with subsequent simulation calculations and provides high-resolution basic data support for the visualization, early warning, and dynamic tracking of numerical fields, thereby enabling quantitative expression and intelligent analysis of complex boundary deformations in goaf areas.

[0039] S3. Establish an oxygen diffusion control trigger set based on the uncertainty level of the boundary morphology, and associate the boundary disturbance state in the oxygen diffusion control trigger set with the diffusion parameters to realize the dynamic adjustment of the simulated boundary oxygen diffusion conditions. In this embodiment, S3 specifically refers to: Based on the uncertainty level of the boundary morphology, an oxygen diffusion regulation trigger set is constructed. The boundary disturbance states are set as low disturbance state, medium disturbance state and high disturbance state respectively. Predefined diffusion regulation parameters are assigned to each boundary disturbance state, including the adjustment amplitude of oxygen concentration diffusion coefficient and the correction value of local airflow velocity. When constructing the oxygen diffusion regulation trigger set, it is necessary to first classify the disturbance characteristics within the propulsion unit into three levels: low disturbance, medium disturbance, and high disturbance, based on the uncertainty level of the boundary morphology. This process can be achieved by extracting the numerical ranges of spatial offset amplitude, disturbance duration, and coupling strength recorded in the boundary morphology uncertainty descriptor for classification. After the classification is completed, a set of corresponding oxygen diffusion regulation parameters is set for each level to adjust the boundary conditions during the simulation. Among them, the adjustment amplitude of the oxygen concentration diffusion coefficient can be set according to the degree of closure corresponding to the actual disturbance level. For example, a lower diffusion coefficient is set in the high disturbance state to simulate the situation where oxygen is difficult to enter. The local airflow velocity correction value is used to reflect the distortion effect of the disturbance on the ventilation path. By scaling down the original wind speed data or correcting the vector direction, the simulation accuracy of nonlinear boundary behavior can be improved. This mechanism enables dynamic adjustment of input parameters, keeping the boundary conditions synchronized with the actual evolution state during the numerical simulation.

[0040] The oxygen diffusion regulation trigger set is a data structure used to manage the mapping relationship between perturbation states and diffusion parameters. Its core purpose is to achieve responsive adjustment of simulated boundary conditions. This structure categorizes perturbation states into low, medium, and high perturbation states and binds them to corresponding diffusion parameters through mapping logic. The predefined diffusion regulation parameters are a set of directly callable numerical variables. The oxygen concentration diffusion coefficient adjustment amplitude controls the diffusion rate of gas in the boundary region, typically set to different reduction factors through a perturbation level mapping table; for example, a 10% reduction for low perturbation, a 30% reduction for medium perturbation, and a 60% reduction for high perturbation. The local airflow velocity correction value is used to correct the wind speed vector, which can be fine-tuned using wind speed sensor data at perturbation points to reflect the true wind field state after the perturbation. The overall goal is to achieve precise dynamic linkage between boundary perturbations and diffusion behavior in the model through this structure.

[0041] A one-to-one correspondence is established between the boundary perturbation states of each oxygen diffusion regulation trigger and the corresponding diffusion regulation parameters. A parameter lookup table is established using the propulsion unit code as an index. The orderly association between the boundary perturbation states and the diffusion parameters is achieved through this parameter lookup table. When establishing a one-to-one correspondence between boundary disturbance states and diffusion regulation parameters, the first step is to bind each disturbance state with a preset combination of diffusion regulation parameters based on the boundary disturbance states already determined by the propulsion unit. This process is achieved by establishing a mapping table, using the classification label of each disturbance state as the primary key and the corresponding oxygen concentration diffusion coefficient adjustment amplitude and local airflow velocity correction value as value items, forming a clear parameter binding structure. Subsequently, a parameter retrieval table is established using the propulsion unit code as the index primary key, enabling the system to quickly retrieve the corresponding boundary disturbance state and call the corresponding regulation parameters based on the current propulsion position during simulation execution. This structured binding method significantly improves the timeliness and accuracy of boundary condition updates during simulation, enabling rapid response and dynamic adaptation to oxygen diffusion conditions under changes in disturbance states.

[0042] Boundary perturbation states are the hierarchical classification of perturbation evolution behavior within a propulsion unit, typically categorized into low, medium, and high perturbation states, representing the degree of influence of boundary geometry changes on oxygen transport. Diffusion regulation parameters include the adjustment amplitude of the oxygen concentration diffusion coefficient and the correction value for local airflow velocity. The former controls the diffusion capacity of oxygen within the simulation region, while the latter corrects the wind speed vector within the perturbed region. The propulsion unit code is a unique identifier for the propulsion location, usually generated jointly through time series and spatial coordinates. The parameter lookup table is a data structure indexed by the propulsion unit code and containing diffusion regulation parameters, supporting efficient data retrieval and updates, and playing a crucial role in rapidly switching simulation boundary conditions in three-dimensional numerical simulations. This design highly couples spatial perturbation information with the oxygen diffusion model, providing structural support for dynamic adjustments.

[0043] Based on the parameter retrieval table, the simulated boundary configuration is dynamically interpolated and updated. When the propulsion unit enters the simulation iteration stage, the boundary disturbance state and diffusion parameters corresponding to the current code are automatically extracted, and the simulated boundary oxygen diffusion conditions are updated in real time to keep them synchronized with the disturbance state.

[0044] In the process of dynamically updating the simulated boundary configuration through interpolation, the boundary disturbance state and diffusion adjustment parameters corresponding to the propulsion unit code can be retrieved in real time based on the established parameter lookup table. When the simulation enters the iteration stage of a certain propulsion unit, the system first identifies the coding information of the current unit, and then extracts the disturbance state category matching the code, the corresponding oxygen concentration diffusion coefficient adjustment amplitude, and the local airflow velocity correction value from the lookup table. To avoid numerical instability caused by abrupt changes in boundary conditions due to changes in the disturbance state level, a time-weighted coefficient can be introduced to dynamically interpolate the disturbance parameters, calculate the transition values ​​of diffusion parameters under disturbance states between consecutive propulsion units, and achieve smooth transition and gradual adjustment of the simulated boundary conditions. This method ensures that the boundary oxygen diffusion conditions always remain dynamically consistent with the disturbance evolution trend during the propulsion process, improving the accuracy of the simulation results in response to locally ventilated areas.

[0045] Dynamic interpolation updates refer to the process of interpolating diffusion parameters between adjacent propulsion units or time steps based on discrete changes in the disturbance state, thereby generating continuous simulated boundary input conditions. This technique can be implemented based on linear interpolation, spline interpolation, or higher-order weighted interpolation algorithms, and the interpolation step size can be adjusted in conjunction with the frequency of disturbance state changes to ensure the stability of the simulation process and the physical continuity of the boundary conditions. The simulated boundary oxygen diffusion conditions typically include the spatial distribution function of the oxygen concentration diffusion coefficient and the airflow velocity vector. Its update frequency is affected by both the simulation step size and the rate of disturbance change. Therefore, implementing parameter updates through a dynamic interpolation mechanism can not only improve the robustness of numerical calculations but also enhance the ability to resolve local disturbance potentials. In regions with strong structural continuity, the interpolation smoothness can be appropriately reduced to improve computational efficiency, while in regions with severe disturbances, a high-precision interpolation strategy is required to ensure both simulation stability and accuracy.

[0046] S4. The dynamically adjusted simulated boundary oxygen diffusion conditions are introduced into the coal oxidation reaction calculation process. The oxygen diffusion changes and oxidation reaction parameters are linked by constructing an oxygen consumption linkage factor set. In this embodiment, S4 specifically refers to: The oxygen concentration diffusion coefficient and local airflow velocity parameters contained in the dynamically adjusted simulated boundary oxygen diffusion conditions are extracted. An oxygen diffusion input sequence is established according to the coding order of the propulsion unit and the simulation time step. The oxygen diffusion input sequence is then spatially mapped to the boundary nodes in the coal oxidation reaction calculation process. To ensure accurate input of dynamically adjusted simulated boundary oxygen diffusion conditions into the coal oxidation reaction calculation process, two key parameters—oxygen concentration diffusion coefficient and local airflow velocity—must first be extracted from the simulated boundary. The oxygen concentration diffusion coefficient reflects the ability of oxygen molecules to diffuse in porous media, while the local airflow velocity characterizes the flow rate of external air entering the model through the boundary. These two parameters dynamically adjust with changes in boundary perturbation states. Therefore, the parameters need to be progressively organized according to the spatial coding order of the propulsion units and the simulation time steps to form a time-series oxygen diffusion input sequence. Next, the parameter values ​​at each time step in this input sequence are spatially matched one-to-one with the corresponding boundary nodes in the simulation mesh, ensuring that each boundary node has a unique diffusion and airflow input at each time step, thus enabling real-time driving of the oxidation reaction model by oxygen diffusion conditions. For example, when a propulsion unit experiences a high-level perturbation state at time step 20, its boundary nodes will receive corresponding low diffusion coefficient and low airflow velocity values, thereby suppressing oxygen input in that region.

[0047] The dynamically adjusted simulated boundary oxygen diffusion conditions mainly consist of two technical features: the oxygen concentration diffusion coefficient and the local airflow velocity parameter. The oxygen concentration diffusion coefficient reflects the rate of oxygen propagation in the coal seam or residual coal deposits, and is influenced by factors such as porous structure, coal quality, and temperature. It can be estimated through experimental fitting or gas transport equations. The local airflow velocity parameter represents the velocity of air entering the model at the boundary, and is typically set in conjunction with ventilation monitoring and boundary disturbance levels. The coal oxidation reaction calculation process is the core step in solving for the oxygen consumption rate based on mass conservation and a kinetic reaction model. This process involves multiple reaction variables, such as oxygen concentration, temperature, and reaction order, and the solution directly determines the spatial distribution of spontaneous combustion risk. After the oxygen diffusion input sequence is spatially mapped to the boundary nodes, the oxygen diffusion conditions can be sequentially input into the reaction model step-by-step, enabling dynamic control of the oxidation reaction by oxygen supply conditions, thereby improving simulation accuracy and the ability to identify local heat accumulation risks.

[0048] In the calculation of coal oxidation reaction, the oxygen concentration diffusion coefficient and local gas flow velocity parameters corresponding to each time step in the oxygen diffusion input sequence are input into the reaction region model. The rate of change of oxygen content, the rate of change of oxygen partial pressure gradient and the time derivative of oxygen concentration in a unit volume region are calculated to form a set of basic reaction input factors characterizing the changes in oxygen diffusion. In the calculation of coal oxidation reaction, to achieve the linkage control between simulated boundary oxygen diffusion conditions and oxidation reaction parameters, key parameters in the oxygen diffusion input sequence need to be introduced into the reaction region model step by step over time. Specifically, the oxygen concentration diffusion coefficient and local gas flow velocity parameters for each time step can be assigned to the corresponding grid cells in the reaction region. These parameters are then used to calculate the rate of change of oxygen content per unit volume. The rate of change of oxygen content can be solved using the oxygen mass conservation equation, while the rate of change of oxygen partial pressure gradient is derived from the ratio of the oxygen partial pressure difference to the distance between adjacent grid cells. The time derivative of oxygen concentration represents the trend of oxygen concentration change over time. These physical quantities all reflect the specific impact of boundary perturbations on oxygen transport behavior. Based on this, a set of oxygen diffusion response values ​​generated at each time step is collected into a time series to support the subsequent construction of the linkage mechanism. For example, under high perturbation conditions, the rate of change of oxygen content per unit volume may decrease significantly, and the corresponding oxygen partial pressure gradient tends to flatten, reflecting a local environment with limited oxygen diffusion, thereby guiding the reaction model to lower the oxidation rate.

[0049] The reaction region model is a spatial computational framework for simulating the physical processes of coal oxidation reactions. Its core involves the simultaneous solution of the temperature field, gas concentration field, and reaction rate field. The rate of change of oxygen content per unit volume represents the total change in oxygen concentration per unit time due to diffusion and reaction, serving as a direct indicator of changes in oxygen supply levels. The rate of change of oxygen partial pressure gradient is used to assess the driving force of oxygen in space, influenced by both the boundary input oxygen partial pressure and the internal consumption rate. The time derivative of oxygen concentration reflects the trend of oxygen concentration change over time; combined with time series analysis, it can be used to identify whether fluctuations in oxygen input exhibit a trend of decline or increase. These parameters together form the basic reaction input factor set, which is the core link connecting oxygen diffusion and reaction behavior in the simulation. By calculating and continuously updating this factor set in real time, the reaction model can dynamically respond to external disturbances, achieving more accurate predictions of oxidation reaction rates and heat release trends.

[0050] Based on the basic reaction input factor set, an oxygen consumption linkage factor set is constructed. The oxygen content change rate, oxygen partial pressure gradient change rate, and oxygen concentration time derivative are dynamically mapped to the reaction rate, reaction order, and oxygen consumption function parameters in the coal oxidation reaction model, respectively. Through the oxygen consumption linkage factor set, the linkage processing of oxygen diffusion changes and oxidation reaction parameters is realized, and the oxygen consumption behavior in the reaction process is updated synchronously.

[0051] To achieve dynamic adjustment of coal oxidation reaction parameters by changes in oxygen diffusion, an oxygen consumption linkage factor set needs to be constructed based on the basic reaction input factor set. This process first extracts the oxygen content change rate, oxygen partial pressure gradient change rate, and oxygen concentration time derivative as independent variables, establishing a mapping relationship between these variables and the core control parameters in the coal oxidation reaction model. The reaction rate is usually directly related to the local oxygen concentration; therefore, the oxygen content change rate can be linearly or nonlinearly mapped to the local reaction rate, so that an increase in oxygen concentration corresponds to an increased oxidation rate. The oxygen partial pressure gradient change rate is mainly used to adjust the dynamic parameters of the reaction order, thereby simulating the nonlinear influence of oxygen-limited conditions on the reaction path. The oxygen concentration time derivative is linked to the oxygen consumption function parameters, which can correct the coupling strength or critical value setting in the oxygen consumption process. When these input variables change, the corresponding reaction model parameters are automatically updated, allowing the model to synchronously reflect the actual diffusion state in the perturbed region. For example, in regions where the oxygen partial pressure decreases rapidly, the reaction order in the model is automatically increased to simulate higher-order reaction paths under diffusion-limited conditions, thus more accurately predicting the reaction trend.

[0052] The oxygen consumption linkage factor set is a dynamic mapping set used to synchronize changes in input boundary conditions with reaction parameters in a coal oxidation reaction model. Within this factor set, the oxygen content change rate reflects the oxygen supply intensity and is the primary control index for regulating the reaction rate; the oxygen partial pressure gradient change rate measures the spatial variation of the oxygen diffusion driving force and can be used to dynamically adjust the series or control terms in the reaction kinetic model; the oxygen concentration time derivative reflects the temporal trend of oxygen input and determines the response rate of the reaction system to external oxygen disturbances. In numerical implementation, these factors are typically established by associating them with the target reaction parameters through lookup tables, fitting functions, or neural network models, forming a reversible control mechanism. After linkage processing, the reaction model can automatically update the oxidation reaction behavior at each time step based on the current linkage factor state, thereby maintaining high-precision oxygen consumption prediction capabilities under multiple disturbance conditions. This technical structure ensures a closed-loop information system between the oxygen diffusion state and the coal oxidation reaction process in the simulated environment, improving the model's adaptability and controllability to disturbance responses.

[0053] S5. Based on the changes in the oxygen consumption linkage factor set, a simulated control equilibrium sequence is constructed, and the calculation frequency, spatial resolution, and calculation step size of the propulsion unit are dynamically controlled according to the response change rate in the simulated control equilibrium sequence.

[0054] In this embodiment, S5 specifically refers to: Based on the changes in the oxygen consumption linkage factors, the numerical differences of the oxygen content change rate, oxygen partial pressure gradient change rate, and oxygen concentration time derivative in continuous time steps are extracted. The change rate of each factor is calculated, and a response change rate sequence is constructed in units of propulsion units. The response change rate sequence in each time step is normalized to form a simulated control equilibrium sequence for dynamic control analysis. To construct the simulated control equilibrium sequence, it is necessary to extract continuous values ​​of three key factors—the rate of change of oxygen content, the rate of change of oxygen partial pressure gradient, and the time derivative of oxygen concentration—from the oxygen consumption linkage factor set at different time steps. The rate of change for each factor is obtained by calculating the numerical difference between two adjacent time steps. For example, in a certain propulsion unit, if the rate of change of oxygen content is 0.12 and 0.18 in two consecutive simulation time steps, its rate of change is 0.06. Similarly, the rate of change of the three factors for each propulsion unit is calculated throughout the entire simulation period, resulting in a three-dimensional time series. The rate of change corresponding to each time step in this series is arranged in order of propulsion unit, forming a response rate of change sequence at the propulsion unit level. This sequence is then normalized to ensure that comparisons and analyses can be performed at a uniform scale across different propulsion units. The normalization method can employ min-max normalization, mapping different physical quantities to a uniform interval between 0 and 1, thereby forming a simulated control equilibrium sequence reflecting the trend of oxygen consumption response changes, which is used for subsequent dynamic adjustment of calculation parameters.

[0055] The oxygen content change rate represents the rate of change of oxygen concentration within a unit volume area per unit time, and is an important indicator for measuring the degree of dynamic response of oxygen diffusion. The oxygen partial pressure gradient change rate describes the magnitude of the change of oxygen partial pressure gradient over time over a unit distance, revealing the changing trend of oxygen transport direction and speed. The oxygen concentration time derivative is the derivative of oxygen concentration with time, used to characterize the instantaneous trend of oxygen concentration change. The propulsion unit is the basic computational unit for spatially dividing the simulation region, and its encoding order determines the execution path of the simulation calculation. The response change rate sequence is a collection of rates of various oxygen consumption factors arranged in time, used to reflect the dynamic evolution trend of oxygen consumption in a local area. The simulation regulation equilibrium sequence is a unified index sequence formed by normalizing the response change rate sequence, serving as the basis for judging the dynamic regulation strategy in the simulation system, and driving the model to finely allocate computational resources at different evolution stages.

[0056] Based on the response change rate in the simulated control equilibrium sequence, a control mapping relationship is established between the response change rate and the calculation frequency, spatial resolution and calculation step size of the propulsion unit. The propulsion units are classified according to the range of response change rate, and corresponding calculation frequency adjustment value, spatial resolution adjustment value and calculation step size adjustment value are assigned to different levels of propulsion units, thus forming a set of control parameters that match the simulated control equilibrium sequence. To achieve dynamic control of the propulsion unit's computational process, a control mapping relationship needs to be established based on the response change rate of each propulsion unit in the simulation control equilibrium sequence. Specifically, the distribution range of the response change rate across all propulsion units is first statistically analyzed, and the propulsion units are then classified according to a predefined grading strategy (e.g., low-speed, medium-speed, high-speed). For example, units with a response change rate below 0.2 are classified as Level 1, the intermediate range as Level 2, and those above 0.6 as Level 3. Subsequently, a set of control parameters is preset for each level, including computation frequency adjustment values, spatial resolution adjustment values, and calculation step size adjustment values. The computation frequency adjustment value controls the time interval between propulsion unit participation in the simulation calculation; for example, Level 1 units can be set to calculate once every 5 steps, while Level 3 units participate in the calculation every step. The spatial resolution adjustment value controls the mesh density, with high-response units using finer spatial meshes. The calculation step size adjustment value determines the accuracy of the simulation time progression in each iteration, reflecting the model's responsiveness to instantaneous changes. By mapping the response change rate of each propulsion unit to these control parameters, a complete set of control parameters can be constructed, realizing a dynamic computational control mechanism that matches the simulation control equilibrium sequence.

[0057] The response rate of change is an index calculated from the difference in oxygen consumption linkage factors within continuous time steps, used to measure the intensity of the dynamic coupling between oxygen diffusion and coal oxidation reaction within a propulsion unit. The computation frequency represents the time step interval at which the model updates the state of a propulsion unit, and is a key parameter controlling the allocation of computational resources. Spatial resolution refers to the level of detail in the 3D mesh generation of the propulsion unit; a higher value indicates a finer mesh within the unit, which helps improve simulation accuracy. The calculation step size is the smallest unit of time advancement in each iteration of the simulation; a smaller step size makes the model more sensitive to instantaneous responses. Adjusting the mapping relationship is the process of establishing a mathematical relationship between the response rate of change and these simulation parameters, typically achieved through interval-level parameter binding. The control parameter set is a collection storing the computation frequency, spatial resolution, and calculation step size in each propulsion unit, used to dynamically control the operation of the entire simulation system, thereby improving computational efficiency while maintaining accuracy.

[0058] The set of control parameters is applied to the simulation system. Using the propulsion unit as an index, the corresponding control parameters are automatically selected based on the rate of change of the response in the simulation control equilibrium sequence corresponding to the current time step. The calculation frequency, spatial resolution and calculation step size of the propulsion unit are dynamically controlled, and adaptive allocation of computing resources and synchronous adjustment of simulation accuracy are realized during the simulation iteration process.

[0059] To achieve dynamic adaptation of computational resources and simulation accuracy during the propulsion process, the set of control parameters needs to be applied in real time to the simulation configuration of each propulsion unit. Specifically, before the simulation system enters the computation phase at each time step, the response rate of the propulsion unit is extracted from the simulation control equilibrium sequence corresponding to the current time step, using the propulsion unit code as an index. The system matches this value with a preset control mapping relationship and automatically selects the computational frequency, spatial resolution, and calculation step size corresponding to the current response rate from the control parameter set. For example, when the response rate of a propulsion unit is in a high-level range, the system will automatically enable high-frequency computation, high-density spatial resolution, and short-step calculation configuration, allowing that unit to receive more resource investment and accuracy assurance during simulation iterations; while units with lower response rates will be assigned lower computational frequencies and moderate or coarse resolutions and step sizes. Through this mechanism, the simulation system can dynamically adjust the computational parameters according to the state of each unit during parallel computation of multiple propulsion units, thereby saving unnecessary computational resources while ensuring accuracy in critical areas, achieving optimal resource utilization and a response-driven computation strategy across the entire domain.

[0060] The control parameter set is a collection of multiple control parameter combinations generated based on the response change rate. Each set of parameters includes computation frequency, spatial resolution, and calculation step size, used to guide the specific behavior of the propulsion unit during the simulation process. The computation frequency represents the frequency of the propulsion unit's participation in the computational process over time, controlling whether the unit updates its state at each time step or executes at set intervals. Spatial resolution determines the fineness of the simulation mesh within the propulsion unit, directly impacting the local accuracy of the coal oxidation reaction and oxygen diffusion model. The calculation step size controls the granularity of time progression; a smaller value indicates more detailed simulation changes per unit time. The simulation control equilibrium sequence is a dual-indexed sequence built based on time and unit number, recording the response change rate of each propulsion unit at any time step, used for dynamically retrieving and adapting parameters. Dynamically applying the control parameter set to the simulation system achieves the goal of real-time allocation and adjustment of computational resources and simulation accuracy according to the reaction state, providing a technical foundation for the precise coupled modeling of oxygen diffusion and coal oxidation processes.

[0061] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented using software, the above embodiments can be implemented, in whole or in part, as a computer program product. A computer program product includes one or more computer instructions or computer programs. When the computer instructions or computer programs are loaded or executed on a computer, all or part of the processes or functions according to the embodiments of this application are generated. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable device. Computer instructions can be stored in a computer-readable storage medium or transmitted from one computer-readable storage medium to another. For example, computer instructions can be transmitted from one website, computer, server, or data center to another website, computer, server, or data center via wired or wireless means (e.g., infrared, wireless, microwave, etc.). A computer-readable storage medium can be any available medium that a computer can access or a data storage device such as a server or data center that includes one or more sets of available media. Available media can be magnetic media (e.g., floppy disks, hard disks, magnetic tapes), optical media (e.g., DVDs), or semiconductor media. Semiconductor media can be solid-state drives.

[0062] It should be understood that in the various embodiments of this application, the order of the above-mentioned processes does not imply the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of this application.

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

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

[0065] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment according to actual needs.

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

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

Claims

1. A three-dimensional dynamic numerical simulation method for spontaneous combustion of coal in the goaf during working face advancement, characterized in that, Specifically, the following steps are included: S1. By collecting the micro-seismic response value, ventilation reversal rate and stress release fluctuation value of compacted residual coal in the propulsion unit, a boundary disturbance identification structure is constructed. Based on the abnormal co-occurrence threshold in the boundary disturbance identification structure, it is determined whether the local subsidence of the goaf boundary space structure occurs during the propulsion process, and the area where the goaf boundary space structure has local subsidence is coded and marked. S2. Based on the coded and identified areas, and combined with the three-dimensional disturbance offset trajectory and ventilation reverse distribution characteristics within the boundary disturbance identification structure, a boundary morphology uncertainty description body is generated. Based on the spatial offset amplitude, disturbance duration, and coupling strength in the boundary morphology uncertainty description body, the boundary morphology uncertainty level is determined when a local subsidence of the goaf boundary spatial structure occurs during the advancement process. S3. Establish an oxygen diffusion control trigger set based on the uncertainty level of the boundary morphology, and associate the boundary disturbance state in the oxygen diffusion control trigger set with the diffusion parameters to realize the dynamic adjustment of the simulated boundary oxygen diffusion conditions. S4. The dynamically adjusted simulated boundary oxygen diffusion conditions are introduced into the coal oxidation reaction calculation process. The oxygen diffusion changes and oxidation reaction parameters are linked by constructing an oxygen consumption linkage factor set. S5. Based on the changes in the oxygen consumption linkage factor set, a simulated control equilibrium sequence is constructed, and the calculation frequency, spatial resolution, and calculation step size of the propulsion unit are dynamically controlled according to the response change rate in the simulated control equilibrium sequence.

2. The three-dimensional dynamic numerical simulation method for spontaneous combustion of coal in the goaf during working face advancement according to claim 1, characterized in that, S1 specifically includes the following steps: S101. Collect microseismic response values, ventilation reversal rate, and compacted residual coal stress release fluctuation values ​​within the propulsion unit, and establish microseismic response value sequences, ventilation reversal rate sequences, and compacted residual coal stress release fluctuation value sequences respectively. Synchronously organize the three types of sequences according to the spatial location of the propulsion unit to form the basic dataset for boundary disturbance identification structure. S102. The three types of sequences in the basic dataset are jointly arranged according to the advancement time sequence to form a boundary disturbance identification structure. The boundary disturbance identification structure uses the time overlap relationship between the sudden increase point of microseismic response value, the change point of ventilation reverse rate and the abnormal point of stress release fluctuation value of compacted residual coal as an abnormal co-occurrence reference index. By comparing the abnormal co-occurrence threshold, it is determined whether the local inward subsidence of the boundary space structure of the goaf occurs during the advancement process. S103. After determining that a local subsidence of the boundary space structure of the goaf has occurred, the main identifier is the advancing unit number, and the auxiliary identifier is the corresponding abnormal co-occurrence location in the boundary disturbance identification structure. The main identifier and the auxiliary identifier are combined to form a coded identifier, and the coded identifier is recorded as the area marking result of the local subsidence of the boundary space structure of the goaf during the advancing process.

3. The three-dimensional dynamic numerical simulation method for spontaneous combustion of coal in the goaf during working face advancement according to claim 2, characterized in that, S102 specifically refers to: Extract the timestamp information of the microseismic response value sequence, ventilation reverse rate sequence and compacted residual coal stress release fluctuation value sequence from the basic dataset, establish a unified time axis according to the advancement sequence, and arrange all observation point data of the three types of sequences onto the unified time axis. By screening out abrupt increases in the microseismic response value sequence, changes in the ventilation reverse rate sequence, and anomalies in the stress release fluctuation value sequence of compacted residual coal, the overlap of the three types of characteristic points on a unified time axis is used as an abnormal co-occurrence reference index for boundary disturbance identification structures according to the spatial coordinates of each propulsion unit. Based on a pre-set abnormal co-occurrence threshold, when the three types of feature points in the propulsion unit coincide on a unified time axis, it is determined that the propulsion unit has a local indentation of the boundary space structure of the goaf during the propulsion process, and the judgment result is written into the indentation marker of the boundary disturbance identification structure.

4. The three-dimensional dynamic numerical simulation method for spontaneous combustion of coal in the goaf during working face advancement according to claim 1, characterized in that, S2 specifically includes the following steps: S201. Extract the three-dimensional disturbance offset trajectory within the area corresponding to the coded identifier, construct a disturbance trajectory vector group through the three-axis displacement sequence, and perform time-by-time offset analysis on each vector group to calculate the spatial offset amplitude and disturbance duration; simultaneously extract ventilation reverse distribution features, spatially match the disturbance trajectory with the ventilation reverse section, and construct a fusion structure of three-dimensional disturbance offset trajectory and ventilation reverse distribution features. S202. Based on the disturbance amplitude, disturbance duration and ventilation reversal degree of each spatial point in the fusion structure, calculate the corresponding coupling strength according to the angle between the spatial overlap ratio and the changing trend. Combine the maximum displacement amplitude and minimum duration of the three-dimensional disturbance offset trajectory to generate a boundary morphological uncertainty description body with structural integrity and quantitative characteristics. S203. Based on the spatial offset amplitude, disturbance duration and coupling strength extracted from the boundary morphology uncertainty description, a scoring model is constructed. Feature weights are assigned to each and weighted calculations are performed to obtain the comprehensive score value of boundary morphology uncertainty. The comprehensive score value is compared with the pre-set uncertainty level range and classified into high-level, medium-level or low-level uncertainty level results, and bound to the coded identification area.

5. The three-dimensional dynamic numerical simulation method for spontaneous combustion of coal in the goaf during working face advancement according to claim 4, characterized in that, S201 specifically refers to: Based on the spatial region corresponding to the coded identifier, the sequence of displacement observation points arranged along the three axes is extracted, and integrated into a three-axis displacement sequence in chronological order. Based on the displacement changes of each axis in adjacent time periods, a disturbance trajectory vector group is constructed. The three-dimensional displacement changes are calculated vector by vector to obtain the spatial offset amplitude of the disturbance trajectory. The time periods in which the disturbance exists continuously are recorded to obtain the disturbance duration. Based on the spatial coordinates of the disturbance trajectory, the ventilation reversal rate data corresponding to the disturbance point in the ventilation monitoring network is called to extract the segment boundary of the ventilation direction change, and the spatial path of the disturbance trajectory is matched with the ventilation reversal segment to identify the overlap range and matching degree between the trajectory point and the ventilation reversal segment. The matching results are normalized, and the processed matching degree is jointly mapped with the three-axis perturbation trajectory. Based on the spatial overlap relationship between the trajectory spatial path and the ventilation reverse distribution, a fusion structure of the three-dimensional perturbation offset trajectory and ventilation reverse distribution features is constructed for subsequent generation of boundary morphology uncertainty descriptor.

6. The three-dimensional dynamic numerical simulation method for spontaneous combustion of coal in the goaf during working face advancement according to claim 4, characterized in that, S202 specifically refers to: Traverse each spatial point in the fusion structure of the three-dimensional disturbance offset trajectory and ventilation reverse distribution characteristics, extract the disturbance amplitude, disturbance duration and ventilation reverse degree data of each spatial point, construct a trend vector based on the disturbance direction and the ventilation reverse change direction, calculate the spatial change trend angle between the disturbance direction and the ventilation reverse change direction, and determine the spatial overlap ratio based on the overlap ratio of the angle and the spatial point in the ventilation change section. Based on the spatial overlap ratio and the angle of change trend of each spatial point, a spatial coupling strength calculation function is constructed according to the product of cosine similarity and overlap, outputting the coupling strength value of the corresponding spatial point, and summarizing the coupling strength set of all spatial points to quantitatively characterize the reverse coupling characteristics of disturbance and ventilation. The range of global maximum and mean values ​​is extracted from the coupling strength set, and combined with the maximum offset amplitude and minimum disturbance duration in the three-dimensional disturbance offset trajectory, and uniformly mapped to the boundary disturbance spatial domain. The boundary morphology uncertainty description body with closed structure and complete numerical data is output in the form of a spatial grid.

7. The three-dimensional dynamic numerical simulation method for spontaneous combustion of coal in the goaf during working face advancement according to claim 1, characterized in that, S3 specifically refers to: Based on the uncertainty level of the boundary morphology, an oxygen diffusion regulation trigger set is constructed. The boundary disturbance states are set as low disturbance state, medium disturbance state and high disturbance state respectively. Predefined diffusion regulation parameters are assigned to each boundary disturbance state, including the adjustment amplitude of oxygen concentration diffusion coefficient and the correction value of local airflow velocity. A one-to-one correspondence is established between the boundary perturbation states of each oxygen diffusion regulation trigger and the corresponding diffusion regulation parameters. A parameter lookup table is established using the propulsion unit code as an index. The orderly association between the boundary perturbation states and the diffusion parameters is achieved through this parameter lookup table. Based on the parameter retrieval table, the simulated boundary configuration is dynamically interpolated and updated. When the propulsion unit enters the simulation iteration stage, the boundary disturbance state and diffusion parameters corresponding to the current code are automatically extracted, and the simulated boundary oxygen diffusion conditions are updated in real time to keep them synchronized with the disturbance state.

8. The three-dimensional dynamic numerical simulation method for spontaneous combustion of coal in the goaf during working face advancement according to claim 1, characterized in that, S4 specifically refers to: The oxygen concentration diffusion coefficient and local airflow velocity parameters contained in the dynamically adjusted simulated boundary oxygen diffusion conditions are extracted. An oxygen diffusion input sequence is established according to the coding order of the propulsion unit and the simulation time step. The oxygen diffusion input sequence is then spatially mapped to the boundary nodes in the coal oxidation reaction calculation process. In the calculation of coal oxidation reaction, the oxygen concentration diffusion coefficient and local gas flow velocity parameters corresponding to each time step in the oxygen diffusion input sequence are input into the reaction region model. The rate of change of oxygen content, the rate of change of oxygen partial pressure gradient and the time derivative of oxygen concentration in a unit volume region are calculated to form a set of basic reaction input factors characterizing the changes in oxygen diffusion. Based on the basic reaction input factor set, an oxygen consumption linkage factor set is constructed. The oxygen content change rate, oxygen partial pressure gradient change rate, and oxygen concentration time derivative are dynamically mapped to the reaction rate, reaction order, and oxygen consumption function parameters in the coal oxidation reaction model, respectively. Through the oxygen consumption linkage factor set, the linkage processing of oxygen diffusion changes and oxidation reaction parameters is realized, and the oxygen consumption behavior in the reaction process is updated synchronously.

9. The three-dimensional dynamic numerical simulation method for spontaneous combustion of coal in the goaf during working face advancement according to claim 1, characterized in that, S5 specifically refers to: Based on the changes in the oxygen consumption linkage factors, the numerical differences of the oxygen content change rate, oxygen partial pressure gradient change rate, and oxygen concentration time derivative in continuous time steps are extracted. The change rate of each factor is calculated, and a response change rate sequence is constructed in units of propulsion units. The response change rate sequence in each time step is normalized to form a simulated control equilibrium sequence for dynamic control analysis. Based on the response change rate in the simulated control equilibrium sequence, a control mapping relationship is established between the response change rate and the calculation frequency, spatial resolution and calculation step size of the propulsion unit. The propulsion units are classified according to the range of response change rate, and corresponding calculation frequency adjustment value, spatial resolution adjustment value and calculation step size adjustment value are assigned to different levels of propulsion units, thus forming a set of control parameters that match the simulated control equilibrium sequence. The set of control parameters is applied to the simulation system. Using the propulsion unit as an index, the corresponding control parameters are automatically selected based on the rate of change of the response in the simulation control equilibrium sequence corresponding to the current time step. The calculation frequency, spatial resolution and calculation step size of the propulsion unit are dynamically controlled, and adaptive allocation of computing resources and synchronous adjustment of simulation accuracy are realized during the simulation iteration process.

Images