Cosmic ray muon imaging based coal mine underground reservoir water storage performance dynamic monitoring method and device

Abstract

This invention discloses a method and device for dynamic monitoring of the water storage performance of underground coal mine reservoirs based on cosmic ray muon imaging, belonging to the technical field of monitoring the water storage performance of underground coal mine reservoirs. The method includes: acquiring first muon flux data of a target area; determining the time-varying density distribution of each voxel in the target area at each time series based on the first muon flux data; determining the time-varying porosity of each voxel at each time series based on the time-varying density distribution and a preset object density set; calculating the porosity loss rate based on the time-varying porosity of each voxel at each time series; constructing a multidimensional feature vector for each voxel based on the porosity loss rate; performing cluster analysis on the multidimensional feature vector to obtain clustering results, and monitoring the target area based on the clustering results. This method can reflect the spatial distribution and time-varying patterns of internal water storage capacity, enabling reasonable reservoir scheduling and accurate assessment of the reservoir's remaining lifespan.

Description

Technical Field

[0001] This invention belongs to the technical field of monitoring the water storage performance of underground reservoirs in coal mines, and specifically relates to a method and device for dynamic monitoring of the water storage performance of underground reservoirs in coal mines based on cosmic ray muon imaging. Background Technology

[0002] Coal mine underground reservoirs are massive underground water conservancy projects that utilize the voids in the rock mass of the goaf formed during coal mining as water storage space, providing a vital guarantee for water supply for mining operations and residential use. The storage coefficient is a core indicator for measuring the water storage capacity of an underground reservoir, directly determining its usable capacity and regulation ability. However, underground reservoirs are buried hundreds of meters underground, and the internal structure of the goaf is complex and inaccessible, making them typical "black box" systems whose water storage status cannot be directly observed through conventional methods.

[0003] During long-term operation, the water storage capacity of underground reservoirs will continuously decrease. The main reasons include: under the stress of the overlying strata, the collapsed rock mass in the goaf undergoes creep compaction, leading to a gradual reduction in effective porosity; fine particulate matter such as coal slime carried in the mine water gradually accumulates in the porosity, blocking the water storage channels. The above processes are spatially heterogeneous, and the decay rate may vary considerably in different areas.

[0004] Currently, the assessment of the water storage capacity of underground reservoirs mainly relies on pumping tests. Pumping tests are time-consuming and costly, and conducting a single overall scheduling operation severely impacts the reservoir's operation. Furthermore, they can only obtain average parameters of the entire reservoir, failing to reflect the spatial distribution and temporal variation of its internal water storage capacity. This leads to unreasonable scientific scheduling of the reservoir and inaccurate assessment of its remaining lifespan. Summary of the Invention

[0005] This application provides a method and apparatus for dynamic monitoring of the water storage performance of underground coal mine reservoirs based on cosmic ray muon imaging, to solve at least one of the aforementioned problems in the prior art. This application can reflect the spatial distribution and temporal variation of internal water storage capacity, enabling reasonable reservoir scheduling and accurate assessment of the reservoir's remaining lifespan.

[0006] On one hand, embodiments of this application provide a method for dynamic monitoring of the water storage performance of underground reservoirs in coal mines based on cosmic ray muon imaging, the method comprising: Obtain the first muon flux data for the target region; Based on the first muon flux data, the time-varying density distribution of each voxel in the target region at each time series is determined; Based on the time-varying density distribution and the preset object density set, the time-varying porosity of each voxel in each time series is determined. The porosity loss rate is calculated based on the time-varying porosity of each voxel at each time series. Based on the porosity loss rate, a multidimensional feature vector is constructed for each voxel. The multidimensional feature vector includes the average loss rate, the slope of the loss rate change trend, the standard deviation of the loss rate fluctuation, the trend fitting determination coefficient, and the cumulative immersion time. Cluster analysis is performed on the multidimensional feature vectors to obtain clustering results, which are then used to monitor the target region.

[0007] Optionally, determining the time-varying density distribution of each voxel in the target region at each time series based on the first muon flux data includes: After filtering and denoising the first muon flux data, the incident zenith angle and azimuth angle of each muon are reconstructed according to the hit position of each muon event in the multi-layer structure of the detector. Muons with the same zenith angle and azimuth angle are divided into the same path group. Based on the principle of muon transmission imaging, the muon survival rate of each path group is converted into equivalent mass length. The target region is discretized into a voxel grid. A first geometric length matrix is ​​established based on the geometric length of each muon path in each voxel. The first matrix equation is constructed with the equivalent mass length as the observed value and the voxel density as the variable to be solved. The time-varying density distribution of each voxel at each time series is obtained by solving the first matrix equation using an inversion algorithm.

[0008] Optionally, before acquiring the first muon flux data of the target region, the method further includes: Obtain the spatial distribution of the target area, including the flux distribution of the underground reservoir space and the feasible spatial locations of tunnels and chambers; With detector coverage and reception efficiency as optimization objectives, a multi-objective optimization algorithm is used to solve for the optimal number of muon detector arrays, the placement position of each muon detector, the zenith angle, and the azimuth angle. The acquisition of the first muon flux data of the target region includes: Acquire the first muon flux data of the muon detector array in the target region.

[0009] Optionally, before determining the time-varying porosity of each voxel at each time series based on the time-varying density distribution and the preset object density set, the method further includes: Acquire the flux data of the second muon in the drained state and the flux data of the third muon in the saturated state of the target area, respectively; Based on the second muon flux data and the third muon flux data, the initial density distribution and saturated density distribution of each voxel in each time series under different states of the target region are determined. The initial density distribution is used to characterize the density distribution of the target region under the empty state, and the saturated density distribution is used to characterize the density distribution of the target region under the saturated state. Calculate the initial porosity based on the initial density distribution and the saturated density distribution; The density of the rock solid skeleton is determined based on the initial porosity and the initial density distribution. Obtain geological materials and mining engineering plans of the target area; Based on the rock solid skeleton density, initial porosity, geological materials, and mining engineering plan, a baseline water storage capacity distribution map of the target area is constructed. In order to determine the time-varying porosity of each voxel in each time series according to the baseline water storage capacity distribution map, the time-varying density distribution, and the preset object density set.

[0010] Optionally, determining the time-varying porosity of each voxel at each time series based on the baseline water storage capacity distribution map, the time-varying density distribution, and the preset object density set includes: Obtain the water level data of the current time series, and divide the voxel grid into submerged voxels and non-submerged voxels according to the water level data; For voxels in non-submerged areas, the rock solid skeleton density of the voxel is read from the baseline water storage capacity distribution map, and the time-varying porosity is calculated based on the relationship between the time-varying density of the voxel and the rock solid skeleton density. For a voxel in the submerged area, the rock solid skeleton density of the voxel is read from the baseline water storage capacity distribution map, and the time-varying porosity is calculated based on the relationship between the time-varying density of the voxel, the rock solid skeleton density, and the density of water.

[0011] Optionally, constructing a multidimensional feature vector for each voxel based on the porosity loss rate includes: Based on the time-series data of the porosity loss rate of each voxel in each monitoring period, the least squares method is used to perform linear fitting to obtain the slope of the loss rate change trend of each voxel. The slope of the loss rate change trend is used to characterize the acceleration or deceleration trend of the porosity loss of the voxel over time. The arithmetic mean of the time-series data of the loss rate of each voxel is calculated to obtain the average loss rate of each voxel. The standard deviation of the loss rate time series data for each voxel is calculated to obtain the standard deviation of the loss rate fluctuation for each voxel. Calculate the coefficient of determination of the linear fit, which is used to evaluate the reliability of the linear trend.

[0012] Optionally, the step of performing cluster analysis on the multidimensional feature vectors to obtain clustering results includes: The multidimensional feature vector of each voxel is standardized to eliminate the influence of differences in size and order of magnitude between different features on cluster analysis. The standardized multidimensional feature vectors are input into the clustering algorithm to divide all voxels into two clusters; Calculate the average loss rate and mean trend slope for each cluster, and determine the region category corresponding to each cluster based on the average loss rate and mean trend slope. The region category includes normal creep region and abnormal decay region.

[0013] Optionally, after performing cluster analysis on the multidimensional feature vectors to obtain the clustering results, the method further includes: The time-varying void ratio distributions of each time series are stored in time series to construct a void ratio time series database; A training sample set is constructed using the porosity time series data of each voxel in the porosity time series database for multiple consecutive historical monitoring periods, the water level data for the corresponding time period, the water injection method code, and the stress data of the overlying strata. Using the input features of historical periods in the training sample set as the model input, and using the porosity of the corresponding voxel in a future preset time window as the prediction label, the time series prediction model is trained so that the model learns the time series mapping relationship between porosity evolution and water level, water injection method and overlying strata stress. The water injection mode codes corresponding to the candidate scheduling schemes and the stress sequence of the overlying strata are input into the trained time-series prediction model to predict the evolution trajectory of the porosity of each voxel under each candidate scheduling scheme, calculate the change curve of the overall water storage coefficient under each candidate scheme, and output the scheme comparison results.

[0014] Optionally, the accumulation yields the overall water storage coefficient variation curves for each candidate scheme, and the scheme comparison results are output, including: For each candidate scheduling scheme, the overall water storage coefficient is calculated step by step by the predicted porosity of each voxel at each time step, and the curve of the overall water storage coefficient changing with time under each candidate scheme is generated. The overall water storage coefficient change curves of different candidate schemes are superimposed on the same coordinate graph, and the predicted water storage coefficient values ​​at preset time nodes and the percentage decrease relative to the initial water storage coefficient are marked for engineers to compare and select the best water replenishment scheduling strategy.

[0015] On the other hand, this application provides a dynamic monitoring device for the water storage performance of underground coal mine reservoirs based on cosmic ray muon imaging. The device includes: The acquisition module is used to acquire the first muon flux data of the target region; The first determining module is used to determine the time-varying density distribution of each voxel in the target region at each time sequence based on the first muon flux data. The second determining module is used to determine the time-varying porosity of each voxel in each time series based on the time-varying density distribution and the preset object density set. The calculation module is used to calculate the porosity loss rate based on the time-varying porosity of each voxel at each time series. The construction module is used to construct a multidimensional feature vector for each voxel based on the porosity loss rate. The multidimensional feature vector includes the average loss rate, the slope of the loss rate change trend, the standard deviation of the loss rate fluctuation, the trend fitting determination coefficient, and the cumulative immersion time. The analysis module is used to perform cluster analysis on the multidimensional feature vectors to obtain clustering results, so as to monitor the target region based on the clustering results.

[0016] In another aspect, embodiments of this application provide an electronic device, the device comprising: a processor and a memory storing computer program instructions; When the processor executes the computer program instructions, it implements the dynamic monitoring method for water storage performance of underground coal mine reservoirs based on cosmic ray muon imaging as described in the first aspect.

[0017] In another aspect, embodiments of this application provide a computer storage medium storing computer program instructions, which, when executed by a processor, implement the dynamic monitoring method for water storage performance of underground coal mine reservoirs based on cosmic ray muon imaging as described in the first aspect.

[0018] In another aspect, embodiments of this application provide a computer program product in which the instructions are executed by the processor of an electronic device, causing the electronic device to perform the dynamic monitoring method for water storage performance of underground coal mine reservoirs based on cosmic ray muon imaging as described in the first aspect.

[0019] Compared with existing technologies, the dynamic monitoring method and apparatus for water storage performance of underground coal mine reservoirs based on cosmic ray muon imaging in this application can determine the time-varying density distribution and time-varying porosity of each voxel in the target area at each time series by acquiring the first muon flux data of the target area, thereby quantifying the porosity loss rate and constructing a multidimensional feature vector for each voxel. Through cluster analysis of the multidimensional feature vector, the water storage performance decay characteristics of different areas inside the underground reservoir can be identified. This overcomes the limitations of traditional pumping tests that cannot reflect spatial distribution and time-varying patterns, and provides a refined and real-time basis for the scientific scheduling and remaining life assessment of the reservoir, so as to reasonably schedule the reservoir and accurately assess its remaining life. Attached Figure Description

[0020] To more clearly illustrate the technical solutions in the embodiments of this specification or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in the embodiments of this specification. For those skilled in the art, other drawings can be obtained based on these drawings. Figure 1 This is a flowchart illustrating a method for dynamic monitoring of water storage performance of underground coal mine reservoirs based on cosmic ray muon imaging, provided in one embodiment of this application. Figure 2 This is a field schematic diagram of a method for dynamic monitoring of the water storage performance of underground reservoirs in goaf areas, provided in one embodiment of this application. Figure 3 This is a field schematic diagram of a method for dynamic monitoring of the water storage performance of underground reservoirs in goaf areas, provided in another embodiment of this application. Figure 4 This is a schematic diagram of the structure of a dynamic monitoring device for the water storage performance of underground coal mine reservoirs based on cosmic ray muon imaging, provided in another embodiment of this application. Figure 5 This is a schematic diagram of the structure of an electronic device provided in another embodiment of this application.

[0021] Figure label: 1-Müller detector; 2-Non-submerged area of ​​reservoir; 3-Submerged area of ​​reservoir; 4-Rock strata outside reservoir; 5-Coal pillar dam; 6-Artificial dam; 7-Adjacent roadway; 8-Breakhole layout point; 9-Müller transmission track. Detailed Implementation

[0022] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0023] To facilitate understanding of the embodiments of this application, further explanation and description will be provided below with reference to the accompanying drawings and specific embodiments. These embodiments do not constitute a limitation on the embodiments of this application. In the drawings, the dimensions and relative dimensions of components may be exaggerated for clarity and / or descriptive purposes. When exemplary embodiments can be implemented differently, a specific process sequence may be performed in a different order than that described. For example, two consecutively described processes may be performed substantially simultaneously or in the reverse order of their description. Furthermore, the same reference numerals denote the same components.

[0024] The terminology used herein is for the purpose of describing particular embodiments and is not intended to be limiting. As used herein, unless the context clearly indicates otherwise, the singular forms “a” and “the” are intended to include the plural forms as well. Furthermore, when the terms “comprising” and / or “including” and variations thereof are used in this specification, it indicates the presence of the stated features, integrals, steps, operations, parts, components, and / or groups thereof, but does not exclude the presence or addition of one or more other features, integrals, steps, operations, parts, components, and / or groups thereof. It should also be noted that, as used herein, the terms “substantially,” “about,” and other similar terms are used as approximate terms rather than as terms of degree, thus explaining the inherent biases in measurements, calculated values, and / or provided values ​​that would be recognized by one of ordinary skill in the art.

[0025] To address the problems in existing technologies, this application provides a method and apparatus for dynamic monitoring of the water storage performance of underground coal mine reservoirs based on cosmic ray muon imaging. In this application embodiment, by acquiring the first muon flux data of the target area, the time-varying density distribution and time-varying porosity of each voxel at different time series in the target area are determined, thereby quantifying the porosity loss rate and constructing a multidimensional feature vector for each voxel. Through cluster analysis of this multidimensional feature vector, the water storage performance degradation characteristics of different areas within the underground reservoir can be identified. This overcomes the limitations of traditional pumping tests in reflecting spatial distribution and time-varying patterns, providing a refined and real-time basis for the scientific scheduling and remaining life assessment of the reservoir, enabling reasonable reservoir scheduling and accurate assessment of the reservoir's remaining life.

[0026] The following section first introduces the method for dynamic monitoring of water storage performance of underground coal mine reservoirs based on cosmic ray muon imaging, as provided in the embodiments of this application.

[0027] Figure 1 This illustration shows a flowchart of a method for dynamic monitoring of water storage performance in underground coal mine reservoirs based on cosmic ray muon imaging, according to an embodiment of this application. Figure 1 As shown, the dynamic monitoring method for the water storage performance of underground coal mine reservoirs based on cosmic ray muon imaging can include S101-S106: S101, Obtain the first muon flux data of the target region.

[0028] In this embodiment, the first muon flux data of the target region can be obtained by deploying a muon detector array at a specific location in the target region. After the deployment is completed, the muon detector array needs to be run continuously for a period of time to record cosmic ray muon events passing through the detector. These muon events can contain information such as the incident direction and energy of the muons. By collecting and preliminarily processing these raw data, the first muon flux data of the target region can be obtained.

[0029] S102, based on the first muon flux data, determine the time-varying density distribution of each voxel in the target region at each time series.

[0030] In this embodiment, after acquiring the first muon flux data, density inversion can be performed using the basic principle of muon transmission imaging. Specifically, the target region can be divided into a voxel grid of a preset size. By statistically analyzing the changes in muon flux passing through the measurement region and combining this with the attenuation patterns of muons in different density media, the average density of the voxel can be preliminarily estimated. By repeating this process at different time points, the density distribution of each voxel at each time series can be obtained.

[0031] S103, based on the time-varying density distribution and the preset object density set, determines the time-varying porosity of each voxel at each time series.

[0032] In some embodiments, after obtaining the time-varying density distribution of each voxel, the porosity can be calculated by combining it with a preset object density set, which typically includes the main materials constituting the underground reservoir, such as the density of the rock skeleton and the density of water. Assuming that the voxel contains only the rock skeleton and water, the voxel density can be expressed as the product of the rock skeleton density and its volume fraction, plus the product of the water density and its volume fraction. Using simple principles of mass or volume conservation and the known densities of the rock skeleton and water, the time-varying porosity can be deduced from the voxel's time-varying density. For example, a fixed rock skeleton density and water density can be set, assuming that the voxel contains only these two substances.

[0033] S104, calculate the porosity loss rate based on the time-varying porosity of each voxel at each time series.

[0034] In this embodiment of the application, the porosity loss rate can be calculated by comparing the porosity values ​​of two adjacent time series.

[0035] Specifically, for each voxel, the average rate of change of porosity over that time period can be obtained by subtracting the porosity of the current time period from the porosity of the previous time period and then dividing by the time interval. If this rate is positive, it indicates a loss of porosity. For example, the difference in porosity between two consecutive monitoring periods can be simply calculated.

[0036] S105, based on the porosity loss rate, constructs a multidimensional feature vector for each voxel.

[0037] In this embodiment of the application, after calculating the porosity loss rate of each voxel, multiple features can be further extracted to construct a multidimensional feature vector, wherein the multidimensional feature vector includes the average loss rate, the slope of the loss rate change trend, the standard deviation of the loss rate fluctuation, the trend fitting determination coefficient, and the cumulative immersion time.

[0038] It is worth noting that the average loss rate is the arithmetic mean of the porosity loss rate of the voxel over all monitoring periods, which is used to characterize the overall rate of porosity decay of the voxel. The slope of the loss rate trend was obtained by linearly fitting the loss rate time series using the least squares method, and was used to characterize the direction of change of the porosity loss of the voxel over time. The standard deviation of the loss rate fluctuation is the standard deviation of the time series data of the loss rate, which is used to characterize the stability of the voxel porosity decay process; The trend fitting determination coefficient is the same as the linear fitting determination coefficient, used to characterize the reliability of the monotonic trend of the voxel porosity loss. The cumulative immersion time is calculated by adding the duration of water submersion of each voxel based on real-time water level monitoring data. This is used to characterize the magnitude of the water-rock interaction and sedimentation potential of that voxel.

[0039] S106, perform cluster analysis on the multidimensional feature vectors to obtain clustering results, and monitor the target area based on the clustering results.

[0040] In this embodiment, the K-means algorithm can be used. The number of clusters is preset, and voxels are iteratively assigned to the nearest cluster centers until the cluster centers no longer change significantly. After clustering, each voxel is assigned to a specific cluster. Based on these clustering results, the target area can be divided into different regions. For example, regions with rapidly declining water storage performance or relatively stable regions can be identified based on the characteristics of the clusters, thereby enabling dynamic monitoring of the target area.

[0041] In this embodiment, by acquiring the first muon flux data of the target area, the time-varying density distribution and time-varying porosity of each voxel in each time series of the target area are determined, and then the porosity loss rate is quantified. A multidimensional feature vector of each voxel is constructed. Through cluster analysis of the multidimensional feature vector, the water storage performance decay characteristics of different areas inside the underground reservoir can be identified. This overcomes the limitation that traditional pumping tests cannot reflect spatial distribution and time-varying patterns, and provides a refined and real-time basis for the scientific scheduling and remaining life assessment of the reservoir, so as to reasonably schedule the reservoir and accurately assess the remaining life of the reservoir.

[0042] In some embodiments, in order to enable accurate monitoring of the underground reservoir, the method may further include, prior to S101: Obtain the spatial distribution of the target area, including the flux distribution of underground water reservoir space and the feasible spatial locations of tunnels and chambers; With detector coverage and reception efficiency as optimization objectives, a multi-objective optimization algorithm is used to solve for the optimal number of muon detector arrays, the placement position of each muon detector, the zenith angle, and the azimuth angle. S101 may include: Acquire the first muon flux data of the muon detector array in the target region.

[0043] In this embodiment, in order to make the obtained first muon flux data more accurately reflect the water storage capacity of the underground reservoir, it is necessary to determine the optimal layout scheme.

[0044] For details, see Figure 2 and Figure 3 Based on mine geological data and mining engineering plans, a three-dimensional geological model can be constructed, including the underground reservoir goaf, non-reservoir strata 4, coal pillar dam 5, artificial dam 6, and adjacent roadways 7. The underground reservoir goaf can include the non-submerged area 2 and the submerged area 3. Corresponding prior density reference values ​​are assigned to each stratum according to the mine geological report. Subsequently, using the Monte Carlo simulation method based on the Geant4 toolkit, a large number of incident muon events conforming to the actual energy spectrum and angular distribution characteristics of cosmic ray muons are generated at the surface. The underground reservoir storage is varied, with five cases: 0%, 25%, 50%, 75%, and 100% storage capacity, and four variations are performed. The transport process of each event is tracked in the constructed three-dimensional geological model and the changes in underground reservoir density, and its flux is recorded.

[0045] To achieve optimal detection results, the number, attitude, and spatial location of the detectors need to be optimized. In this embodiment, the muon detector array is divided into a detection group and a reference group. The detection group is planned to be deployed in tunnels, chambers, or specially constructed boreholes near the underground reservoir to receive attenuated muon signals that penetrate the reservoir area; the reference group is deployed on the surface to synchronously record the same-direction incident muon flux that does not penetrate the reservoir area, thereby eliminating errors in muon flux fluctuations caused by environmental factors such as atmospheric pressure and temperature changes.

[0046] When determining the optimal deployment scheme, multiple objectives must be comprehensively considered, including detection sensitivity, effective track rate, and deployment cost. Let the set of candidate deployment points be... The decision variables to be solved are the number of detectors n and the deployment location of each detector. Zenith angle and azimuth A complete deployment plan is denoted as The objective function F is used to evaluate the merits of the entire deployment scheme. Optionally, F is defined as:

[0047] in: The overall detection sensitivity of the entire detection array to density changes in the reservoir area is determined by statistical results from Monte Carlo simulations, and can be selected as the average of the relative change rates of muon flux during four water level changes in the reservoir area. , The effective muon flux received; The effective track rate during emptying is the proportion of track flux penetrated by muons within the reservoir area to the total received muon flux. The normalized total cost of equipment and construction; The construction difficulty level of the deployment scheme is scored by engineering experts based on factors such as borehole inclination angle, rock strata stability, and chamber size; α, β, γ, and δ are the weight coefficients of each sub-objective, satisfying α+β+γ+δ=1 (e.g., α=0.4, β=0.3, γ=0.2, δ=0.1). The candidate space is searched using algorithms such as genetic algorithms or particle swarm optimization to... Largest deployment plan This allows us to obtain the optimal number, attitude, and spatial position of the detectors.

[0048] As an example of a deployment result, such as Figure 2 and Figure 3 As shown, after optimization calculations, it was determined that four drilling points 8 should be constructed downwards in the tunnel 7 near the reservoir, and a total of four muon detectors 1 should be installed.

[0049] After determining the deployment scheme of the muon detector array, the water storage performance of the underground reservoir can be monitored by acquiring data from the muon detector array 1.

[0050] In some embodiments, S102 may include: After filtering and denoising the first muon flux data, the incident zenith angle and azimuth angle of each muon are reconstructed based on the hit position of each muon event in the detector's multi-layer structure. Muons with the same zenith angle and azimuth angle are divided into the same path group. Based on the principle of muon transmission imaging, the muon survival rate of each path group is converted into equivalent mass length. The target region is discretized into a voxel grid. The first geometric length matrix is ​​established based on the geometric length of each muon path in each voxel. The first matrix equation is constructed with the equivalent mass length as the observed value and the voxel density as the variable to be solved. The time-varying density distribution of each voxel at each time series is obtained by solving the first matrix equation using an inversion algorithm.

[0051] In this embodiment, the first muon flux data is first filtered and denoised. The raw data received by the muon detector may contain background noise, scattered muons, and other ineffective signals, so filtering and denoising are necessary to improve data quality. The filtering and denoising methods can include various techniques based on energy thresholds, time coincidence, and trajectory fitting goodness of fit. For example, the response signals of each detection unit in the multi-layer structure of the detector can be used to identify and remove abnormal trajectories or noise points through a trajectory reconstruction algorithm.

[0052] Subsequently, by analyzing the coordinates of the Mun's impact point on different detection layers, the zenith angle (the angle between the Mun and the vertical direction) and azimuth angle (the projection angle in the horizontal plane) of the Mun entering the detector can be accurately calculated using geometric principles or trajectory fitting algorithms.

[0053] Next, since the attenuation of muons when passing through the target area is related to the density of the matter they pass through and the path length, and muons with the same incident direction can be considered to have passed through similar paths, muon events with the same or similar incident zenith angle and azimuth angle can be classified into the same path group.

[0054] Based on this, the muon survival rate of each path group is converted into equivalent mass length according to the principle of muon transmission imaging. When muons pass through matter, they lose energy due to mechanisms such as ionization loss and radiation loss, resulting in flux attenuation. The degree of attenuation is positively correlated with the density of the matter and the length of the path traversed. The muon survival rate refers to the ratio of the number of muons received by the detector after passing through the target region in a specific path group to the number of muons incident on the target region. This ratio directly reflects the degree of attenuation experienced by muons during penetration. The equivalent mass length is obtained by inverting the muon survival rate through the muon attenuation model, and it represents the total mass thickness of the matter traversed by muons along a specific path.

[0055] Furthermore, the target region is discretized into a voxel mesh. To achieve a fine description of the density distribution within the target region, the continuous target region needs to be divided into a series of discrete voxels with fixed geometry and size.

[0056] Then, a first geometric length matrix is ​​constructed based on the geometric length of each muon path within each voxel. For each muon path group, its muons will traverse multiple voxels in the target region. The first geometric length matrix records the geometric length traversed by each muon path group within each voxel.

[0057] The first matrix equation is constructed using the equivalent mass length as the observed value and voxel density as the variable to be solved. Based on the principle of muon transmission imaging, the equivalent mass length (observed value) of each path group can be expressed as the sum of the products of the densities of all voxels traversed by the path group and their geometric lengths within the corresponding voxels. Combining the equations of all path groups forms the first matrix equation, where the observed value is the equivalent mass length and the variable to be solved is the voxel density.

[0058] Finally, the first matrix equation is solved using an inversion algorithm to obtain the time-varying density distribution of each voxel at each time series. The first matrix equation is typically a large, overdetermined or underdetermined system of linear equations that requires an inversion algorithm for solution. Commonly used inversion algorithms include least squares methods (such as Tikhonov regularization), iterative algorithms (such as ART, SIRT, CGLS, etc.), or statistical methods. Since muon flux data is acquired continuously over time, the above inversion calculations can be performed periodically. Each calculation is based on muon flux data within a specific time period (one time series), thus obtaining the density distribution of each voxel within the target region at that time series. By repeating this process at different time series, the time-varying density distribution of each voxel can be obtained, reflecting the change in the water storage capacity of the underground reservoir over time.

[0059] By filtering and denoising the first muon flux data and accurately reconstructing the muon incident direction, muons with similar paths are grouped into path groups, laying the data foundation for subsequent density inversion. Based on the principle of muon transmission imaging, muon survival rate is converted into equivalent mass length, providing key observational values ​​reflecting the mass thickness of the material. Furthermore, by discretizing the target region into a voxel grid and establishing a first geometric length matrix and a first matrix equation, the complex physical problem is transformed into a computable mathematical model. Finally, by solving this matrix equation through an inversion algorithm, the time-varying density distribution of each voxel at each time series can be obtained stably and accurately, thus providing reliable and refined basic data for subsequent porosity calculation and water storage performance assessment, significantly improving the accuracy and practicality of the monitoring method.

[0060] As an example, during the monitoring process, the initial muon flux measurement is initiated by the detection and reference groups during the initial construction phase of the underground reservoir or when it is in a state of being emptied due to artificial scheduling. The data acquisition cycle is determined based on the detector's receiving efficiency and the required statistical accuracy, typically ranging from several days to several weeks, with a possible cycle of 15 days. After data acquisition, the received muon events are filtered and denoised, retaining only the valid muon events.

[0061] Muon detectors record the impact points of muons through a multi-layered structure, and then analyze the impact points from two layers. Given the interlayer spacing H, the direction vector based on the local coordinate system of the muon detector is calculated. The transformation between the detector coordinate system (X′, Y′, Z′) and the geographic coordinate system (X, Y, Z) is achieved by a rotation matrix R(α, β, γ) consisting of the angles between the two coordinate systems. The direction vector in the geographic coordinate system is obtained through this coordinate transformation. Then calculate the true zenith angle of the muon. and azimuth The muon transmission track 9 was obtained, and muons with the same incident direction were grouped together.

[0062] Based on the mature muon transmission imaging principle, the underground modeling space is discretized into A voxel units, where the effective water storage space of the underground reservoir contains N voxels. Local optimization algorithms, global optimization algorithms, or machine learning algorithms are used to solve the matrix equations, and the time-varying density distribution of each voxel under the condition of the underground reservoir being emptied is obtained through inversion.

[0063] In other embodiments, after reconstructing the incident zenith angle and azimuth angle of each muon based on the impact location of each muon event in the detector's multilayer structure, the method further includes: Acquire multiple muon events recorded by the muon detector, each muon event containing incident direction information and energy information; Based on energy information, multiple muon events are divided into K energy spectrum channels, where K ≥ 3; Transmission imaging inversion was performed on muon events in each energy spectral channel to obtain the equivalent density distribution of each voxel in the target region under each energy spectral channel. For each voxel, its equivalent density in K energy spectral channels is combined to form the energy spectral fingerprint vector of that voxel. The energy spectrum fingerprint vector of each voxel is matched with the pre-generated standard energy spectrum fingerprint library to determine the water storage medium change type corresponding to each voxel based on the matching results. Output the identification results of the water storage medium change type for each voxel.

[0064] In this embodiment, when recording muons passing through its multi-layered structure, the muon detector can not only reconstruct the incident direction (zenith angle and azimuth angle) of the muons, but also measure the energy lost by the muons in the detector or their total energy.

[0065] Specifically, since the interaction between muons and matter (such as ionization loss and radiation loss) is energy-dependent, muons of different energies exhibit significantly different attenuation characteristics when penetrating matter of different densities and atomic numbers. Therefore, it is possible to monitor the water storage performance of underground reservoirs while also distinguishing matter of different energies.

[0066] To analyze the interaction between muons and matter more precisely, this embodiment divides the detected muon events into K different energy spectrum channels based on their energy information. These K different energy spectrum channels provide sufficient resolution to distinguish the differences in the responses of muons of different energies to matter. By setting a series of energy thresholds or energy ranges, muon events can be categorized into corresponding energy spectrum channels; for example, they can be divided into low-energy, medium-energy, and high-energy channels, or even more detailed energy ranges.

[0067] After dividing muon events into different energy spectrum channels, inversion is performed independently for muon events within each channel using the principle of muon transmission imaging. Since muons of different energies have different attenuation patterns when penetrating matter, the voxel equivalent density distribution obtained from the inversion in each energy spectrum channel will reflect the response characteristics of the target region's matter to muons of that specific energy.

[0068] Subsequently, the equivalent density values ​​obtained by inversion of the same voxel under K different energy spectral channels are combined to form a K-dimensional energy spectral fingerprint vector. This energy spectral fingerprint vector uniquely characterizes the material composition and state of the voxel because it integrates the response information of matter to muons of different energies. For example, if there are K energy spectral channels, the energy spectral fingerprint vector of each voxel can be represented as [ρ eq1 ,ρ eq2 ,...,ρ eqK ], where ρ eqK This represents the equivalent density of the voxel in the Kth energy channel.

[0069] To identify the material type or variation type of a voxel, this embodiment pre-establishes a standard energy spectrum fingerprint database. This database contains theoretical or experimental energy spectrum fingerprints of known substances (such as different rocks, coal, water, air, etc.) under the influence of muons at different energies. By performing pattern matching between the measured energy spectrum fingerprint vectors of each voxel and the standard fingerprints in the database (e.g., by calculating Euclidean distance, cosine similarity, or using machine learning algorithms), the material type or variation type of the voxel can be identified, for example, determining whether it is water intrusion into rock pores and fissures, coal seam deterioration, or cavity formation.

[0070] Finally, the type of water storage medium change identified for each voxel (e.g., "water-saturated rock," "dry coal seam," "water-coal mixture," "cavity," etc.) will be output. These results can be presented in the form of visualization maps, data reports, or 3D model annotations, providing engineers with intuitive and detailed information on underground reservoir medium changes.

[0071] Through the above technical solution, this application overcomes the limitation of accurately identifying the type and changing properties of water storage media by relying solely on a single density value. By performing energy resolution on muon events and conducting independent transmission imaging inversion for different energy spectral channels, a unique energy spectral fingerprint vector can be constructed for each voxel. This energy spectral fingerprint vector contains comprehensive response information of matter to muons of different energies, providing richer and more refined material composition and state characteristics than a single density value. By performing pattern matching of this energy spectral fingerprint vector with a pre-established standard energy spectral fingerprint library, the changing type of water storage media corresponding to each voxel in the target area can be accurately identified. This enables the monitoring system not only to sense changes in density but also to diagnose the specific physical processes causing these changes, such as distinguishing between increased water saturation, rock skeleton degradation, or intrusion of other substances. Therefore, this application significantly improves the accuracy and specificity of dynamic monitoring of the water storage performance of underground reservoirs in coal mines, providing more reliable and in-depth technical support for the refined management, risk assessment, and formulation of water replenishment scheduling strategies for underground reservoirs.

[0072] In some other embodiments, prior to S103, the method may further include: Acquire the flux data of the second muon in the drained state and the flux data of the third muon in the saturated state of the target area, respectively; Based on the second muon flux data and the third muon flux data, the initial density distribution and saturated density distribution of each voxel in the target region at different times under different states are determined. The initial density distribution is used to characterize the density distribution of the target region under the empty state, and the saturated density distribution is used to characterize the density distribution of the target region under the saturated state. The initial porosity is calculated based on the initial density distribution and the saturated density distribution. The density of the rock solid skeleton is determined based on the initial porosity and initial density distribution; Obtain geological materials and mining engineering plans for the target area; Based on the density of the rock solid skeleton, initial porosity, geological materials, and mining engineering plan, a baseline water storage capacity distribution map of the target area is constructed. Based on the baseline water storage capacity distribution map, time-varying density distribution, and preset object density set, the time-varying porosity of each voxel at each time series is determined.

[0073] In this embodiment, the water storage capacity of the underground reservoir can be assessed at the initial stage of its construction, serving as the basis for monitoring the initial water storage performance.

[0074] Specifically, the initial density distribution can be determined using the methods described above for determining time-varying density distribution. and saturated density distribution Based on the density superposition principle, the density in the saturated state can be considered as a result of the combined contribution of the solid framework and pore water. Therefore, an equation can be established to calculate the initial porosity on a voxel-by-voxel basis: ,in, This represents the average density of the mine water. From this, the density of the rock solid framework for each voxel can be deduced. .

[0075] Then, by using the acquired geological materials and mining engineering plan, a three-dimensional geological model was obtained. The density of the solid rock skeleton and the initial porosity were then assigned to the three-dimensional geological model to construct a baseline water storage capacity distribution map.

[0076] Furthermore, formulas can be used Calculate the initial storage coefficient of the entire reservoir, where Let be the volume of the i-th voxel.

[0077] After the reservoir enters its long-term operation phase, a time window where the water level is relatively stable is selected. In this embodiment, submersible water level gauges are installed in the reservoir area to monitor water level data in real time. The data collection frequency is once per hour. The system continuously calculates the maximum water level change over the past 24 hours. .when When the water level reaches a certain value (e.g., 0.25 meters), the system determines that the water level has entered a stable state and automatically triggers a 7-day detection cycle. If the water level changes by more than 0.25 meters during the 7-day collection period, the collection is interrupted, and the system waits for the next stable period to restart. After each collection, the time-varying density of each voxel in the reservoir area is inverted using multi-directional muon attenuation data. .

[0078] In other embodiments, the time-varying porosity of each voxel at each time series is determined based on a baseline water storage capacity distribution map, a time-varying density distribution, and a preset object density set, including: Obtain the current time series water level data, and divide the voxel grid into submerged voxels and non-submerged voxels based on the water level data; For voxels in non-submerged areas, the rock solid skeleton density of the voxel is read from the baseline water storage capacity distribution map, and the time-varying porosity is calculated based on the relationship between the time-varying density of the voxel and the rock solid skeleton density. For a voxel in the submerged area, the rock solid skeleton density of the voxel is read from the baseline water storage capacity distribution map, and the time-varying porosity is calculated based on the relationship between the time-varying density of the voxel, the rock solid skeleton density, and the density of water.

[0079] In this embodiment, to accurately determine the changed porosity, the solution domain needs to be divided into a non-submerged region 2 and a submerged region 3 based on the current water level boundary conditions of the reservoir. For the voxels in the non-submerged region located above the current water level, their pores are filled with air, and the air density is negligible, thus satisfying the condition... For voxels in the submerged area located below the current water level, their density is composed of both the framework and water, satisfying the following conditions: Due to the density of the solid skeleton By solving the above equations simultaneously, with the invariants known, we can obtain different time series. Time-varying porosity distribution of voxels in the lower reservoir area By obtaining the results of each calculation By comparing the calculations, the porosity loss rate of each voxel can be calculated simultaneously.

[0080] In some embodiments, in S105, the porosity field of each time voxel at the previous time step can be... With porosity field Perform voxel-by-voxel difference operations to obtain the porosity loss rate of each voxel. All timing sequences , The corresponding time, water injection method, and water level boundary conditions are stored together to form a multi-dimensional holographic time series database containing spatial coordinates, time series, and changes in physical quantities.

[0081] In some other embodiments, S105 may include: Based on the time-series data of the porosity loss rate of each voxel in each monitoring period, the least squares method is used to perform linear fitting to obtain the slope of the loss rate change trend of each voxel. The slope of the loss rate change trend is used to characterize the acceleration or deceleration trend of the porosity loss of the voxel over time. The arithmetic mean of the time-series data of the loss rate of each voxel is calculated to obtain the average loss rate of each voxel. The standard deviation of the loss rate time series data for each voxel is calculated to obtain the standard deviation of the loss rate fluctuation for each voxel. Calculate the coefficient of determination for the linear fit. The coefficient of determination is used to evaluate the reliability of the linear trend.

[0082] In this embodiment, firstly, based on the time-series data of the porosity loss rate of each voxel in each monitoring period, a linear fit is performed using the least squares method to obtain the slope of the loss rate change trend for each voxel. This slope is a key indicator; its sign and magnitude can intuitively characterize the acceleration or deceleration trend of the voxel's porosity loss over time. For example, a positive slope indicates that the loss rate is accelerating, a negative slope indicates that the loss rate is decelerating, and a slope close to zero indicates that the loss rate is relatively stable. The least squares method can effectively extract this linear trend from discrete time-series data.

[0083] Secondly, the arithmetic mean of the time-series loss rate data for each voxel is calculated to obtain the average loss rate for each voxel. The average loss rate provides the overall average level of porosity loss for that voxel throughout the entire monitoring period and is a fundamental indicator for measuring its macroscopic attenuation.

[0084] Furthermore, the standard deviation of the time-series data on the loss rate of each voxel is calculated to obtain the standard deviation of the loss rate fluctuation for each voxel. The standard deviation is a statistic that measures the degree of dispersion of data. The larger the standard deviation of the loss rate fluctuation, the stronger the fluctuation of the porosity loss rate between different monitoring periods, which may indicate more complex physicochemical changes inside the water storage medium or external environmental disturbances.

[0085] Finally, the coefficient of determination (R-squared) for the linear fit is calculated. The R-squared is a value between 0 and 1, used to evaluate the reliability of the linear trend. The closer the R-squared is to 1, the stronger the explanatory power of the linear trend line fitted by the least squares method for the actual data points; that is, the closer the change trend of the porosity loss rate of the voxel is to linear, and the higher the reliability of its trend prediction. Conversely, a low R-squared may mean that the change in the porosity loss rate is not a simple linear relationship and may be influenced by more nonlinear factors.

[0086] Through the aforementioned technical solution, this application no longer relies solely on the single porosity loss rate. Instead, it constructs a multi-dimensional feature vector by introducing the slope of the loss rate change trend, the average loss rate, the standard deviation of the loss rate fluctuation, and the trend fitting determination coefficient. This multi-dimensional analysis method enables a more refined and comprehensive characterization of the degradation process of the water storage medium in coal mine underground reservoirs. The average loss rate provides a quantification of the overall degradation, the slope of the loss rate change trend reveals the dynamic acceleration or deceleration characteristics of the degradation process, the standard deviation of the loss rate fluctuation quantifies the stability or complexity of the degradation process, and the trend fitting determination coefficient assesses the reliability of these trend judgments. The comprehensive application of these features greatly enhances the depth and accuracy of understanding the dynamic changes in the water storage performance of underground reservoirs, enabling more effective identification of degradation patterns in different areas and prediction of their future evolution trends. This provides more scientific and reliable data support for risk assessment, maintenance strategy formulation, and water replenishment scheduling optimization of underground reservoirs.

[0087] In some other embodiments, S106 may include: The multidimensional feature vector of each voxel is standardized to eliminate the influence of differences in size and order of magnitude between different features on cluster analysis. The standardized multidimensional feature vectors are input into the clustering algorithm to divide all voxels into two clusters; Calculate the average loss rate and mean trend slope for each cluster, and determine the region category corresponding to each cluster based on the average loss rate and mean trend slope. The region category includes normal creep region and abnormal decay region.

[0088] In this embodiment, the constructed multidimensional feature vectors characterize the porosity evolution behavior of each voxel from five dimensions: attenuation intensity, trend, stability, reliability, and water-rock interaction time. To eliminate the influence of different feature dimensions and orders of magnitude, Z-score standardization is performed on each feature of all voxels, so that the mean of each feature is 0 and the standard deviation is 1.

[0089] Specifically, Z-score standardization is achieved by subtracting the mean of each feature value from its mean and then dividing by the standard deviation of that feature, resulting in a mean of 0 and a standard deviation of 1. Min-Max standardization, on the other hand, scales the data to a specific range of [0,1] or [-1,1] by subtracting the minimum value of each feature value from its minimum value and then dividing by the difference between the maximum and minimum values. The choice of a suitable standardization method depends on the distribution characteristics of the data and the requirements of the clustering algorithm.

[0090] The K-means clustering algorithm (K=2) is used to divide voxels within a region into two clusters. The specific process of K-means is as follows: two initial cluster centers are randomly selected, and the assignment and update are performed iteratively—each voxel is assigned to the nearest cluster center (the distance can be standardized Euclidean distance), and then the center of each cluster is recalculated until the change in cluster centers is less than a threshold or the maximum number of iterations is reached. Since K-means is sensitive to the initial centers, it can be run 20 times to select the result with the minimum sum of squared errors within the cluster.

[0091] After clustering to obtain two clusters, it is necessary to determine which cluster represents the normal creep region and which represents the abnormal decay region. The determination rule is as follows: calculate the average loss rate and the mean trend slope of each cluster. The cluster with a smaller average loss rate and a trend slope close to zero is determined to be the normal creep region; the cluster with a larger average loss rate and a positive trend slope is determined to be the abnormal decay region. If the average loss rates of the two clusters are similar, their trend slopes are further compared; the one with a significantly positive slope is the abnormal region, and the one with a slope close to zero is the normal region. The physical basis of this determination rule is that abnormal decay is mainly caused by siltation or water-rock weakening, manifested as a continuous acceleration of porosity loss, while the loss rate of normal compaction, whether dry or wet compaction, tends to stabilize or slow down.

[0092] After completing the above clustering and determination, the category label of each voxel and its corresponding immersion area type are mapped back to three-dimensional spatial coordinates.

[0093] It is worth noting that the specific thresholds in the above judgment rules can be dynamically adjusted according to the actual monitoring data, and K-means clustering can be replaced by other data analysis methods.

[0094] By standardizing the multidimensional feature vectors using the above technical solution, the influence of differences in the dimensions and orders of magnitude of different features on clustering analysis is eliminated, ensuring that all features have fair weights during the clustering process, thereby improving the accuracy and reliability of the clustering results. Inputting the standardized feature vectors into the clustering algorithm and dividing them into two clusters effectively and automatically distinguishes voxels in the target area into two main performance states. Furthermore, by calculating the average loss rate and mean trend slope of each cluster, the region is classified as either a normal creep region or an abnormal decay region, making the monitoring results more interpretable and instructive. This classification method can clearly identify areas in underground reservoirs where water storage performance is evolving normally and areas at risk of accelerated decay, providing a scientific basis for the refined management and early warning of coal mine underground reservoirs, and helping to take timely intervention measures to prevent further deterioration of water storage performance.

[0095] In some other embodiments, after S106, the method may further include: The time-varying void ratio distributions of each time series are stored in time series to construct a void ratio time series database; A training sample set was constructed using porosity time series data of each voxel in the porosity time series database for multiple consecutive historical monitoring periods, water level data for the corresponding time periods, water injection method codes, and stress data of the overlying strata. Using the input features of historical periods in the training sample set as the model input, and the porosity of the corresponding voxel in the future preset time window as the prediction label, the time series prediction model is trained so that the model learns the time series mapping relationship between porosity evolution and water level, water injection method and overlying strata stress. The water injection mode codes and overlying strata stress sequences corresponding to the candidate scheduling schemes are input into the trained time-series prediction model to predict the evolution trajectory of the porosity of each voxel under each candidate scheduling scheme. The overall water storage coefficient change curve under each candidate scheme is calculated, and the scheme comparison results are output.

[0096] In this embodiment, the time-varying porosity distribution of each time series is stored in time sequence to construct a porosity time series database. The aim is to systematically record and archive the porosity value of each voxel at different time points. This can be achieved by storing the porosity distribution of each time series (e.g., the porosity value of a three-dimensional mesh) along with the corresponding timestamp in a structured database, such as a relational database, a NoSQL database, or a dedicated time series database.

[0097] Based on this, a training sample set is constructed using porosity time series data of each voxel in the porosity time series database over multiple consecutive historical monitoring periods, water level data for the corresponding time periods, water injection mode codes, and overlying strata stress data. Each sample in the training sample set typically contains a porosity sequence of a voxel over consecutive historical monitoring periods, as well as water level data synchronized with these porosity data, water injection mode codes representing different water injection strategies (e.g., continuous water injection, intermittent water injection, specific flow rate, etc.), and overlying strata stress data (e.g., obtained through sensor measurement or numerical simulation).

[0098] Subsequently, using the input features from historical periods in the training sample set as model input, and the porosity of the corresponding voxels in a future preset time window as the prediction label, the time-series prediction model is trained. This allows the model to learn the temporal mapping relationship between porosity evolution and water level, injection method, and overlying strata stress. This process utilizes machine learning algorithms (e.g., Recurrent Neural Networks (RNNs), Long Short-Term Memory Networks (LSTMs), Gated Recurrent Units (GRUs), Transformer models, or statistical models such as ARIMA) to learn the complex nonlinear relationship between dynamic changes in porosity and external influencing factors. By analyzing patterns in historical data, the model establishes a mapping relationship from input features (historical porosity, water level, injection method, and overlying strata stress) to future porosity, thereby acquiring predictive capabilities.

[0099] Finally, the injection method codes and overlying strata stress sequences corresponding to the candidate scheduling schemes are input into the trained time-series prediction model to predict the evolution trajectory of the porosity of each voxel under each candidate scheduling scheme. The overall storage coefficient variation curve under each candidate scheme is calculated, and the scheme comparison results are output. For each candidate scheduling scheme to be evaluated (e.g., different injection rates, injection cycles, or shutdown times), its corresponding future water level, injection method code, and overlying strata stress sequence are generated and provided as input to the trained model. The model will output the predicted porosity value of each voxel in the future time period. By calculating the predicted porosity of all voxels, the overall storage coefficient variation curve of the entire underground reservoir under different candidate schemes can be calculated over time. Finally, these curves and related performance indicators (such as the total storage capacity and decay rate at a future point in time) will be presented to engineers so that they can intuitively compare the advantages and disadvantages of different schemes and make the optimal decision.

[0100] As an example, a reservoir capacity evolution prediction model is established using accumulated time-series data, employing an LSTM network as the prediction model. Input features include porosity, water level, injection pressure, injection method encoding, and overlying stratum stress from a predetermined number of monitoring periods. The output is the predicted porosity value for the next 30 days. The injection method uses discrete encoding, for example, encoding "continuous injection" as 0, "intermittent injection" as 1, and "stopped injection" as 2. The effective stress of the overlying stratum is calculated based on the mining plan and the thickness of the overlying stratum. All input features are normalized before being input into the network, mapping their values ​​to the [0,1] interval. The network structure is a two-layer LSTM, with 64 units per layer. The training set is based on 12 months of historical monitoring data, using mean squared error as the loss function, the Adam optimizer, a learning rate of 0.001, a batch size of 32, and 200 training epochs. Training is terminated early when the validation loss does not decrease for 10 consecutive epochs.

[0101] This predictive model can accept different boundary condition inputs, thereby assessing the impact of different engineering schemes on the reservoir's water storage performance. For example, changes in the mining plan will alter the stress field of the overlying strata: assuming that mining is carried out on the adjacent coal seam above the reservoir within the next 12 months, the effective stress of the voxel directly above the goaf will change. By inputting the stress change of each voxel as a dynamic boundary condition into the model, the evolution trend of porosity under stress changes can be predicted. Another example is different water replenishment scheduling schemes: Scheme A is "continuous low-flow water replenishment," with a constant daily replenishment of 100 m³. 3 Scheme A corresponds to a stable injection pressure of 0.2 MPa, with the injection mode coded as 0. Scheme B is "intermittent high-flow water replenishment," with 500 m³ of water replenished every 5 days, and no water replenishment at other times. The corresponding injection pressure on replenishment days is 0.5 MPa, and on non-replenishment days it is 0 MPa, with the injection mode coded as 1. Using the two schemes as input sequences for future periods, the model can predict the evolution trajectory of the porosity of each voxel under the two schemes, and then calculate the change curve of the overall water storage coefficient, providing a quantitative basis for selecting the optimal water replenishment scheduling strategy.

[0102] In this embodiment, by constructing a porosity time-series database and training a time-series prediction model using historical data, this application can learn the complex mapping relationship between porosity evolution and water level, water injection method, and overlying strata stress. Furthermore, by simulating different candidate scheduling schemes and predicting the evolution trajectory of porosity of each voxel, the overall water storage coefficient variation curve can be calculated. This provides a scientific basis for the refined management and optimization decision-making of underground water reservoirs in coal mines, effectively avoiding water storage capacity decay caused by improper management, extending the service life of the reservoir, improving resource utilization efficiency, and reducing operational risks.

[0103] In other embodiments, the variation curves of the overall water storage coefficient under each candidate scheme are calculated, and the scheme comparison results are output, including: For each candidate scheduling scheme, the overall water storage coefficient is calculated step by step by the predicted porosity of each voxel at each time step, and the curve of the overall water storage coefficient changing with time under each candidate scheme is generated. The overall water storage coefficient change curves of different candidate schemes are superimposed on the same coordinate graph, and the predicted water storage coefficient values ​​at preset time nodes and the percentage decrease relative to the initial water storage coefficient are marked for engineers to compare and select the best water replenishment scheduling strategy.

[0104] In this embodiment, after obtaining the predicted void ratio values ​​of each voxel under each candidate scheduling scheme within a preset future time window using the time-series prediction model, it is necessary to integrate these discrete prediction data distributed across different voxels and time steps. For each candidate scheduling scheme, at each prediction time step, the predicted void ratio values ​​of all voxels within that time step are calculated. By repeating this process for all prediction time steps, a curve reflecting the dynamic change of the overall water storage coefficient over time under that candidate scheme can be generated. This curve visually demonstrates the evolution trend of the total water storage capacity of the underground reservoir under a specific water replenishment scheduling strategy.

[0105] To facilitate intuitive comparisons by engineers, the overall water storage coefficient variation curves generated for different candidate scheduling schemes are plotted on the same two-dimensional coordinate system. Typically, the horizontal axis represents time, and the vertical axis represents the overall water storage coefficient. Each curve represents a specific candidate scheme and can be distinguished by different colors, line styles, or markers. This overlay display allows engineers to clearly observe the differences in water storage capacity evolution among different schemes, such as which schemes better maintain water storage capacity, which schemes lead to faster decay, and the relative performance of different schemes at specific time points.

[0106] To provide more precise quantitative information, each curve is labeled at key "preset time points" on the overlaid graphs. These labels include the specific predicted value of the overall water storage coefficient for that scheme at that time point. Furthermore, the percentage decrease of this predicted value relative to an "initial water storage coefficient" (e.g., the water storage coefficient at the start of the prediction or a baseline value) is calculated and labeled. The formula for calculating the percentage decrease is typically: ((current predicted value - initial water storage coefficient) / initial water storage coefficient) × 100%. These numerical labels provide engineers with quantitative performance indicators, allowing them not only to see trends but also to accurately understand the specific degree of change in water storage capacity at key time points.

[0107] The core purpose of generating, overlaying, and labeling the overall water storage coefficient variation curves, as well as key data, is to provide engineers with a comprehensive, intuitive, and quantitative decision support tool. Through this visual comparison, engineers can conduct in-depth analysis and evaluation of different water replenishment scheduling strategies based on comprehensive factors such as the actual operational needs of the underground reservoir, safety standards, and economic benefits. For example, they can identify schemes that maintain a high water storage coefficient in the long term, or schemes that exhibit the slowest decay within a specific time period, thereby "selecting the optimal" water replenishment scheduling strategy that best meets current or future needs and optimizes the operation and management of the underground reservoir.

[0108] The following is a specific example to illustrate this: The current porosity of each voxel in the reservoir area With initial porosity The ratio is defined as the water storage capacity retention index. Voxels are chromatographically colored from blue to red according to the water storage capacity retention index from 1 to 0. This can intuitively show the degree of water storage performance degradation in different areas of the entire reservoir. Abnormal degradation areas will be marked prominently on the spectrum.

[0109] Based on the porosity loss rate and predicted trends in various areas of the reservoir, a tiered early warning system is established: a Level III warning is issued when the water storage capacity maintenance index of a certain area is below 0.7; a Level II warning is issued when it is below 0.5; and a Level I warning is issued when it is below 0.3. Combining the simulation results of different scheduling schemes using the prediction model, auxiliary decision support is provided for reservoir management: by comparing the predicted evolution of water storage performance under different water injection schemes using the S5 prediction model, a scheduling scheme that can mitigate the expansion of water storage coefficient decay is recommended; based on the water storage coefficient decay trend of the reservoir area, the overall water storage coefficient of the reservoir area is predicted to change over time. Engineers can manually select when the water storage coefficient is less than 0.1 to reach the decommissioning requirement, thereby estimating the remaining effective service life and providing a scientific basis for the safe and efficient operation of coal mine underground reservoirs throughout their entire life cycle.

[0110] The above technical solution calculates the overall water storage coefficient by interpreting the predicted porosity values ​​of each voxel under each candidate scheduling scheme, and presents this as a visually intuitive time curve, solving the problem of difficulty in making intuitive comparisons using only raw prediction data. Furthermore, by overlaying the overall water storage coefficient change curves of different schemes onto the same coordinate graph and marking the predicted values ​​and attenuation percentages at key time points, engineers can clearly and quantitatively compare the long-term effects of different water replenishment scheduling strategies. This combination of visualization and quantification significantly improves decision-making efficiency and accuracy, helping engineers quickly identify and select the optimal water replenishment scheduling strategy, thereby effectively managing the water storage performance of underground reservoirs and ensuring their safe and stable operation.

[0111] Based on the method for dynamic monitoring of water storage performance of underground coal mine reservoirs based on cosmic ray muon imaging provided in the above embodiments, this application also provides a specific implementation of a device for dynamic monitoring of water storage performance of underground coal mine reservoirs based on cosmic ray muon imaging. Please refer to the following embodiments.

[0112] First see Figure 4 The dynamic monitoring device 400 for the water storage performance of underground coal mine reservoirs based on cosmic ray muon imaging provided in this application embodiment may include: The acquisition module 401 is used to acquire the first muon flux data of the target area; The first determining module 402 is used to determine the time-varying density distribution of each voxel in the target region at each time sequence based on the first muon flux data. The second determining module 403 is used to determine the time-varying porosity of each voxel in each time series based on the time-varying density distribution and the preset object density set. The calculation module 404 is used to calculate the porosity loss rate based on the time-varying porosity of each voxel at each time sequence. Module 405 is used to construct a multi-dimensional feature vector for each voxel based on the porosity loss rate. The multi-dimensional feature vector includes the average loss rate, the slope of the loss rate change trend, the standard deviation of the loss rate fluctuation, the trend fitting determination coefficient, and the cumulative immersion time. Analysis module 406 is used to perform cluster analysis on multidimensional feature vectors to obtain clustering results, so as to monitor the target area based on the clustering results.

[0113] Figure 5 A schematic diagram of the hardware structure of the electronic device provided in an embodiment of this application is shown.

[0114] An electronic device may include a processor 501 and a memory 502 storing computer program instructions.

[0115] Specifically, the processor 501 may include a central processing unit (CPU), an application-specific integrated circuit (ASIC), or one or more integrated circuits that can be configured to implement the embodiments of this application.

[0116] Memory 502 may include mass storage for data or instructions. For example, and not limitingly, memory 502 may include a hard disk drive (HDD), floppy disk drive, flash memory, optical disk, magneto-optical disk, magnetic tape, or Universal Serial Bus (USB) drive, or a combination of two or more of these. In one instance, memory 502 may include removable or non-removable (or fixed) media, or memory 502 may be a non-volatile solid-state memory. Memory 502 may be internal or external to an electronic device.

[0117] In one instance, memory 502 may be read-only memory (ROM). In one instance, the ROM may be a mask-programmed ROM, a programmable ROM (PROM), an erasable PROM (EPROM), an electrically erasable PROM (EEPROM), an electrically rewritable ROM (EAROM), or flash memory, or a combination of two or more of these.

[0118] Memory 502 may include read-only memory (ROM), random access memory (RAM), disk storage media device, optical storage media device, flash memory device, electrical, optical, or other physical / tangible memory storage device. Therefore, generally, memory includes one or more tangible (non-transitory) computer-readable storage media (e.g., memory devices) encoded with software including computer-executable instructions, and when the software is executed (e.g., by one or more processors), it is operable to perform the operations described with reference to the method for dynamic monitoring of water storage performance of underground coal mine reservoirs based on cosmic ray muon imaging according to the first aspect of this disclosure.

[0119] The processor 501 reads and executes computer program instructions stored in the memory 502 to achieve... Figure 1 The embodiment shown is a method for dynamic monitoring of the water storage performance of underground reservoirs in coal mines based on cosmic ray muon imaging.

[0120] In one example, the electronic device may also include a communication interface 503 and a bus 504. For example, Figure 5 As shown, the processor 501, memory 502, and communication interface 503 are connected through bus 504 and complete communication with each other.

[0121] The communication interface 503 is mainly used to realize communication between various modules, devices, units and / or equipment in the embodiments of this application.

[0122] Bus 504 includes hardware, software, or both, that couples components of an electronic device together. For example, and not as a limitation, the bus may include an Accelerated Graphics Port (AGP) or other graphics bus, an Extended Industry Standard Architecture (EISA) bus, a Front Side Bus (FSB), a HyperTransport (HT) interconnect, an Industry Standard Architecture (ISA) bus, an Infinite Bandwidth Interconnect, a Low Pin Count (LPC) bus, a memory bus, a Microchannel Architecture (MCA) bus, a Peripheral Component Interconnect (PCI) bus, a PCI-Express (PCI-X) bus, a Serial Advanced Technology Attachment (SATA) bus, a Video Electronics Standards Association Local (VLB) bus, or other suitable buses, or combinations of two or more of these. Where appropriate, bus 504 may include one or more buses. Although specific buses are described and illustrated in embodiments of this application, this application contemplates any suitable bus or interconnect.

[0123] This electronic device can execute the dynamic monitoring method for the water storage performance of underground coal mine reservoirs based on cosmic ray muon imaging as described in this application embodiment, thereby achieving a combination of... Figures 1-4 The invention describes a method and apparatus for dynamic monitoring of water storage performance of underground reservoirs in coal mines based on cosmic ray muon imaging.

[0124] Furthermore, in conjunction with the dynamic monitoring method for the water storage performance of underground coal mine reservoirs based on cosmic ray muon imaging in the above embodiments, this application embodiment can provide a computer storage medium for implementation. This computer storage medium stores computer program instructions; when these computer program instructions are executed by a processor, they implement any of the dynamic monitoring methods for the water storage performance of underground coal mine reservoirs based on cosmic ray muon imaging in the above embodiments.

[0125] In an optional embodiment, in conjunction with the dynamic monitoring method for the water storage performance of underground coal mine reservoirs based on cosmic ray muon imaging in the above embodiments, this application embodiment can provide a computer program product to implement it. The instructions in the computer program product are executed by the processor of an electronic device, enabling the electronic device to implement any of the dynamic monitoring methods for the water storage performance of underground coal mine reservoirs based on cosmic ray muon imaging in the above embodiments.

[0126] It should be clarified that this application is not limited to the specific configurations and processes described above and shown in the figures. For the sake of brevity, detailed descriptions of known methods are omitted here. In the above embodiments, several specific steps are described and shown as examples. However, the method process of this application is not limited to the specific steps described and shown. Those skilled in the art can make various changes, modifications, and additions, or change the order of steps, after understanding the spirit of this application.

[0127] The functional blocks shown in the above block diagram can be implemented as hardware, software, firmware, or a combination thereof. When implemented in hardware, they can be, for example, electronic circuits, application-specific integrated circuits (ASICs), appropriate firmware, plug-ins, function cards, etc. When implemented in software, the elements of this application are programs or code segments used to perform the required tasks. Programs or code segments can be stored on a machine-readable medium or transmitted over a transmission medium or communication link via data signals carried on a carrier wave. "Machine-readable medium" can include any medium capable of storing or transmitting information. Examples of machine-readable media include electronic circuits, semiconductor memory devices, ROM, flash memory, erasable ROM (EROM), floppy disks, CD-ROMs, optical disks, hard disks, fiber optic media, radio frequency (RF) links, etc. Code segments can be downloaded via computer networks such as the Internet, intranets, etc.

[0128] It should also be noted that the exemplary embodiments mentioned in this application describe methods or systems based on a series of steps or apparatus. However, this application is not limited to the order of the above steps; that is, the steps can be performed in the order mentioned in the embodiments, or in a different order, or several steps can be performed simultaneously.

[0129] The aspects of this disclosure have been described above with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this disclosure. It should be understood that each block in the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, or other programmable data processing apparatus to produce a machine such that these instructions, executable via the processor of the computer or other programmable data processing apparatus, enable the implementation of the functions / actions specified in one or more blocks of the flowchart illustrations and / or block diagrams. Such a processor can be, but is not limited to, a general-purpose processor, a special-purpose processor, a special application processor, or a field-programmable logic circuit. It is also understood that each block in the block diagrams and / or flowcharts, and combinations of blocks in the block diagrams and / or flowcharts, can also be implemented by special-purpose hardware performing the specified functions or actions, or can be implemented by a combination of special-purpose hardware and computer instructions.

[0130] The above description is merely a specific implementation of this application. Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the specific working processes of the systems, modules, and units described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here. It should be understood that the protection scope of this application is not limited thereto. Any person skilled in the art can easily conceive of various equivalent modifications or substitutions within the technical scope disclosed in this application, and these modifications or substitutions should all be covered within the protection scope of this application.

Claims

1. A method for dynamic monitoring of water storage performance of underground reservoirs in coal mines based on cosmic ray muon imaging, characterized in that, include: Obtain the first muon flux data for the target region; Based on the first muon flux data, the time-varying density distribution of each voxel in the target region at each time series is determined; Based on the time-varying density distribution and the preset object density set, the time-varying porosity of each voxel in each time series is determined. The porosity loss rate is calculated based on the time-varying porosity of each voxel at each time series. Based on the porosity loss rate, a multidimensional feature vector is constructed for each voxel. The multidimensional feature vector includes the average loss rate, the slope of the loss rate change trend, the standard deviation of the loss rate fluctuation, the trend fitting determination coefficient, and the cumulative immersion time. Cluster analysis is performed on the multidimensional feature vectors to obtain clustering results, which are then used to monitor the target region.

2. The method according to claim 1, characterized in that, The step of determining the time-varying density distribution of each voxel in the target region at each time series based on the first muon flux data includes: After filtering and denoising the first muon flux data, the incident zenith angle and azimuth angle of each muon are reconstructed according to the hit position of each muon event in the multi-layer structure of the detector. Muons with the same zenith angle and azimuth angle are divided into the same path group. Based on the principle of muon transmission imaging, the muon survival rate of each path group is converted into equivalent mass length. The target region is discretized into a voxel grid. A first geometric length matrix is ​​established based on the geometric length of each muon path in each voxel. The first matrix equation is constructed with the equivalent mass length as the observed value and the voxel density as the variable to be solved. The time-varying density distribution of each voxel at each time series is obtained by solving the first matrix equation using an inversion algorithm.

3. The method according to claim 1, characterized in that, Before acquiring the first muon flux data of the target region, the method further includes: Obtain the spatial distribution of the target area, including the flux distribution of the underground reservoir space and the feasible spatial locations of tunnels and chambers; With detector coverage and reception efficiency as optimization objectives, a multi-objective optimization algorithm is used to solve for the optimal number of muon detector arrays, the placement position of each muon detector, the zenith angle, and the azimuth angle. The acquisition of the first muon flux data of the target region includes: Acquire the first muon flux data of the muon detector array in the target region.

4. The method according to claim 1, characterized in that, Before determining the time-varying porosity of each voxel at each time series based on the time-varying density distribution and the preset object density set, the method further includes: Acquire the flux data of the second muon in the drained state and the flux data of the third muon in the saturated state of the target area, respectively; Based on the second muon flux data and the third muon flux data, the initial density distribution and saturated density distribution of each voxel in each time series under different states of the target region are determined. The initial density distribution is used to characterize the density distribution of the target region under the empty state, and the saturated density distribution is used to characterize the density distribution of the target region under the saturated state. Calculate the initial porosity based on the initial density distribution and the saturated density distribution; The density of the rock solid skeleton is determined based on the initial porosity and the initial density distribution. Obtain geological materials and mining engineering plans of the target area; Based on the rock solid skeleton density, initial porosity, geological materials, and mining engineering plan, a baseline water storage capacity distribution map of the target area is constructed. In order to determine the time-varying porosity of each voxel in each time series according to the baseline water storage capacity distribution map, the time-varying density distribution, and the preset object density set.

5. The method according to claim 4, characterized in that, The step of determining the time-varying porosity of each voxel at each time series based on the baseline water storage capacity distribution map, the time-varying density distribution, and the preset object density set includes: Obtain the water level data of the current time series, and divide the voxel grid into submerged voxels and non-submerged voxels according to the water level data; For voxels in non-submerged areas, the rock solid skeleton density of the voxel is read from the baseline water storage capacity distribution map, and the time-varying porosity is calculated based on the relationship between the time-varying density of the voxel and the rock solid skeleton density. For a voxel in the submerged area, the rock solid skeleton density of the voxel is read from the baseline water storage capacity distribution map, and the time-varying porosity is calculated based on the relationship between the time-varying density of the voxel, the rock solid skeleton density, and the density of water.

6. The method according to claim 1, characterized in that, The construction of a multidimensional feature vector for each voxel based on the porosity loss rate includes: Based on the time-series data of the porosity loss rate of each voxel in each monitoring period, the least squares method is used to perform linear fitting to obtain the slope of the loss rate change trend of each voxel. The slope of the loss rate change trend is used to characterize the acceleration or deceleration trend of the porosity loss of the voxel over time. The arithmetic mean of the time-series data of the loss rate of each voxel is calculated to obtain the average loss rate of each voxel. The standard deviation of the loss rate time series data for each voxel is calculated to obtain the standard deviation of the loss rate fluctuation for each voxel. Calculate the coefficient of determination of the linear fit, which is used to evaluate the reliability of the linear trend.

7. The method according to claim 1, characterized in that, The clustering analysis of the multidimensional feature vectors to obtain the clustering results includes: The multidimensional feature vector of each voxel is standardized to eliminate the influence of differences in size and order of magnitude between different features on cluster analysis. The standardized multidimensional feature vectors are input into the clustering algorithm to divide all voxels into two clusters; Calculate the average loss rate and mean trend slope for each cluster, and determine the region category corresponding to each cluster based on the average loss rate and mean trend slope. The region category includes normal creep region and abnormal decay region.

8. The method according to claim 1, characterized in that, After performing cluster analysis on the multidimensional feature vectors to obtain the clustering results, the method further includes: The time-varying void ratio distributions of each time series are stored in time series to construct a void ratio time series database; A training sample set is constructed using the porosity time series data of each voxel in the porosity time series database for multiple consecutive historical monitoring periods, the water level data for the corresponding time period, the water injection method code, and the stress data of the overlying strata. Using the input features of historical periods in the training sample set as the model input, and using the porosity of the corresponding voxel in a future preset time window as the prediction label, the time series prediction model is trained so that the model learns the time series mapping relationship between porosity evolution and water level, water injection method and overlying strata stress. The water injection mode codes corresponding to the candidate scheduling schemes and the stress sequence of the overlying strata are input into the trained time-series prediction model to predict the evolution trajectory of the porosity of each voxel under each candidate scheduling scheme, calculate the change curve of the overall water storage coefficient under each candidate scheme, and output the scheme comparison results.

9. The method according to claim 8, characterized in that, The calculation yields the overall water storage coefficient variation curves for each candidate scheme, and outputs the scheme comparison results, including: For each candidate scheduling scheme, the overall water storage coefficient is calculated step by step by the predicted porosity of each voxel at each time step, and the curve of the overall water storage coefficient changing with time under each candidate scheme is generated. The overall water storage coefficient change curves of different candidate schemes are superimposed on the same coordinate graph, and the predicted water storage coefficient values ​​at preset time nodes and the percentage decrease relative to the initial water storage coefficient are marked for engineers to compare and select the best water replenishment scheduling strategy.

10. A dynamic monitoring device for the water storage performance of underground reservoirs in coal mines based on cosmic ray muon imaging, characterized in that, The device includes: The acquisition module is used to acquire the first muon flux data of the target region; The first determining module is used to determine the time-varying density distribution of each voxel in the target region at each time sequence based on the first muon flux data. The second determining module is used to determine the time-varying porosity of each voxel in each time series based on the time-varying density distribution and the preset object density set. The calculation module is used to calculate the porosity loss rate based on the time-varying porosity of each voxel at each time series. The construction module is used to construct a multidimensional feature vector for each voxel based on the porosity loss rate. The multidimensional feature vector includes the average loss rate, the slope of the loss rate change trend, the standard deviation of the loss rate fluctuation, the trend fitting determination coefficient, and the cumulative immersion time. The analysis module is used to perform cluster analysis on the multidimensional feature vectors to obtain clustering results, so as to monitor the target region based on the clustering results.

Images