Mine gas concentration dynamic analysis method and system

Abstract

This invention relates to the field of gas concentration analysis technology, specifically a method and system for dynamic analysis of mine gas concentration, comprising the following steps: collecting airflow fluctuation parameters, extracting the wind speed change vector, mapping it with the concentration frequency electrical signal to construct a disturbance topology map, extracting the flow resistance attenuation coefficient to calculate the lag time constant, performing phase compensation translation on the electrical signal to generate a waveform leading-edge time series, using the spot distortion parameter to perform turbidity restoration to calculate the actual optical path attenuation parameter, and comparing and outputting dynamic concentration values. In this invention, by introducing a concentration time-series inversion logic based on airflow coupling characteristics, the dynamic airflow parameters are fitted with interference fringe parameters to construct a pressure coupling matrix. Gradient compensation is performed based on the cross-sectional wind resistance coefficient to eliminate gas chamber response delay errors. Micro-dust spot distortion parameters are extracted to establish an inverse analytical relationship between turbidity and attenuation index. The optical path detection benchmark is compensated and calibrated to output continuous concentration distribution coordinates.

Description

Technical Field

[0001] This invention relates to the field of gas concentration analysis technology, and in particular to a method and system for dynamic analysis of mine gas concentration. Background Technology

[0002] The field of gas concentration analysis technology involves the quantitative determination and continuous monitoring of gas composition and volume fraction within a specific space using physical or chemical principles. Traditional methods for dynamic analysis of mine gas concentration rely on a single sensor array for periodic sampling and the delivery of gas samples from pipelines to a chromatograph for reaction.

[0003] Traditional gas chamber sampling and measurement mechanisms, based on fixed-point timed sampling, are unable to capture the transient characteristics of gas diffusion caused by changes in the flow field of the underground ventilation network. When encountering transient large-volume gas outflows generated by coal mining operations and airflow pressure fluctuations in multi-branch roadways, the physical response delay inside the fixed gas chamber causes a serious lag in the measurement results. A single temperature compensation circuit is unable to eliminate the interference of high-concentration dust particles scattering and attenuating the optical path, and the output electrical signal is unable to accurately reflect the spatiotemporal variation of concentration distribution in deep working faces. Summary of the Invention

[0004] The purpose of this invention is to overcome the shortcomings of existing technologies and to propose a method and system for dynamic analysis of mine gas concentration.

[0005] To achieve the above objectives, the present invention adopts the following technical solution: a method for dynamic analysis of mine gas concentration, comprising the following steps: S1: Collect airflow fluctuation parameters of the roadway cross section using a mine microwave anemometer, extract the wind speed change vector from the airflow fluctuation parameters, and perform nonlinear spatiotemporal mapping operation on the wind speed change vector and the original concentration frequency electrical signal collected by the gas detection probe to construct an environmental flow field disturbance topology map. S2: Extract the absolute pressure distribution gradient of discrete nodes in the environmental flow field disturbance topology map, calculate the air flow resistance attenuation coefficient based on the absolute pressure distribution gradient, substitute the air flow resistance attenuation coefficient into the transport mass conservation partial differential equation for parameter iteration calculation, and calculate the transient lag time constant. S3: Based on the transient lag time constant, the original concentration frequency electrical signal is phase-reverse compensated and shifted to generate a waveform leading edge time series. The dust turbidity of the waveform leading edge time series is restored and analyzed using the light spot geometric distortion parameter collected by the scattering dust meter, and the actual optical path attenuation parameter is calculated. S4: Obtain the pre-configured baseline dust-free absorbance difference value, perform differential calculation and comparison between the actual optical path attenuation parameter and the baseline dust-free absorbance difference value, and output the dynamic gas concentration value based on the differential calculation and comparison result. As a further aspect of the present invention, step S1 specifically comprises: S11: Obtain the airflow fluctuation parameters at the roadway cross-section collected by the mine microwave anemometer, perform high-frequency filtering and noise reduction processing on the airflow fluctuation parameters to filter out the periodic mechanical vibration interference noise caused by the mine ventilation fan, and extract the wind speed change vector in the three-dimensional spatial direction that reflects the local wind field turbulence characteristics to generate a wind speed feature sequence. S12: Acquire the original concentration frequency electrical signal collected by the gas detection probe, use the wind speed feature sequence and the original concentration frequency electrical signal as multivariate time series input items, perform the nonlinear spatiotemporal mapping operation on the multivariate time series input items, calculate the cross-information entropy between different physical measuring points at different time delay scales, and generate a spatiotemporal correlation matrix. S13: Extract the matrix elements in the spatiotemporal correlation matrix whose values ​​exceed the preset structure retention threshold as the weight values ​​of the directed edges, map each physical monitoring point in the well to an independent connected node in the network, and establish a spatial connectivity network between the independent connected nodes based on the weight values ​​of the directed edges to construct the environmental flow field disturbance topology map. As a further aspect of the present invention, step S2 specifically includes: S21: Extract the absolute pressure monitoring data and three-dimensional spatial coordinate information of each discrete node in the environmental flow field disturbance topology map, calculate the ratio of the absolute pressure difference between adjacent discrete nodes to the Euclidean distance in physical space, and generate the absolute pressure distribution gradient. S22: Obtain the pre-determined frictional resistance coefficient of the tunnel wall and the fluid dynamic viscosity parameters, calculate the degree of local energy loss of the gas during the flow process based on the absolute pressure distribution gradient, and calculate the air flow resistance attenuation coefficient in combination with the equivalent hydraulic diameter of the tunnel. S23: Obtain the partial differential equation of transport mass conservation that reflects the comprehensive law of gas molecule diffusion and convection transport, substitute the air flow resistance attenuation coefficient as the updated damping variable into the partial differential equation of transport mass conservation, perform the parameter iteration operation until the calculated rate of change of time response parameter reaches a stable convergence state, and calculate the transient lag time constant. As a further aspect of the present invention, step S3 specifically comprises: S31: Obtain the transient lag time constant, determine the physical time delay span required for the gas to diffuse from the release source to the optical detection path based on the transient lag time constant, perform the phase reverse compensation shift on the original concentration frequency electrical signal based on the physical time delay span, filter out the phase deviation caused by the shift delay, and generate the waveform leading edge time series; S32: Obtain the geometric distortion parameters of the light spot on the detection surface of the optical receiver collected by the scattering dust meter, analyze the spatial scattering and blocking effect of suspended particulate matter on the detection beam characterized by the geometric distortion parameters of the light spot, and extract optical interference features; S33: Using the optical interference characteristics, perform baseline compensation and dust turbidity reduction analysis on the abnormal amplitude depression region caused by dust scattering in the waveform leading edge time series, deduct non-absorbent optical loss error, and calculate the actual optical path attenuation parameter. As a further aspect of the present invention, step S4 specifically comprises: S41: Obtain the pre-configured reference dust-free absorbance difference value calibrated in an ideal dust-free and zero-gas environment, input the actual optical path attenuation parameter and the reference dust-free absorbance difference value into the differential calculation module to perform the differential calculation comparison, deduct the background drift caused by the aging of the optical sensor device, and generate a dynamic differential absorption spectrum. S42: Perform peak location processing on the specific gas absorption peak band in the dynamic differential absorption spectrum, and use a polynomial integration algorithm to perform area integration on the envelope below the specific gas absorption peak band to calculate the characteristic absorption area. S43: Obtain the nonlinear concentration inversion mapping curve established in advance through standard concentration methane gas calibration, substitute the characteristic absorption area into the nonlinear concentration inversion mapping curve for interpolation inversion calculation, and output the dynamic methane concentration value. As a further aspect of the present invention, the process of performing the nonlinear spatiotemporal mapping operation specifically includes: The wind speed feature sequence and the original concentration frequency electrical signal are subjected to dynamic time warping to compensate for the asynchronous deviation between the two in the data acquisition frequency. The aligned data sequence is divided into multiple time sliding windows of equal length to generate segmented aligned data windows. Obtain preset phase space reconstruction delay parameters and embedding dimensions, and use the phase space reconstruction delay parameters and embedding dimensions to perform high-dimensional phase space projection transformation on the data stream in the segmented aligned data window to reconstruct the system dynamic evolution characteristics and generate a high-dimensional trajectory matrix; The phase space reconstruction delay parameter is a time interval step determined based on the data sequence mutual information decay rate; The embedding dimension is a spatial dimension value determined based on the requirement of fully unfolding the system's dynamic characteristics. The Lyapunov exponent and mutual information of the high-dimensional trajectory matrix under different spatial node combinations are calculated. The Lyapunov exponent and mutual information are used as indicators of nonlinear coupling strength for normalization processing to generate the spatiotemporal correlation matrix. As a further aspect of the present invention, the process of performing the parameter iteration calculation specifically includes: A continuous computational domain grid containing time and spatial partial derivative terms is established. The initial boundary conditions in the transport mass conservation partial differential equation and the air resistance attenuation coefficient are distributed to each grid node of the continuous computational domain grid to construct a grid parameter distribution model. The convection and diffusion terms in the grid parameter distribution model are discretized and decomposed using an implicit finite difference scheme, transforming the complex continuous differential equations into a set of discretized equations that can be solved by algebraic matrices, thus generating discrete difference algebraic equations. Configure the initial time step and error tolerance threshold, and use the successive over-relaxation iterative algorithm to solve the discrete difference algebraic equation in a loop. After each solution loop, compare the absolute value of the residual between the current time step solution and the previous time step solution. When the absolute value of the residual is less than the error tolerance threshold, terminate the loop and calculate the transient lag time constant. As a further aspect of the present invention, the process of performing the dust turbidity reduction analysis specifically includes: The optical interference characteristics are fitted and compared with the preset Mie scattering theory attenuation model to calculate the equivalent particle size distribution range of suspended dust in the current tunnel environment and the corresponding optical extinction cross-sectional area parameters, and generate multi-band attenuation compensation coefficients. The multi-band attenuation compensation coefficient is used to perform corresponding amplitude amplification and baseline stretching on different frequency components in the waveform leading edge time series to reconstruct the original infrared absorption spectrum profile features and generate a restored spectrum profile model. Extract the ratio of actual transmitted light intensity to incident light intensity at the target characteristic absorption peak of the reduced spectral profile model, perform a negative logarithmic transformation on the ratio of actual transmitted light intensity to incident light intensity to obtain a pure absorbance parameter, and calculate the actual optical path attenuation parameter. As a further aspect of the present invention, the process of performing the differential calculation and comparison specifically includes: Logarithmic scaling transformation is performed on the actual optical path attenuation parameter and the difference in the reference dust-free absorbance, and the logarithmic amplitude envelope of the two in the frequency domain is extracted to generate a smooth absorption envelope. Perform first-order derivative operations with respect to wavelength on the smooth absorption envelope to convert the slowly varying spectral baseline background into tiny fluctuations near zero, highlighting the peak structure of narrowband gas absorption characteristics and generating first-order derivative spectral curves. Calculate the amplitude difference values ​​of the two sets of first derivative spectral curves at corresponding wavelength positions, and filter out broadband absorption interference caused by ambient temperature and light source fluctuations by the extreme value distribution of the difference values ​​to generate the dynamic differential absorption spectrum. A dynamic analysis system for mine gas concentration, the system being used to implement the aforementioned dynamic analysis method for mine gas concentration, the system comprising: The disturbance topology construction module is used to collect airflow fluctuation parameters of the roadway cross section through a mine microwave anemometer, extract the wind speed change vector from the airflow fluctuation parameters, and perform nonlinear spatiotemporal mapping operation on the wind speed change vector and the original concentration frequency electrical signal collected by the gas detection probe to construct an environmental flow field disturbance topology map. The time constant calculation module is used to extract the absolute pressure distribution gradient of discrete nodes in the environmental flow field disturbance topology map, calculate the air flow resistance attenuation coefficient based on the absolute pressure distribution gradient, substitute the air flow resistance attenuation coefficient into the transport mass conservation partial differential equation for parameter iteration calculation, and calculate the transient lag time constant. The optical path attenuation analysis module is used to perform phase reverse compensation shift on the original concentration frequency electrical signal based on the transient lag time constant, generate a waveform leading edge time series, and use the light spot geometric distortion parameter collected by the scattering dust meter to perform dust turbidity restoration analysis on the waveform leading edge time series to calculate the actual optical path attenuation parameter. The concentration value calculation module is used to obtain the pre-configured baseline dust-free absorbance difference value, perform differential calculation and comparison between the actual optical path attenuation parameter and the baseline dust-free absorbance difference value, and output the dynamic gas concentration value based on the differential calculation and comparison result.

[0006] Compared with the prior art, the advantages and positive effects of the present invention are as follows: In this invention, by introducing concentration time-series inversion logic based on spatiotemporal multidimensional airflow coupling characteristics, dynamic airflow feedback parameters output by roadway wind speed transmitters and wind pressure probes are collected in real time. The dynamic airflow feedback parameters are then fitted with nonlinear surfaces using interference fringe displacement parameters generated by optical interferometry detection elements to construct a flow field pressure coupling disturbance matrix. Based on the branch section drag coefficient, time-varying gradient compensation is performed on diffusion transient characteristics to eliminate gas chamber response delay errors caused by pressure fluctuations. Geometric distortion parameters of dust scattering light spots in the detection network are extracted, and an inverse analytical relationship between dust turbidity and optical attenuation index is established. The optical path detection benchmark is dynamically compensated and calibrated, and continuous concentration distribution coordinates reflecting the physical laws of gas emission are output. Attached Figure Description

[0007] Figure 1 This is the overall flow chart of the dynamic analysis method for mine gas concentration of the present invention; Figure 2 Flowchart for constructing an environmental flow field disturbance topology diagram for this invention; Figure 3 This is a flowchart illustrating the calculation of the transient lag time constant in this invention. Figure 4This is a flowchart illustrating the calculation of the actual optical path attenuation parameter in this invention. Figure 5 This invention outputs a flowchart of dynamic gas concentration values. Detailed Implementation

[0008] To make the objectives, technical solutions, and advantages of this invention clearer, the software-based technical solution is described in detail below with reference to system architecture diagrams and embodiments. It should be understood that the specific embodiments described herein are only for explaining the technical solutions of this invention and do not constitute a limitation on the scope of protection.

[0009] In the description of this invention, the system architecture relationships or data processing flows indicated by terms such as "layer," "module," "interface," "data flow," "client," and "server" are all defined based on the architecture diagram or flowchart corresponding to the embodiments. This way of describing is only used to clearly illustrate the logical relationships between the elements in the technical solution, and not to limit the physical deployment form. The term "multiple" includes two or more technical units, including but not limited to multiple data nodes, processing threads, service instances, or functional components and other scalable elements. The specific number is determined according to the actual business scenario and needs to be specifically specified.

[0010] Please see Figure 1 and Figure 2 This invention provides a technical solution: a method for dynamic analysis of mine gas concentration, comprising the following steps: S1: Collect airflow fluctuation parameters of the roadway cross section using a mine microwave anemometer, extract the wind speed change vector from the airflow fluctuation parameters, and perform nonlinear spatiotemporal mapping calculations on the wind speed change vector and the original concentration frequency electrical signal collected by the gas detection probe to construct an environmental flow field disturbance topology map.

[0011] The specific steps of S1 are as follows: S11: Obtain the airflow fluctuation parameters at the roadway cross-section collected by the mine microwave anemometer, perform high-frequency filtering and noise reduction on the airflow fluctuation parameters, filter out the periodic mechanical vibration interference noise caused by the mine ventilation fan, and extract the three-dimensional spatial direction wind speed change vector reflecting the local wind field turbulence characteristics to generate a wind speed feature sequence.

[0012] S12: Acquire the original concentration frequency electrical signal collected by the gas detection probe, use the wind speed feature sequence and the original concentration frequency electrical signal as multivariate time series inputs, perform nonlinear spatiotemporal mapping operation on the multivariate time series inputs, calculate the cross-information entropy between different physical measuring points at different time delay scales, and generate a spatiotemporal correlation matrix.

[0013] The process of performing nonlinear spacetime mapping operations specifically includes: Dynamic time warping is performed on the wind speed feature sequence and the original concentration frequency electrical signal to compensate for the asynchronous deviation between the two in the data acquisition frequency. The aligned data sequence is divided into multiple time sliding windows of equal length to generate segmented aligned data windows.

[0014] Obtain the preset phase space reconstruction delay parameters and embedding dimension, and use the phase space reconstruction delay parameters and embedding dimension to perform high-dimensional phase space projection transformation on the data stream in the segmented aligned data window to reconstruct the dynamic evolution characteristics of the system and generate a high-dimensional trajectory matrix.

[0015] The phase space reconstruction delay parameter is the time interval step determined based on the mutual information decay rate of the data sequence.

[0016] The embedding dimension is a spatial dimension value determined based on the requirements for fully unfolding the system's dynamic characteristics.

[0017] The Lyapunov exponent and mutual information of the high-dimensional trajectory matrix under different spatial node combinations are calculated. The Lyapunov exponent and mutual information are used as indicators of nonlinear coupling strength for normalization processing to generate a spatiotemporal correlation matrix.

[0018] S13: Extract matrix elements in the spatiotemporal correlation matrix whose values ​​exceed the preset structure retention threshold as weight values ​​of directed edges, map each physical monitoring point in the well to an independent connected node in the network, and establish a spatial connectivity network between independent connected nodes based on the weight values ​​of directed edges to construct an environmental flow field disturbance topology map.

[0019] A sequence of 1000 discrete wind speed measurement points was continuously collected by a mine microwave anemometer installed 3 meters above the tunnel cross-section at a sampling frequency of 100 Hz, generating a wind flow fluctuation parameter sequence. This sequence was then subjected to high-frequency filtering and noise reduction. The sliding data window length was set to 5 consecutive sampling points. The wind speed measurement data within the window were summed by arithmetic and then divided to output a smoothed basic wind flow data sequence. This basic wind flow data sequence was input into a band-stop digital filter with a center stopband frequency of 50 Hz. Frequency domain signal blocking was performed to remove the mechanical vibration frequency band corresponding to 50 Hz. The longitudinal wind speed variation components along the tunnel axis, the transverse wind speed variation components along the tunnel cross-section, and the vertical wind speed variation components perpendicular to the tunnel floor were extracted. These three spatial dimensions of wind speed variation components were then matrix-concatenated according to the chronological order of the collection timestamps to generate a wind speed feature sequence matrix containing three columns of data. For example, a microwave anemometer in a mine can collect data showing a longitudinal wind speed change of 2.5 m / s, a lateral wind speed change of 0.2 m / s, and a vertical wind speed change of 0.1 m / s at a given moment. These three values ​​can be directly combined to form a 3D wind speed feature vector corresponding to that moment. The raw concentration frequency electrical signal is obtained by acquiring gas detection probes placed on the sidewall of the roadway at an equivalent sampling frequency of 100 Hz. The wind speed feature sequence matrix and the original concentration frequency electrical signal are used as inputs for a multivariate time series and dynamic time warping is performed. The operation logic is to calculate the Euclidean geometric absolute distance between each discrete data vector in the wind speed feature sequence matrix and each discrete data point in the original concentration frequency electrical signal in turn. All the calculated Euclidean geometric absolute distances are filled into a two-dimensional matrix grid to construct a distance penalty matrix. A dynamic programming search algorithm is used to find the connected node path with the minimum cumulative sum of values ​​from the starting coordinate of the upper left corner to the ending coordinate of the lower right corner in the distance penalty matrix. Based on the connected node path with the minimum cumulative sum of values, the time axis index of the wind speed feature sequence matrix and the time axis index of the original concentration frequency electrical signal are remapped and bound. The aligned and recombined data sequence is forcibly divided into multiple time sliding windows of equal length according to the time span of 5 seconds, generating segmented aligned data windows. The phase space reconstruction delay parameter and embedding dimension value are obtained. The process of setting the phase space reconstruction delay parameter is to calculate the mutual information decay value of the data sequence in the segmented aligned data window at different time intervals in sequence. The time interval step value corresponding to the first drop of the mutual information decay value to below 10% of the mutual information value in the initial zero delay state is found in ascending order of time step size. The time interval step value is directly assigned as the phase space reconstruction delay parameter.For example, the mutual information value in the initial zero-delay state is recorded as 5.0. As the step size increases, the mutual information value is calculated. When the calculation reaches the 8th time step, the mutual information value decays to 0.45. This value is lower than the 0.5 threshold obtained by multiplying the initial value of 5.0 by 10%. Therefore, the phase space reconstruction delay parameter is directly set to 8. The process of setting the embedding dimension uses the false nearest neighbor calculation logic. It calculates the Euclidean distance between two adjacent data points in the current dimension. Then, it adds one dimension and recalculates the Euclidean distance between these two data points in the new dimension. The absolute difference between the two Euclidean distances before and after adding the dimension is obtained. This absolute difference is divided by the Euclidean distance in the original dimension to calculate the spatial dispersion ratio. Data points with a spatial dispersion ratio greater than the judgment benchmark threshold of 15 are marked as false nearest neighbors. The percentage of false nearest neighbors among all data points is counted. When this percentage drops to less than 5%, the current test spatial dimension value is determined as the final embedding dimension. For example, when testing the third dimension, the percentage of false nearest neighbors is 12%. When the fourth dimension is added, this percentage decreases to 3%, meeting the condition of less than 5%, so the embedding dimension is set to 4. Using the phase space reconstruction delay parameter value of 8 and the embedding dimension value of 4, a high-dimensional phase space projection transformation is performed on the data stream in the segmented aligned data window. The original one-dimensional time series is truncated into data points at intervals of 8 time steps, and vector combinations are performed in a 4-dimensional spatial coordinate system to generate a high-dimensional trajectory matrix. The Lyapunov exponent of mutual evolution between spatially discrete nodes in the high-dimensional trajectory matrix and the mutual information value between node states are calculated. The calculated Lyapunov exponent and mutual information value are substituted into a linear weighted summation operation logic for normalization. This linear weighted summation operation logic multiplies the Lyapunov exponent by the first weight coefficient to obtain the dynamic evolution parameter, multiplies the mutual information value by the second weight coefficient to obtain the static correlation parameter, and performs an arithmetic addition operation on the dynamic evolution parameter and the static correlation parameter to generate the spatiotemporal correlation value. The first and second weighting coefficients were set using an arithmetic progression test based on previous downhole measured data calibration sets. The first weighting coefficient was set to 0.6 and the second weighting coefficient to 0.4. For example, when the Lyapunov exponent of a combination of two spatial nodes in the high-dimensional trajectory matrix is ​​calculated to be 0.35 and the mutual information value is 1.2, these values ​​are substituted into the calculation logic. 0.35 multiplied by 0.6 yields the dynamic evolution parameter 0.21, and 1.2 multiplied by 0.4 yields the static correlation parameter 0.48. Adding 0.21 and 0.48 gives the spatiotemporal correlation degree of the node combination, which is 0.69. The spatiotemporal correlation degree values ​​of all node combinations are calculated and arranged to generate a multi-row, multi-column spatiotemporal correlation degree matrix. Matrix elements in the spatiotemporal correlation degree matrix whose values ​​exceed a preset structure retention threshold are extracted as the weight values ​​of directed edges. The process of setting the structure retention threshold is determined using the experimental data listed in Table 1 below.

[0020] Table 1. Experimental parameters for candidate structure retention threshold

[0021] Table 1 lists the experimental data of network topology attributes under different candidate structure retention thresholds. Based on the connectivity constraint standard established by the network structure, the number of isolated unconnected nodes in the connected network must be 0, and the overall edge density percentage of the network must be within the range of 15% to 25%. The candidate structure retention threshold that meets this constraint standard is determined as the optimal threshold. As shown in Table 1, when the candidate structure retention threshold is set to 0.60, the number of isolated unconnected nodes is 0, and the overall edge density percentage of the network is 18.6%, which meets the above conditions. Therefore, the structure retention threshold is set to 0.60. All matrix elements with values ​​greater than 0.60 in the spatiotemporal correlation matrix are extracted, and their corresponding spatiotemporal correlation values ​​are directly recorded as the weight values ​​of directed edges. Each physical monitoring point actually distributed underground is mapped into an independent connected node in the network model according to its coordinate position. Based on the extracted weight values ​​of directed edges, connecting line segments with numerical labels are drawn between the corresponding independent connected nodes to construct an environmental flow field disturbance topology map containing node positions and edge weights.

[0022] The aforementioned band-stop digital filter refers to an electronic filtering device or algorithm model that allows signals of most frequencies to pass through, but attenuates signals of certain specific frequency bands to an extremely low level, in order to eliminate interference noise at specific frequencies.

[0023] Please see Figure 1 and Figure 3 S2: Extract the absolute pressure distribution gradient of discrete nodes in the topology of the environmental flow field disturbance, calculate the air flow resistance attenuation coefficient based on the absolute pressure distribution gradient, substitute the air flow resistance attenuation coefficient into the partial differential equation of transport mass conservation for parameter iteration calculation, and calculate the transient lag time constant.

[0024] The specific steps of S2 are as follows: S21: Extract the absolute pressure monitoring data and three-dimensional spatial coordinate information of each discrete node in the environmental flow field disturbance topology map, calculate the ratio of the absolute pressure difference between adjacent discrete nodes to the Euclidean distance in physical space, and generate the absolute pressure distribution gradient.

[0025] S22: Obtain the pre-determined frictional resistance coefficient of the roadway wall and the fluid dynamic viscosity parameters, calculate the degree of local energy loss of gas during the flow process based on the absolute pressure distribution gradient, and calculate the air flow resistance attenuation coefficient in combination with the equivalent hydraulic diameter of the roadway.

[0026] S23: Obtain the partial differential equation of transport mass conservation that reflects the comprehensive law of gas molecule diffusion and convection transport. Substitute the air flow resistance attenuation coefficient as the updated damping variable into the partial differential equation of transport mass conservation to perform parameter iterative calculation until the calculated rate of change of time response parameters reaches a stable convergence state. Calculate the transient lag time constant.

[0027] The process of performing parameter iteration calculations specifically includes: A continuous computational domain grid containing time and spatial partial derivative terms is established. The initial boundary conditions and air resistance attenuation coefficients in the transport mass conservation partial differential equation are distributed to each grid node of the continuous computational domain grid to construct a grid parameter distribution model.

[0028] An implicit finite difference scheme is used to discretize and decompose the convection and diffusion terms in the grid parameter distribution model, transforming the complex continuous differential equations into a set of discretized equations that can be solved by algebraic matrices, thus generating discrete difference algebraic equations.

[0029] Configure the initial time step and error tolerance threshold, and use the successive over-relaxation iterative algorithm to solve the discrete difference algebraic equation in a loop. After each solution loop, compare the absolute value of the residual between the current time step solution and the previous time step solution. When the absolute value of the residual is less than the error tolerance threshold, terminate the loop and calculate the transient lag time constant.

[0030] The absolute pressure monitoring data of each discrete node contained within the environmental flow field disturbance topology map and the 3D spatial coordinate information determined based on mine engineering drawings are extracted. The absolute pressure difference between two adjacent discrete nodes is calculated, and the physical spatial Euclidean distance between the 3D spatial coordinates of these two discrete nodes is calculated. The absolute pressure difference is divided by the physical spatial Euclidean distance to calculate the absolute pressure distribution gradient. For example, if the absolute pressure monitoring data of discrete node A is 101325 Pascals and the absolute pressure monitoring data of discrete node B is 101300 Pascals, the absolute pressure difference between the two is calculated to be 25 Pascals. Using the coordinates, the physical spatial Euclidean distance between node A and node B is calculated to be 50 meters. Dividing 25 Pascals by 50 meters yields an absolute pressure distribution gradient of 0.5 Pascals per meter between node A and node B. The system obtains the roadway wall friction resistance coefficient, pre-determined through actual measurements using an anemometer and barometer in the corresponding roadway section; the fluid dynamic viscosity parameter obtained from a standard fluid property table at the current ambient temperature; and the equivalent hydraulic diameter of the roadway calculated based on the relationship between the roadway cross-section perimeter and area. Based on these parameters, the system executes the calculation logic for the air resistance attenuation coefficient. This logic involves multiplying the roadway wall friction resistance coefficient by the previously calculated absolute pressure distribution gradient to obtain the basic resistance parameter; dividing the fluid dynamic viscosity parameter by the equivalent hydraulic diameter of the roadway to calculate the viscous loss parameter; combining the basic resistance parameter and the viscous loss parameter to obtain the comprehensive resistance value; and dividing the comprehensive resistance value by the square root of the standard air density to calculate the air resistance attenuation coefficient. For example, the obtained tunnel wall friction resistance coefficient is 0.015, the obtained absolute pressure distribution gradient is 0.5 Pascals per meter, the obtained fluid dynamic viscosity parameter is 0.000018 Pascals per second, the obtained equivalent hydraulic diameter of the tunnel is 3.0 meters, and the standard air density is taken as 1.2 kg per cubic meter. Multiplying 0.015 by 0.5 Pascals per meter yields the basic resistance parameter of 0.0075. Dividing 0.000018 Pascals per second by 3.0 meters yields the viscous loss parameter of 0.000006. Adding 0.0075 and 0.000006 yields the comprehensive resistance value of 0.007506. The square root of the standard air density of 1.2 is approximately 1.0954. Dividing 0.007506 by 1.0954 finally yields an air flow resistance attenuation coefficient of approximately 0.00685.A transport mass conservation partial differential equation characterizing the relationship between methane concentration and spatial coordinates and time was obtained. A continuous computational domain grid covering the three-dimensional space of the tunnel was established. The continuous computational domain grid was divided into multiple hexahedral grid cells with a grid spacing of 0.5 meters. The initial concentration value and boundary wind speed input value in the transport mass conservation partial differential equation were used as initial boundary conditions. Together with the air flow resistance attenuation coefficient of 0.00685 obtained above, they were assigned to the parameter matrix of each hexahedral grid node in the continuous computational domain grid to construct a grid parameter distribution model. An implicit finite difference scheme is used to discretize the grid parameter distribution model. For the time partial derivative terms in the transport mass conservation partial differential equation, a forward difference replacement is performed by dividing the difference in grid node concentration between the current time step and the next time step by the time step size. For the spatial partial derivatives contained in the spatial convection and diffusion terms in the transport mass conservation partial differential equation, a central difference replacement is performed by dividing the difference in concentration between the two adjacent spatial nodes of the current grid node by twice the grid spacing. After the above difference replacement operations, the continuous partial differential equations are transformed into a system of linear algebraic equations containing multiple unknown node concentrations, generating a discrete difference algebraic equation matrix. The initial time step is configured to be 0.1 seconds, and an adaptive time step adjustment mechanism is introduced to ensure computational stability. The error tolerance threshold is configured to be 0.001. A successive over-relaxation iterative algorithm is used to solve the discrete difference algebraic equation matrix iteratively. The execution logic of the successive over-relaxation iterative algorithm is as follows: a relaxation factor term is added to the right side of the conventional Gauss-Seidel iterative formula, with the relaxation factor value set to 1.2. In each iteration, the predicted solution for the current node is calculated. The difference between this predicted solution and the retained solution for the same node in the previous iteration is calculated. This difference is multiplied by the relaxation factor value of 1.2 and accumulated above the retained solution to generate the updated solution for the current time step. After each solution cycle, the difference between the current time step solution and the previous time step solution for each grid node is calculated one by one, and the absolute value is taken. The maximum value among all the absolute differences is taken as the current absolute value of the residual. Table 2 below records the data of the iteration process.

[0031] Table 2. Residual Convergence Table for Successive Over-Relaxation Iterations

[0032] Table 2 lists the residual convergence trend of the successive over-relaxation iterative algorithm during the solution process. The system continuously compares the absolute value of the residual with the error tolerance threshold of 0.001. Referring to Table 2, in the 12th iteration, the absolute value of the residual is 0.0050, which is greater than 0.001, and the system instructs the iteration process to continue. When the 15th iteration is executed, the absolute value of the residual drops to 0.0008, which is less than the error tolerance threshold of 0.001, and the system determines that the convergence condition has been met and terminates the iterative solution operation. The cumulative total value of the time step calculated based on the adaptive step size when the convergence state is achieved is extracted as the transient lag time constant. For example, based on the above iteration results, the final determined transient lag time constant is 2.65 seconds.

[0033] The aforementioned implicit finite difference scheme refers to a discretization method that uses the values ​​of nodal variables at unknown times to construct difference equations in numerical solutions of partial differential equations, thereby improving the stability of numerical computation.

[0034] Please see Figure 1 and Figure 4 S3: Based on the transient lag time constant, the original concentration frequency electrical signal is phase-reverse compensated and shifted to generate a waveform leading edge time series. The dust turbidity of the waveform leading edge time series is restored and analyzed using the light spot geometric distortion parameter collected by the scattering dust meter, and the actual optical path attenuation parameter is calculated.

[0035] The specific steps for S3 are as follows: S31: Obtain the transient lag time constant, determine the physical time delay span required for the gas to diffuse from the release source to the optical detection path based on the transient lag time constant, perform phase inverse compensation translation on the original concentration frequency electrical signal based on the physical time delay span, filter out the phase deviation caused by the transport delay, and generate the waveform leading edge time series.

[0036] S32: Obtain the geometric distortion parameters of the light spot on the detection surface of the optical receiver collected by the scattering dust meter, analyze the spatial scattering and blocking effect of suspended particles on the detection beam characterized by the geometric distortion parameters of the light spot, and extract optical interference features.

[0037] S33: Utilize optical interference characteristics to perform baseline compensation and dust turbidity restoration analysis on the abnormal amplitude depression region caused by dust scattering in the waveform leading edge time series, deduct non-absorbent optical loss error, and calculate the actual optical path attenuation parameter.

[0038] The specific steps involved in performing dust turbidity reduction analysis include: By fitting and comparing the optical interference characteristics with the preset Mie scattering theory attenuation model, the equivalent particle size distribution range of suspended dust in the current tunnel environment and the corresponding optical extinction cross-sectional area parameters are calculated, and multi-band attenuation compensation coefficients are generated.

[0039] By using multi-band attenuation compensation coefficients, the amplitude of different frequency components in the waveform leading edge time series is amplified and the baseline is stretched accordingly, the original infrared absorption spectrum profile features are reconstructed, and the restored spectrum profile model is generated.

[0040] The ratio of actual transmitted light intensity to incident light intensity at the target characteristic absorption peak is extracted and reduced from the spectral profile model. A negative logarithmic transformation is performed on the ratio of actual transmitted light intensity to incident light intensity to obtain the pure absorbance parameter, and the actual optical path attenuation parameter is calculated.

[0041] The transient lag time constant value of 2.65 seconds, calculated in step 2, is obtained. Based on this transient lag time constant, the physical time delay span consumed by the gas as it diffuses from the source of the roadway and is carried by the airflow to the detection optical path inside the optical sensor is determined to be 2.65 seconds. A phase-reverse compensation translation operation is performed on the original concentration frequency electrical signal on the time axis based on this physical time delay span value. Specifically, combining the previously set signal sampling frequency of 100 Hz, the number of discrete sampling point offsets corresponding to the 2.65-second physical time delay span is calculated. That is, 2.65 seconds is multiplied by 100 Hz to obtain 265 data points. The original concentration frequency electrical signal is then shifted 265 data points towards the starting point of the time, generating a waveform leading-edge time series. The geometric distortion parameter data of the light spot transmitted by the photosensitive detection panel array of the scattering dust meter arranged on the same roadway cross-section is obtained. This geometric distortion parameter data includes the eccentricity value of the light spot edge and the diffusion rate value of the light spot area. The eccentricity values ​​at the edge of the light spot and the diffuser values ​​of the light spot area are extracted as optical interference feature vectors characterizing the spatial scattering and blocking effect of environmental suspended particulate matter on the detection beam. The Mie scattering theory attenuation model pre-stored in the environmental dust monitoring system is obtained. The aforementioned optical interference feature vectors are input into the fitting equation of the Mie scattering theory attenuation model. The unknown parameters of the fitting equation are iteratively approximated using the Gauss-Newton nonlinear least squares method to calculate the median equivalent particle size distribution of suspended dust in the current tunnel air and the corresponding optical extinction cross-sectional area parameter. For example, the fitting yields an equivalent particle size distribution median of 2.5 micrometers, with a corresponding optical extinction cross-sectional area parameter of 0.00035 square micrometers. The calculated optical extinction cross-sectional area parameter of 0.00035 square micrometers is substituted into the derivation formula of the Lambert-Beer scattering law for multiplication and conversion to calculate the light intensity attenuation ratio at different frequency bands, generating a corresponding multi-band attenuation compensation coefficient array. The multi-band attenuation compensation coefficient array is used to perform frequency domain restoration on the waveform leading-edge time series data. For the abnormal amplitude dips caused by dust scattering in the waveform leading-edge time series, the amplitude of each frequency component of the waveform leading-edge time series is directly multiplied by the compensation coefficient value of the corresponding frequency band in the multi-band attenuation compensation coefficient array, performing amplitude amplification. Zero-mean baseline stretching is then performed in the non-characteristic absorption bands to reconstruct and output the restored spectral profile model. The actual transmitted light intensity at the target characteristic absorption peak of 3.31 micrometers and the incident light intensity under gas-free absorption conditions are extracted from the restored spectral profile model. The ratio between the actual transmitted light intensity and the incident light intensity is calculated. For example, if the actual transmitted light intensity at the target characteristic absorption peak is 0.8 milliwatts and the incident light intensity is 1.0 milliwatts, the ratio 0.8 milliwatts divided by 1.0 milliwatts is 0.8.Perform a negative logarithmic transformation operation with the natural constant as the base on the ratio 0.8 to calculate the negative natural logarithm function value. For example, the natural logarithm of the ratio 0.8 is approximately negative 0.223. After taking the negative sign, we get 0.223. Use this value of 0.223 as the actual optical path attenuation parameter for extracting the target absorbance characteristics.

[0042] The aforementioned Mie scattering theory attenuation model refers to an analytical calculation model based on Maxwell's electromagnetic wave theory to solve the scattering effect of spherical particles on plane light waves, which is used to quantitatively evaluate the blocking effect of dust particles on the detection beam.

[0043] Please see Figure 1 and Figure 5 S4: Obtain the pre-configured baseline dust-free absorbance difference value, perform differential calculation and comparison between the actual optical path attenuation parameter and the baseline dust-free absorbance difference value, and output the dynamic gas concentration value based on the differential calculation and comparison results.

[0044] The specific steps for S4 are as follows: S41: Obtain the pre-configured baseline dust-free absorbance difference obtained under ideal dust-free and zero-gas environment, input the actual optical path attenuation parameter and the baseline dust-free absorbance difference into the differential calculation module to perform differential calculation and comparison, subtract the background drift caused by the aging of optical sensor devices, and generate a dynamic differential absorption spectrum.

[0045] The specific process of performing differential calculation and comparison includes: Logarithmic scaling transformations are performed on the actual optical path attenuation parameter and the difference in the baseline dust-free absorbance, and the logarithmic amplitude envelope of the two in the frequency domain is extracted to generate a smooth absorption envelope.

[0046] By performing first-order derivative operations with respect to wavelength on the smooth absorption envelope, the slowly varying spectral baseline background is transformed into tiny fluctuations near zero, highlighting the peak structure of narrowband gas absorption characteristics and generating first-order derivative spectral curves.

[0047] Calculate the amplitude difference values ​​of the two sets of first derivative spectral curves at corresponding wavelength positions, and filter out broadband absorption interference caused by ambient temperature and light source fluctuations by analyzing the extreme value distribution of the difference values ​​to generate a dynamic differential absorption spectrum.

[0048] S42: Perform peak location processing on specific gas absorption peak bands in the dynamic differential absorption spectrum, and use a polynomial integration algorithm to perform area integration on the envelope surface below the specific gas absorption peak band to calculate the characteristic absorption area.

[0049] S43: Obtain the nonlinear concentration inversion mapping curve established in advance through standard concentration methane gas calibration, substitute the characteristic absorption area into the nonlinear concentration inversion mapping curve for interpolation inversion calculation, and output the dynamic methane concentration value.

[0050] Obtain the difference in absorbance of the optical sensor under no-load conditions by continuously testing for 1 hour in an ideal, dust-free, zero-gas environment constructed in a pre-constructed, temperature-controlled, sealed gas calibration chamber filled with high-purity nitrogen. For example, obtain a baseline absorbance difference of 0.005. The actual optical path attenuation parameter 0.223 calculated in step 3 and the difference between the reference dust-free absorbance and 0.005 are input into the difference calculation node. A scaling transformation based on a common logarithmic function is performed on these two parameters to generate corresponding logarithmic data point sets. An envelope detection algorithm is used to connect the local maxima in the logarithmic data point sets, extracting the logarithmic amplitude envelope in the frequency domain. First-order derivative operations with respect to the scanning wavelength are performed on the logarithmic amplitude envelope corresponding to the actual optical path attenuation parameter and the logarithmic amplitude envelope corresponding to the reference dust-free absorbance difference, respectively. The numerical sequence of the tangent slope of the logarithmic amplitude envelope at each wavelength coordinate point is extracted, generating a set of first-order derivative spectral curves containing a peak structure. The amplitudes of the first-order derivative spectral curve of the actual optical path attenuation parameter and the first-order derivative spectral curve of the reference dust-free absorbance difference are subtracted at the same wavelength position to obtain the difference numerical sequence. The extreme value distribution coordinates in the differential numerical sequence are searched, and broadband absorption interference data segments with extreme value widths exceeding 50 nm wavelength ranges are removed, retaining narrowband data sequences with extreme value widths within 20 nm to generate a dynamic differential absorption spectrum. In the dynamic differential absorption spectrum, the first derivative values ​​are scanned point by point along the wavelength increasing direction to find the wavelength coordinate position corresponding to the point where the first derivative value changes from positive to negative and crosses zero. This wavelength coordinate position is locked as the center point of the specific gas absorption peak. The wavelength span of 10 nm on both sides of the center point is determined as the specific gas absorption peak band. The Simpson polynomial integration algorithm is used to perform area integration on the envelope contour below the specific gas absorption peak band. The specific execution logic of the Simpson polynomial integral algorithm is as follows: the integration band interval from 10 nm to -10 nm is divided into 10 sub-intervals, resulting in 11 discrete wavelength nodes. A weight of 1 is assigned to the amplitude data at the first and last two wavelength nodes, a weight of 4 is assigned to the amplitude data at odd-numbered wavelength nodes, and a weight of 2 is assigned to the amplitude data at even-numbered wavelength nodes. All weighted amplitude data are then summed. The sum is multiplied by the sub-interval step size and divided by 3 to calculate the characteristic absorption area. For example, after various weighted summations and multiplication / division operations, the calculated characteristic absorption area is 0.045. A parameter model for the nonlinear concentration inversion mapping curve, pre-established through calibration with multiple sets of standard concentration methane gas, is obtained. Table 3 below shows the correlation of the basic experimental data for this inversion calibration.

[0051] Table 3. Standard Inversion Calibration Table for Gas Concentration

[0052] Table 3 shows the data from multiple discrete test points obtained through the gradient concentration ventilation calibration experiment. The system uses the integral value of the characteristic absorption area in Table 3 as the input variable and the percentage of standard gas concentration as the output variable. Based on the least squares principle, it fits and generates a quadratic polynomial nonlinear concentration inversion mapping curve equation. The characteristic absorption area value of 0.045, calculated by the integral in the previous step, is directly substituted into the independent variable position of this nonlinear concentration inversion mapping curve equation to perform mathematical multiplication and addition interpolation inversion operations, and the corresponding result is calculated based on the polynomial coefficients. Combining the data distribution range in Table 3, 0.045 falls between 0.038 and 0.055. Substituting this value into the curve equation yields a corresponding dependent variable value of approximately 1.20%. This value of 1.20% is then output as the final dynamic gas concentration value. The experimental results show that the absolute deviation between the numerical result of 1.20% of this technical solution and the actual concentration of 1.21% obtained by the on-site gas sampling and testing by the portable chromatograph carried by the mine gas inspector is small. Compared with the traditional infrared direct reading method without wind speed and air pressure compensation, the deviation is reduced by more than 80% with a static error of up to 0.15%.

[0053] The aforementioned Simpson polynomial integral algorithm refers to a numerical analysis method that approximates the integrand curve piecewise using a quadratic polynomial and then approximates the definite integral value by performing area-weighted summation on the regions below each parabola segment.

[0054] A dynamic analysis system for mine gas concentration, used to execute the aforementioned dynamic analysis method for mine gas concentration, the system comprising: The disturbance topology construction module is used to collect airflow fluctuation parameters of the roadway cross section through a mine microwave anemometer, extract the wind speed change vector in the airflow fluctuation parameters, and perform nonlinear spatiotemporal mapping operation on the wind speed change vector and the original concentration frequency electrical signal collected by the gas detection probe to construct an environmental flow field disturbance topology map. The time constant calculation module is used to extract the absolute pressure distribution gradient of discrete nodes in the topology of the environmental flow field disturbance, calculate the air flow resistance attenuation coefficient based on the absolute pressure distribution gradient, substitute the air flow resistance attenuation coefficient into the partial differential equation of transport mass conservation for parameter iteration calculation, and calculate the transient lag time constant. The optical path attenuation analysis module is used to perform phase reverse compensation shift on the original concentration frequency electrical signal based on the transient lag time constant, generate the waveform leading edge time series, and use the light spot geometric distortion parameter collected by the scattering dust meter to perform dust turbidity restoration analysis on the waveform leading edge time series to calculate the actual optical path attenuation parameter. The concentration value calculation module is used to obtain the pre-configured baseline dust-free absorbance difference value, perform differential calculation and comparison between the actual optical path attenuation parameter and the baseline dust-free absorbance difference value, and output the dynamic gas concentration value based on the differential calculation and comparison results.

[0055] The above embodiments illustrate preferred embodiments of the present invention. Any equivalent adjustments to the technical solution based on software engineering methods are within the scope of protection, including but not limited to: implementing algorithm logic using different programming languages, refactoring functional modules into services, adjusting data interaction protocols, and optimizing resource scheduling strategies. Any implementation scheme derived from reasonable modifications to the data processing flow, service call chain, or system architecture layer without departing from the core technology of the present invention should be considered within the protection scope defined by the technical solution of the present invention.

Claims

1. A method for dynamic analysis of mine gas concentration, characterized in that, Includes the following steps: S1: Collect airflow fluctuation parameters of the roadway cross section using a mine microwave anemometer, extract the wind speed change vector from the airflow fluctuation parameters, and perform nonlinear spatiotemporal mapping operation on the wind speed change vector and the original concentration frequency electrical signal collected by the gas detection probe to construct an environmental flow field disturbance topology map. S2: Extract the absolute pressure distribution gradient of discrete nodes in the environmental flow field disturbance topology map, calculate the air flow resistance attenuation coefficient based on the absolute pressure distribution gradient, substitute the air flow resistance attenuation coefficient into the transport mass conservation partial differential equation for parameter iteration calculation, and calculate the transient lag time constant. S3: Based on the transient lag time constant, the original concentration frequency electrical signal is phase-reverse compensated and shifted to generate a waveform leading edge time series. The dust turbidity of the waveform leading edge time series is restored and analyzed using the light spot geometric distortion parameter collected by the scattering dust meter, and the actual optical path attenuation parameter is calculated. S4: Obtain the pre-configured baseline dust-free absorbance difference value, perform differential calculation and comparison between the actual optical path attenuation parameter and the baseline dust-free absorbance difference value, and output the dynamic gas concentration value based on the differential calculation and comparison result. The specific steps of S3 are as follows: S31: Obtain the transient lag time constant, determine the physical time delay span required for the gas to diffuse from the release source to the optical detection path based on the transient lag time constant, perform the phase reverse compensation shift on the original concentration frequency electrical signal based on the physical time delay span, filter out the phase deviation caused by the shift delay, and generate the waveform leading edge time series; S32: Obtain the geometric distortion parameters of the light spot on the detection surface of the optical receiver collected by the scattering dust meter, analyze the spatial scattering and blocking effect of suspended particulate matter on the detection beam characterized by the geometric distortion parameters of the light spot, and extract optical interference features; S33: Using the optical interference characteristics, perform baseline compensation and dust turbidity reduction analysis on the abnormal amplitude depression region caused by dust scattering in the waveform leading edge time series, deduct non-absorbent optical loss error, and calculate the actual optical path attenuation parameter. The specific steps of S4 are as follows: S41: Obtain the pre-configured reference dust-free absorbance difference value calibrated in an ideal dust-free and zero-gas environment, input the actual optical path attenuation parameter and the reference dust-free absorbance difference value into the differential calculation module to perform the differential calculation comparison, deduct the background drift caused by the aging of the optical sensor device, and generate a dynamic differential absorption spectrum. S42: Perform peak location processing on the gas absorption peak band in the dynamic differential absorption spectrum, and use a polynomial integration algorithm to perform area integration on the envelope below the gas absorption peak band to calculate the characteristic absorption area. S43: Obtain the nonlinear concentration inversion mapping curve established in advance through standard concentration methane gas calibration, substitute the characteristic absorption area into the nonlinear concentration inversion mapping curve for interpolation inversion calculation, and output the dynamic methane concentration value.

2. The mine gas concentration dynamic analysis method according to claim 1, characterized in that, The specific steps of S1 are as follows: S11: Obtain the airflow fluctuation parameters at the roadway cross-section collected by the mine microwave anemometer, perform high-frequency filtering and noise reduction processing on the airflow fluctuation parameters to filter out the periodic mechanical vibration interference noise caused by the mine ventilation fan, and extract the wind speed change vector in the three-dimensional spatial direction that reflects the local wind field turbulence characteristics to generate a wind speed feature sequence. S12: Acquire the original concentration frequency electrical signal collected by the gas detection probe, use the wind speed feature sequence and the original concentration frequency electrical signal as multivariate time series input items, perform the nonlinear spatiotemporal mapping operation on the multivariate time series input items, calculate the cross-information entropy between different physical measuring points at different time delay scales, and generate a spatiotemporal correlation matrix. S13: Extract the matrix elements in the spatiotemporal correlation matrix whose values ​​exceed the preset structure retention threshold as the weight values ​​of the directed edges, map each physical monitoring point in the well to an independent connected node in the network, and establish a spatial connectivity network between the independent connected nodes based on the weight values ​​of the directed edges to construct the environmental flow field disturbance topology map.

3. The mine gas concentration dynamic analysis method according to claim 1, characterized in that, The specific steps of S2 are as follows: S21: Extract the absolute pressure monitoring data and three-dimensional spatial coordinate information of each discrete node in the environmental flow field disturbance topology map, calculate the ratio of the absolute pressure difference between adjacent discrete nodes to the Euclidean distance in physical space, and generate the absolute pressure distribution gradient. S22: Obtain the pre-determined frictional resistance coefficient of the tunnel wall and the fluid dynamic viscosity parameters, calculate the degree of local energy loss of the gas during the flow process based on the absolute pressure distribution gradient, and calculate the air flow resistance attenuation coefficient in combination with the equivalent hydraulic diameter of the tunnel. S23: Obtain the partial differential equation of transport mass conservation that reflects the comprehensive law of gas molecule diffusion and convection transport, substitute the air flow resistance attenuation coefficient as the updated damping variable into the partial differential equation of transport mass conservation, and perform the parameter iteration operation until the calculated rate of change of the time response parameter reaches a stable convergence state, and calculate the transient lag time constant.

4. The method for dynamic analysis of mine gas concentration according to claim 2, characterized in that, The process of performing the nonlinear spatiotemporal mapping operation specifically includes: The wind speed feature sequence and the original concentration frequency electrical signal are subjected to dynamic time warping to compensate for the asynchronous deviation between the two in the data acquisition frequency. The aligned data sequence is divided into multiple time sliding windows of equal length to generate segmented aligned data windows. Obtain preset phase space reconstruction delay parameters and embedding dimensions, and use the phase space reconstruction delay parameters and embedding dimensions to perform high-dimensional phase space projection transformation on the data stream in the segmented aligned data window to reconstruct the system dynamic evolution characteristics and generate a high-dimensional trajectory matrix; The phase space reconstruction delay parameter is a time interval step determined based on the data sequence mutual information decay rate; The embedding dimension is a spatial dimension value determined based on the requirement of fully unfolding the system's dynamic characteristics. The Lyapunov exponent and mutual information of the high-dimensional trajectory matrix under different spatial node combinations are calculated. The Lyapunov exponent and mutual information are used as indicators of nonlinear coupling strength for normalization processing to generate the spatiotemporal correlation matrix.

5. The method for dynamic analysis of mine gas concentration according to claim 3, characterized in that, The process of performing the parameter iteration calculation specifically includes: A continuous computational domain grid containing time and spatial partial derivative terms is established. The initial boundary conditions in the transport mass conservation partial differential equation and the air resistance attenuation coefficient are distributed to each grid node of the continuous computational domain grid to construct a grid parameter distribution model. The convection and diffusion terms in the grid parameter distribution model are discretized and decomposed using an implicit finite difference scheme, transforming the complex continuous differential equations into a set of discretized equations that can be solved by algebraic matrices, thus generating discrete difference algebraic equations. Configure the initial time step and error tolerance threshold, and use the successive over-relaxation iterative algorithm to solve the discrete difference algebraic equation in a loop. After each solution loop, compare the absolute value of the residual between the current time step solution and the previous time step solution. When the absolute value of the residual is less than the error tolerance threshold, terminate the loop and calculate the transient lag time constant.

6. The method for dynamic analysis of mine gas concentration according to claim 1, characterized in that, The specific steps of performing the dust turbidity reduction analysis include: The optical interference characteristics are fitted and compared with the preset Mie scattering theory attenuation model to calculate the equivalent particle size distribution range of suspended dust in the current tunnel environment and the corresponding optical extinction cross-sectional area parameters, and generate multi-band attenuation compensation coefficients. The multi-band attenuation compensation coefficient is used to perform corresponding amplitude amplification and baseline stretching on different frequency components in the waveform leading edge time series to reconstruct the original infrared absorption spectrum profile features and generate a restored spectrum profile model. Extract the ratio of actual transmitted light intensity to incident light intensity at the target characteristic absorption peak of the reduced spectral profile model, perform a negative logarithmic transformation on the ratio of actual transmitted light intensity to incident light intensity to obtain a pure absorbance parameter, and calculate the actual optical path attenuation parameter.

7. The method for dynamic analysis of mine gas concentration according to claim 1, characterized in that, The specific process of performing the differential calculation and comparison includes: Logarithmic scaling transformation is performed on the actual optical path attenuation parameter and the difference in the reference dust-free absorbance, and the logarithmic amplitude envelope of the two in the frequency domain is extracted to generate a smooth absorption envelope. Perform first-order derivative operations with respect to wavelength on the smooth absorption envelope to convert the slowly varying spectral baseline background into tiny fluctuations near zero, highlighting the peak structure of narrowband gas absorption characteristics and generating first-order derivative spectral curves. The amplitude difference values ​​of the two sets of first derivative spectral curves at corresponding wavelength positions are calculated. The broadband absorption interference caused by ambient temperature and light source fluctuations is filtered out by the extreme value distribution of the difference values ​​to generate the dynamic differential absorption spectrum.

8. A dynamic analysis system for mine gas concentration, characterized in that, The system is used to implement the dynamic analysis method for mine gas concentration according to any one of claims 1-7, and the system comprises: The disturbance topology construction module is used to collect airflow fluctuation parameters of the roadway cross section through a mine microwave anemometer, extract the wind speed change vector from the airflow fluctuation parameters, and perform nonlinear spatiotemporal mapping operation on the wind speed change vector and the original concentration frequency electrical signal collected by the gas detection probe to construct an environmental flow field disturbance topology map. The time constant calculation module is used to extract the absolute pressure distribution gradient of discrete nodes in the environmental flow field disturbance topology map, calculate the air flow resistance attenuation coefficient based on the absolute pressure distribution gradient, substitute the air flow resistance attenuation coefficient into the transport mass conservation partial differential equation for parameter iteration calculation, and calculate the transient lag time constant. The optical path attenuation analysis module is used to perform phase reverse compensation shift on the original concentration frequency electrical signal based on the transient lag time constant, generate a waveform leading edge time series, and use the light spot geometric distortion parameter collected by the scattering dust meter to perform dust turbidity restoration analysis on the waveform leading edge time series to calculate the actual optical path attenuation parameter. The concentration value calculation module is used to obtain the pre-configured baseline dust-free absorbance difference value, perform differential calculation and comparison between the actual optical path attenuation parameter and the baseline dust-free absorbance difference value, and output the dynamic gas concentration value based on the differential calculation and comparison result.

Images