A method and system for coal and rock stress inversion based on multi-material level set topology optimization

Abstract

This application relates to the field of coal mine safety technology, specifically to a method and system for coal and rock stress inversion based on multi-material level set topology optimization. The method includes: setting up an artificial seismic source and multiple detection points; exciting source signals and acquiring seismic wave response signals; constructing a three-dimensional spatial distribution map of the material using at least two level set functions and mapping wave velocity and quality factor; performing forward simulation using the time-fractional-order viscoelastic wave equation, followed by full-waveform inversion to obtain an optimized wave velocity field; constructing a three-parameter coupled model of wave velocity, quality factor, and stress to build a three-dimensional stress field; identifying high-stress concentration areas and anomalous structures, and outputting risk identification results; and integrating the three-dimensional spatial distribution map of the material, the coordinates of the artificial seismic source and detection points, the optimized wave velocity field, the three-dimensional stress field, and the risk identification results into a three-dimensional visualization platform. This application can achieve high-precision wave velocity reconstruction, quantitative stress conversion, and dynamic risk display within coal and rock masses.

Description

Technical Field

[0001] This application relates to the field of coal mine safety technology, specifically to a method and system for coal and rock stress inversion based on multi-material level set topology optimization. Background Technology

[0002] Coal seam rockbursts are sudden and highly destructive dynamic disasters deep within coal and rock masses, posing a serious threat to the safety of underground personnel and mining equipment. As the depth of coal seam mining increases, the stress level of the surrounding rock continues to rise, significantly increasing the probability of rockburst accidents. Therefore, timely acquisition of the three-dimensional stress distribution characteristics within the coal and rock mass is crucial for disaster prevention and control. However, due to the complex underground environment, large-scale direct stress measurement is difficult to achieve.

[0003] In coal and rock media, P-wave velocity is positively correlated with stress level. The three-dimensional P-wave velocity field within the coal and rock mass can indirectly indicate areas of high stress concentration, providing a basis for identifying rockburst risks. Therefore, existing technologies mainly capture rock mass fracturing events through microseismic monitoring and acoustic emission, and combine this with seismic exploration techniques to invert the three-dimensional velocity field within the coal and rock mass, indirectly predicting the stress distribution characteristics within the coal and rock mass. Furthermore, the quality factor Q of the medium, as a characterization index of its stiffness and energy dissipation properties, is used in some existing technologies to infer changes in stress state.

[0004] However, existing methods for predicting stress distribution characteristics within coal and rock masses still have significant limitations: First, the inversion process of three-dimensional wave velocity fields within coal and rock masses generally ignores the viscoelastic dissipation characteristics of the coal and rock masses, leading to distortion in the description of waveform amplitude attenuation and dispersion behavior, thus affecting the accuracy of wave velocity field reconstruction; Second, existing technologies mostly use grid cells to directly invert physical property parameters, resulting in high computational dimensionality and difficulty in clearly characterizing the boundaries of various media, with prominent interface ambiguity issues; Third, existing methods lack a synergistic constraint mechanism between wave velocity, quality factor, and static stress, and have not established a quantitative conversion relationship between dynamic wave field parameters and static stress fields. Summary of the Invention

[0005] To address the problems of waveform distortion, interface ambiguity, and lack of multi-parameter constraints in existing methods for predicting stress within coal and rock masses, which result in insufficient accuracy of wave velocity field reconstruction, unclear characterization of medium boundaries, and inability to quantitatively convert dynamic parameters to static stress fields, this application provides a coal and rock stress inversion method and system based on multi-material level set topology optimization. This method can achieve high-precision wave velocity reconstruction, quantitative stress conversion, and dynamic risk display within coal and rock masses.

[0006] In a first aspect, this application provides a coal and rock stress inversion method based on multi-material level set topology optimization, comprising the following steps:

[0007] S1. Artificial seismic sources and multiple detection points are set up in the area to be explored within the coal and rock mass. Seismic source signals are excited at the artificial seismic sources, and seismic wave response signals at each detection point are collected simultaneously. ;

[0008] in, Spatial location coordinates;

[0009] The wavefield evolution time is calculated from the start of the excitation source signal.

[0010] Number the testing points;

[0011] S2. Obtain empirical geological information of the area to be explored, and based on this empirical geological information, utilize at least two level set functions. Geological structure modeling is performed to obtain the location points of each location in the area to be explored. The level set function values ​​are used to construct a three-dimensional spatial distribution map of the material;

[0012] The level set function values ​​are mapped to wave velocity and quality factor to construct an initial multi-material property parameter field, including an initial three-dimensional wave velocity field and an initial three-dimensional quality factor field.

[0013] in, , The total number of functions in the level set;

[0014] Wave velocity is the propagation speed of P-waves in the medium, and the three-dimensional wave velocity field is the spatial distribution of the propagation speed of P-waves in the medium;

[0015] S3. Based on the source signal and the initial three-dimensional wave velocity field, forward modeling is performed using the time-fractional viscoelastic wave equation, outputting the full-space forward modeling wave field and the forward modeling response signals at each detection point. In the forward modeling process, the short memory fast algorithm is used to solve the fractional derivative terms of the time fractional elastoviscous wave equation. The forward modeling wave field in the whole space is either a sound pressure field or a displacement field.

[0016] S4. Based on seismic wave response signals Forward modeling response signal The full-space forward modeling wave field and the initial wave velocity field are used to perform full-waveform inversion to obtain the optimized wave velocity field;

[0017] S5. A three-parameter coupled model of wave velocity, quality factor and stress is constructed by posterior experimental calibration. Based on the optimized wave velocity field and the initial quality factor field, the stress value of each location point is mapped using the three-parameter coupled model of wave velocity, quality factor and stress to construct a three-dimensional stress field.

[0018] S6. Identify high stress concentration regions and anomalous structures in a three-dimensional stress field, and output risk identification results, including the coordinates of high stress concentration regions and boundary markers of anomalous structures;

[0019] S7. Integrate the material's three-dimensional spatial distribution map, the coordinates of artificial seismic sources and detection points, the optimized wave velocity field, the three-dimensional stress field, and the risk identification results into a three-dimensional visualization platform.

[0020] It should be further noted that in step S2, the empirical geological information includes the distribution of three media: the coal and rock body, the water-bearing weak zone, and the hard interbedded gangue zone.

[0021] Level set function include , ,in:

[0022] and Indicates position The medium at that location is primarily coal and rock;

[0023] and Indicates position The medium at that location is a water-bearing, weak zone;

[0024] Indicates position The medium at that location is a hard coal-containing area.

[0025] It should be further noted that step S2, which involves modeling the geological structure using multi-material level set functions, includes:

[0026] S201. Initialize the level set function ;

[0027] S202. Iteratively update the level set function using a topological evolution method based on the reaction-diffusion equation. The iterative update formula is:

[0028]

[0029] in, Indicates dynamic weighting factor;

[0030] This represents the sensitivity normalization factor;

[0031] Represents the objective function value;

[0032] This represents the diffusion term that controls the smoothness of the boundary.

[0033] Represents the regularization coefficient;

[0034] Indicates the location obtained based on empirical geological information. Place Known value.

[0035] It should be further noted that in step S202, the iterative update will terminate if at least one of the following conditions is met:

[0036] The number of iterations has reached the preset maximum number of iterations.

[0037] The difference between the objective function values ​​updated in two consecutive iterations is less than the preset convergence tolerance;

[0038] In two consecutive iterations, the minimum rate of change of the level set function at the same position is lower than the set threshold.

[0039] It should be further explained that in step S2, the level set values ​​are mapped to wave velocity and quality factor using the SIMP interpolation strategy.

[0040] It should be further noted that in step S3, the time-fractional viscoelastic wave equation is expressed as:

[0041]

[0042] in, Spatial location coordinates;

[0043] For position P-wave propagation speed at that location;

[0044] Represents the forward-modeling wave field across the entire space;

[0045] These are scale parameters related to the relaxation properties of the medium.

[0046] Indicates the order of the Caputo type is The time fractional derivative is defined as:

[0047]

[0048] in , This is the Gamma function.

[0049] It should be further explained that the fractional derivative term of the time-fractional elastoviscous wave equation is specifically solved using the short memory fast algorithm as follows:

[0050] Set memory length as The time step is The current time is Let the fractional derivative term at the current moment be expressed as follows:

[0051]

[0052] in, , ;

[0053] The coefficient of the weighting term is defined as:

[0054] .

[0055] It should be further explained that, .

[0056] It should be further noted that step S4, the full waveform inversion steps, include:

[0057] S401. Using the inversion objective function Calculate the forward modeling response signal and seismic wave response signals Differences between them:

[0058]

[0059] in, To describe the propagation speed of P waves Spatial distribution function;

[0060] Wave field evolution time The total length;

[0061] S402. Calculate the objective function using the adjoint state method. right gradient ;

[0062] S403. Iterative Update The optimized wave velocity field is obtained when the objective function converges.

[0063] It should be further noted that step S403 uses an optimization algorithm for iterative updates. The optimization algorithms include the conjugate gradient method or the L-BFGS method.

[0064] It should be further explained that, in step S5, the specific operation of constructing the three-parameter coupled model of wave velocity, quality factor, and stress through posterior experimental calibration is as follows:

[0065] Standard samples of various media are taken, and different static stresses are applied to each standard sample under laboratory conditions. At the same time, the wave velocity and quality factor of the samples are measured. The empirical curves of wave velocity and quality factor as a function of stress for each medium are fitted by multiple sets of wave velocity and quality factor data under static stress. The formula of the empirical curve is used as the formula of the three-parameter coupling model of wave velocity-quality factor-stress for that medium.

[0066] It should be further explained that in step S6, a method combining threshold criteria and spatial gradient analysis is used to identify high stress concentration regions and anomalous structures in the three-dimensional stress field. The specific operation is as follows:

[0067] Extract the extreme regions in the stress field that exceed a preset multiple of the regional average stress and mark them as high stress concentration regions;

[0068] Analyze the spatial gradients of the wave velocity field and stress field, locate regions where gradient abrupt changes coincide with stress distortion locations, and mark them as anomalous structures.

[0069] Secondly, this application provides a coal and rock stress inversion system based on multi-material level set topology optimization, used to implement the above-mentioned coal and rock stress inversion method, including:

[0070] The data acquisition module is used to set up artificial seismic sources and multiple detection points in the area to be detected within the coal and rock mass. The artificial seismic source generates a source signal and simultaneously acquires the seismic wave response signals of each detection point.

[0071] The geological modeling and initial parameter field construction module is used to acquire empirical geological information of the area to be explored, and based on the empirical geological information, to model the geological structure using at least two level set functions, to obtain the level set function values ​​of each location point in the area to be explored, and to construct a three-dimensional spatial distribution map of materials based on the level set function values; and to map the level set function values ​​to wave velocity and quality factor, and to construct an initial multi-material physical property parameter field, including an initial three-dimensional wave velocity field and an initial three-dimensional quality factor field.

[0072] The forward modeling module is used to perform forward modeling based on the source signal and the initial three-dimensional wave velocity field, using the time fractional viscoelastic wave equation, and outputs the forward modeling wave field in the whole space and the forward modeling response signal of each detection point.

[0073] The full waveform inversion module is used to perform full waveform inversion based on seismic wave response signal, forward modeling response signal, full-space forward modeling wave field and initial wave velocity field to obtain optimized wave velocity field;

[0074] The stress field construction module is used to construct a three-parameter coupled model of wave velocity, quality factor and stress by posterior experimental calibration. Based on the optimized wave velocity field and the initial quality factor field, the three-parameter coupled model of wave velocity, quality factor and stress is used to map the stress value at each location point to construct a three-dimensional stress field.

[0075] The risk identification module is used to identify high stress concentration areas and anomalous structures in a three-dimensional stress field and output the risk identification results.

[0076] The 3D visualization integration module is used to integrate the material's 3D spatial distribution map, the coordinates of artificial seismic sources and detection points, the optimized wave velocity field, the 3D stress field, and the risk identification results into the 3D visualization platform.

[0077] Thirdly, this application provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the steps of the above-described coal and rock stress inversion method.

[0078] Fourthly, this application provides a storage medium storing a computer program, which, when executed by a processor, implements the steps of the above-described coal and rock stress inversion method.

[0079] As can be seen from the above technical solutions, this application has the following advantages:

[0080] 1. This application solves the problems of high calculation dimensionality and blurred interface in the existing grid inversion method by using at least two level set functions to model geological structures, construct a three-dimensional spatial distribution map of materials, map wave velocity and quality factor, and construct an initial multi-material physical property parameter field. It achieves clear division of complex geological structures and continuous mapping of physical property parameters.

[0081] 2. This application utilizes the time-fractional viscoelastic wave equation for forward modeling. During the process, a short-memory fast algorithm is used to solve the fractional derivative terms of the time-fractional viscoelastic wave equation, which solves the waveform distortion problem of traditional elastic models, realizes high-fidelity simulation of seismic wave propagation characteristics, and provides a reliable data foundation for wave velocity inversion.

[0082] 3. This application constructs a three-parameter coupled model of wave velocity, quality factor, and stress by means of a post-hoc experimental calibration, maps the stress value at each location point, and constructs a three-dimensional stress field, which solves the problem of missing multi-parameter collaborative constraints and realizes the quantitative conversion from dynamic parameters to a three-dimensional static stress field.

[0083] 4. This application identifies high-stress areas and abnormal structures based on a three-dimensional stress field and integrates the results into a visualization platform, which solves the problem of insufficient early warning information output and realizes accurate marking and three-dimensional dynamic display of risk sources. Attached Figure Description

[0084] To more clearly illustrate the technical solution of this application, the accompanying drawings used in the description will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0085] Figure 1 This is a flowchart of a coal and rock stress inversion method based on multi-material level set topology optimization in one embodiment of this application.

[0086] Figure 2 This is a schematic block diagram of a coal and rock stress inversion system based on multi-material level set topology optimization in one embodiment of this application.

[0087] Figure 3 This is a schematic diagram of the hardware structure of an electronic device in one embodiment of this application. Detailed Implementation

[0088] To make the purpose, features, and advantages of this application more apparent and understandable, specific embodiments and accompanying drawings will be used to clearly and completely describe the technical solution protected by this application. Obviously, the embodiments described below are only some embodiments of this application, and not all embodiments. Based on the embodiments in this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0089] The coal and rock stress inversion method involved in this application will be described in detail below. Specific details such as particular system structures and technologies are presented for illustrative purposes and not for limitation, in order to provide a thorough understanding of the embodiments of this application. However, those skilled in the art will understand that this application can also be implemented in other embodiments without these specific details.

[0090] In the coal and rock stress inversion method involved in this application, the term "comprising" indicates the presence of the described feature, whole, step, operation, element, and / or component, but does not exclude the presence or addition of one or more other features, wholes, steps, operations, elements, components, and / or sets thereof. The terms "comprising," "including," "having," and variations thereof all mean "including but not limited to," unless otherwise specifically emphasized.

[0091] To facilitate a clear description of the technical solutions of this application, the terms "first" and "second" are used to distinguish identical or similar items with essentially the same function and effect. Those skilled in the art will understand that the terms "first" and "second" do not limit the quantity or execution order, and that the terms "first" and "second" do not necessarily imply that they are different.

[0092] The terms "one embodiment" or "some embodiments" used in this application mean that one or more embodiments of this application include the specific features, structures, or characteristics described in that embodiment. Therefore, the phrases "in one embodiment," "in some embodiments," "in other embodiments," "in still other embodiments," etc., appearing in different parts of this application do not necessarily refer to the same embodiment, but rather mean "one or more, but not all, embodiments," unless otherwise specifically emphasized.

[0093] The technical solutions in the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings.

[0094] The coal and rock stress inversion method provided in this application embodiment is executed by a computer device, and correspondingly, the coal and rock stress inversion system based on multi-material level set topology optimization runs in the computer device.

[0095] Figure 1 This is a flowchart of a coal and rock stress inversion method based on multi-material level set topology optimization, according to an embodiment of this application. Figure 1 The implementing entity can be a coal and rock stress inversion system. Depending on different needs, the order of the steps in this flowchart can be changed, and some steps can be omitted.

[0096] like Figure 1 As shown, the coal and rock stress inversion method based on multi-material level set topology optimization includes:

[0097] Step S1: Set up an artificial seismic source and multiple detection points in the area to be detected within the coal and rock mass. Excite source signals at the artificial seismic source and simultaneously collect seismic wave response signals from each detection point. ;

[0098] in, Spatial location coordinates;

[0099] The wavefield evolution time is calculated from the start of the excitation source signal.

[0100] This is the testing point number.

[0101] By deploying artificial seismic sources and multiple receivers in the area to be explored and simultaneously acquiring seismic wave response signals, comprehensive and high signal-to-noise ratio observation data on the dynamic response of underground structures were obtained, providing a real and reliable input information source for subsequent inversion.

[0102] Step S2: Obtain empirical geological information of the area to be explored. Based on the empirical geological information, utilize at least two level set functions. Geological structure modeling is performed to obtain the location points of each location in the area to be explored. The level set function values ​​are used to construct a three-dimensional spatial distribution map of the material;

[0103] The level set function values ​​are mapped to wave velocity and quality factor to construct an initial multi-material property parameter field, including an initial three-dimensional wave velocity field and an initial three-dimensional quality factor field.

[0104] in, , The total number of functions in the level set;

[0105] Wave velocity is the propagation speed of P-waves in the medium, and the three-dimensional wave velocity field is the spatial distribution of the propagation speed of P-waves in the medium.

[0106] By using multi-material level set functions based on empirical geological information to model geological structures and mapping level set values ​​to wave velocity and quality factor, a reasonable initial three-dimensional physical property parameter field (wave velocity field and quality factor field) for complex coal and rock strata (containing multiple media) is realized, providing an initial physical model that conforms to the geological background for the forward modeling of wave equations.

[0107] In some specific embodiments, empirical geological information includes the distribution of three media: the coal and rock body, the water-bearing weak zone, and the hard interbedded gangue zone;

[0108] Level set function include , ,in:

[0109] and Indicates position The medium at that location is primarily coal and rock;

[0110] and Indicates position The medium at that location is a water-bearing, weak zone;

[0111] Indicates position The medium at that location is a hard coal-containing area.

[0112] By limiting empirical geological information to include three key media types and defining the combination judgment rules of horizontal set functions, a clear distinction and accurate geometric characterization of the coal and rock body, water-bearing weak zone and hard interbedded gangue zone are achieved. This ensures the structural authenticity of the initial geological model and the rationality of the subsequent assigned physical property parameters (wave velocity, quality factor), providing reliable geological background constraints for stress inversion.

[0113] In some specific embodiments, the steps for geological structure modeling using multi-material level set functions include:

[0114] S201. Initialize the level set function ;

[0115] S202. Iteratively update the level set function using a topological evolution method based on the reaction-diffusion equation. The iterative update formula is:

[0116]

[0117] in, Indicates dynamic weighting factor;

[0118] This represents the sensitivity normalization factor;

[0119] Represents the objective function value;

[0120] This represents the diffusion term that controls the smoothness of the boundary.

[0121] Represents the regularization coefficient;

[0122] Indicates the location obtained based on empirical geological information. Place Known value.

[0123] By employing a topological evolution method based on reaction-diffusion equations to iteratively update the level set function, adaptive optimization and smooth boundary transition of the geological structure model are achieved. This can improve the model's ability to characterize complex structures (such as thin interlayers and irregular interfaces) while matching prior geological information, laying an accurate foundation for subsequent initialization of physical property parameter fields.

[0124] In some specific embodiments, in step S202, the iterative update is terminated if at least one of the following conditions is met:

[0125] The number of iterations has reached the preset maximum number of iterations.

[0126] The difference between the objective function values ​​updated in two consecutive iterations is less than the preset convergence tolerance;

[0127] In two consecutive iterations, the minimum rate of change of the level set function at the same position is lower than the set threshold.

[0128] By setting clear iteration termination conditions, including the maximum number of iterations, the tolerance for changes in the objective function, and the threshold for the rate of change of the level set function, intelligent control and computational resource optimization of the level set topology evolution process are achieved. This can effectively prevent the algorithm from getting stuck in invalid iterations or premature convergence, and ensure that geological modeling achieves a balance between efficiency and accuracy.

[0129] In some specific embodiments, the SIMP interpolation strategy is used to map the level set values ​​to wave velocity and quality factor.

[0130] By using an interpolation strategy to map the level set function values ​​to physical property parameters, including wave velocity and quality factor, a smooth transition from a discrete geometric model to a continuous physical parameter field is achieved. This avoids non-physical jumps in physical property parameters at material interfaces, ensuring the physical consistency and numerical stability of the initial wave velocity field and quality factor field, and supporting subsequent forward and inverse models.

[0131] Step S3: Based on the source signal and the initial three-dimensional wave velocity field, perform forward modeling using the time-fractional viscoelastic wave equation, and output the full-space forward modeling wave field and the forward modeling response signals of each detection point. In the forward modeling process, the short memory fast algorithm is used to solve the fractional derivative terms of the time fractional elastoviscous wave equation. The forward modeling wave field in the whole space is either a sound pressure field or a displacement field.

[0132] By performing forward modeling based on the initial wave velocity field using the time fractional viscoelastic wave equation and employing a short-memory fast algorithm to solve the fractional derivative terms, accurate numerical simulation of the propagation process of seismic wave fields with attenuation and dispersion effects was achieved. The theoretical simulation response signals and spatial wave field evolution at each receiver point were obtained, providing forward modeling prediction data for inversion.

[0133] In some specific embodiments, the time-fractional viscoelastic wave equation is expressed as:

[0134]

[0135] in, Spatial location coordinates;

[0136] For position P-wave propagation speed at that location;

[0137] Represents the forward-modeling wave field across the entire space;

[0138] These are scale parameters related to the relaxation properties of the medium.

[0139] Indicates the order of the Caputo type is The time fractional derivative is defined as:

[0140]

[0141] in , This is the Gamma function.

[0142] By employing a viscoelastic wave equation containing a time fractional derivative term to describe wave propagation, a precise mathematical characterization of the dispersion and attenuation effects of seismic waves in coal and rock formations is achieved. This significantly improves the fitting accuracy of forward modeling to the actual observed wave field and provides a high-fidelity forward modeling engine for full waveform inversion.

[0143] In some specific embodiments, the fractional derivative term of the time-fractional elastoviscous wave equation is solved using the short memory fast algorithm as follows:

[0144] Set memory length as The time step is The current time is Let the fractional derivative term at the current moment be expressed as follows:

[0145]

[0146] in, , ;

[0147] The coefficient of the weighting term is defined as:

[0148] .

[0149] By employing a short-memory fast algorithm to approximate the fractional derivative terms (using only historical wave field values ​​with finite time steps and specific weighting coefficients), the computational complexity of solving the fractional viscoelastic wave equation is significantly reduced. This can greatly improve the computational efficiency of forward modeling while ensuring accuracy, making full waveform inversion feasible on an engineering scale.

[0150] In some specific embodiments, .

[0151] By limiting the upper limit of the short memory fast algorithm, further optimization and control of computational accuracy and efficiency are achieved. This effectively suppresses the cumulative error caused by excessive memory length and controls memory consumption, ensuring the numerical stability and efficiency of forward modeling.

[0152] Step S4, based on seismic wave response signal Forward modeling response signal The full-space forward modeling wave field and the initial wave velocity field are used to perform full-waveform inversion to obtain the optimized wave velocity field.

[0153] By performing full waveform inversion based on measured response signals, forward modeling response signals, spatial wave field, and initial wave velocity field, iterative optimization and correction of the initial wave velocity model were achieved, resulting in a high-precision three-dimensional optimized wave velocity field that better matches the actual observation data.

[0154] In some specific embodiments, the full waveform inversion steps include:

[0155] S401. Using the inversion objective function Calculate the forward modeling response signal and seismic wave response signals Differences between them:

[0156]

[0157] in, To describe the propagation speed of P waves Spatial distribution function;

[0158] Wave field evolution time The total length;

[0159] S402. Calculate the objective function using the adjoint state method. right gradient ;

[0160] S403. Iterative Update The optimized wave velocity field is obtained when the objective function converges.

[0161] By defining an inversion objective function based on the sum of squared waveform residuals and applying a gradient optimization algorithm for iterative updates, the optimal estimation of the spatial distribution function of wave velocity in the least squares sense is achieved. This minimizes the difference between the observed wave field and the simulated wave field, thereby obtaining a high-resolution optimized wave velocity field that conforms to actual geological and physical laws.

[0162] In some specific embodiments, step S403 uses an optimization algorithm to iteratively update. The optimization algorithms include the conjugate gradient method or the L-BFGS method.

[0163] By limiting the use of a specific gradient optimization algorithm for wave velocity field updates, the objective function is reduced efficiently and stably during the iteration process. This accelerates the convergence speed of full waveform inversion and enhances the robustness of the algorithm under complex models, ensuring the acquisition of a reliable optimized wave velocity field.

[0164] Step S5: Construct a three-parameter coupled model of wave velocity, quality factor, and stress using a post-hoc experimental calibration method. Based on the optimized wave velocity field and the initial quality factor field, map the stress value at each location point using the three-parameter coupled model of wave velocity, quality factor, and stress to construct a three-dimensional stress field.

[0165] By utilizing a coupled wave velocity-quality factor-stress model calibrated by a post-hoc experiment and combining optimized wave velocity field and initial quality factor field to map stress values, the dynamic physical property parameters obtained from seismic wave inversion were transformed into a static three-dimensional absolute stress field, providing direct mechanical state information for risk assessment.

[0166] In some specific embodiments, the specific operation of constructing the wave velocity-quality factor-stress three-parameter coupled model through posterior experimental calibration is as follows:

[0167] Standard samples of various media are taken, and different static stresses are applied to each standard sample under laboratory conditions. At the same time, the wave velocity and quality factor of the samples are measured. The empirical curves of wave velocity and quality factor as a function of stress for each medium are fitted by multiple sets of wave velocity and quality factor data under static stress. The formula of the empirical curve is used as the formula of the three-parameter coupling model of wave velocity-quality factor-stress for that medium.

[0168] By establishing an empirical coupling model between wave velocity, quality factor and stress in different media under static stress using laboratory calibration methods, we can quantitatively map the absolute stress value inside coal and rock mass using wave velocity field and quality factor information obtained from field inversion. This can transform the seismic wave inversion results into a three-dimensional stress field that can be directly used for engineering risk assessment.

[0169] Step S6: Identify high stress concentration areas and anomalous structures in the three-dimensional stress field, and output the risk identification results, including the coordinates of high stress concentration areas and the boundary markers of anomalous structures.

[0170] By using quantitative methods in a three-dimensional stress field to identify high stress concentration areas and anomalous structures, the system achieves automated location and boundary delineation of potential hazards for coal and rock dynamic disasters (such as rockbursts), and outputs specific risk identification results to guide safety prevention and control.

[0171] In some specific embodiments, a method combining threshold criteria and spatial gradient analysis is used to identify high stress concentration regions and anomalous structures in a three-dimensional stress field. The specific operation is as follows:

[0172] Extract the extreme regions in the stress field that exceed a preset multiple of the regional average stress and mark them as high stress concentration regions;

[0173] Analyze the spatial gradients of the wave velocity field and stress field, locate regions where gradient abrupt changes coincide with stress distortion locations, and mark them as anomalous structures.

[0174] By combining stress threshold criteria with physical property field gradient analysis to identify risk areas, the system achieves automated and quantitative location and boundary extraction of high stress concentration areas (corresponding to abnormally high stress values) and potential abnormal structures (corresponding to areas of abrupt changes in stress / wave velocity gradient). This can accurately identify potential hazards such as rockbursts and guide the formulation of prevention and control measures.

[0175] Step S7: Integrate the material's three-dimensional spatial distribution map, the coordinates of the artificial seismic source and detection points, the optimized wave velocity field, the three-dimensional stress field, and the risk identification results into the three-dimensional visualization platform.

[0176] By integrating geological models, observation systems, wave velocity fields, stress fields, and risk identification results obtained through inversion into a three-dimensional visualization platform, a comprehensive and intuitive display of the stress state and spatial distribution of risks in coal and rock is achieved, providing an intuitive and efficient interactive analysis environment for geological interpretation, engineering decision-making, and disaster early warning.

[0177] In some specific embodiments, the 3D visualization platform is based on the 3D spatial coordinates inside the coal and rock, and displays the coordinates of artificial seismic sources and detection points, optimized wave velocity field, 3D stress field, risk identification results, mine roadway layout, borehole measuring points, etc. in a 3D overlay. The platform software supports functions such as multi-window display, real-time data input, background calculation and result query.

[0178] For example, the main interface can simultaneously display the three-dimensional structural model of the coal seam, the location and frequency of real-time microseismic events, the readings of stress sensors in each monitoring borehole, and a comprehensive risk assessment index. Users can intuitively view the relative location and range of high-stress zones in the roadway by rotating and zooming the three-dimensional view.

[0179] In some specific embodiments, the three-dimensional visualization platform connects to the online monitoring data stream. When the microseismic monitoring detects abnormally frequent energy events or when the stress at a certain monitoring point approaches the warning threshold, the platform will automatically update and optimize the wave velocity field and the three-dimensional stress field, and highlight the dangerous area in the interface. At the same time, the platform assesses the impact hazard level based on comprehensive indicators and gives the real-time risk level of each zone in color or numerical form at the bottom of the interface.

[0180] The 3D visualization platform can provide timely forecasts and early warnings of rockburst hazards caused by high stress concentration in actual production, significantly improving the scientific nature and effectiveness of mine safety management. With the help of the 3D visualization platform, mine workers can monitor the evolution trend of coal seam stress and potential rockburst hazards around the clock, realizing a shift from passive rescue to proactive early warning.

[0181] In one specific embodiment, the steps of the coal and rock stress inversion method based on multi-material level set topology optimization include:

[0182] Step S1: Set up an artificial seismic source and multiple detection points in the area to be detected within the coal and rock mass. Excite source signals at the artificial seismic source and simultaneously collect seismic wave response signals from each detection point. ;

[0183] in, Spatial location coordinates;

[0184] The wavefield evolution time is calculated from the start of the excitation source signal.

[0185] This is the testing point number.

[0186] Step S2: Obtain empirical geological information of the area to be explored. Based on the empirical geological information, utilize at least two level set functions. Geological structure modeling is performed to obtain the location points of each location in the area to be explored. The level set function values ​​are used to construct a three-dimensional spatial distribution map of the material;

[0187] Among them, empirical geological information includes the distribution of three media: the main coal and rock mass, the water-bearing weak zone, and the hard interbedded gangue zone;

[0188] Level set function include , ,in:

[0189] and Indicates position The medium at that location is primarily coal and rock;

[0190] and Indicates position The medium at that location is a water-bearing, weak zone;

[0191] Indicates position The medium at that location is a hard, interbedded area of ​​gangue;

[0192] The steps for geological structure modeling using multi-material level set functions include:

[0193] S201. Initialize the level set function ;

[0194] S202. Iteratively update the level set function using a topological evolution method based on the reaction-diffusion equation. The iterative update formula is:

[0195]

[0196] in, Indicates dynamic weighting factor;

[0197] This represents the sensitivity normalization factor;

[0198] Represents the objective function value;

[0199] This represents the diffusion term that controls the smoothness of the boundary.

[0200] Represents the regularization coefficient;

[0201] Indicates the location obtained based on empirical geological information. Place Known value;

[0202] The iterative update will terminate if at least one of the following conditions is met:

[0203] The number of iterations has reached the preset maximum number of iterations.

[0204] The difference between the objective function values ​​updated in two consecutive iterations is less than the preset convergence tolerance;

[0205] In two consecutive iterations, the minimum rate of change of the level set function at the same position is lower than a set threshold;

[0206] The SIMP interpolation strategy is used to map the level set function values ​​to wave velocity and quality factor to construct an initial multi-material property parameter field, including an initial three-dimensional wave velocity field and an initial three-dimensional quality factor field.

[0207] Wave velocity is the propagation speed of P-waves in the medium, and the three-dimensional wave velocity field is the spatial distribution of the propagation speed of P-waves in the medium.

[0208] Step S3: Based on the source signal and the initial three-dimensional wave velocity field, perform forward modeling using the time-fractional viscoelastic wave equation, and output the full-space forward modeling wave field and the forward modeling response signals of each detection point. In the forward modeling process, the short memory fast algorithm is used to solve the fractional derivative terms of the time fractional elastoviscous wave equation. The forward modeling wave field in the whole space is either a sound pressure field or a displacement field.

[0209] The time-fractional order viscoelastic wave equation is expressed as:

[0210]

[0211] in, Spatial location coordinates;

[0212] For position P-wave propagation speed at that location;

[0213] Represents the forward-modeling wave field across the entire space;

[0214] These are scale parameters related to the relaxation properties of the medium.

[0215] Indicates the order of the Caputo type is The time fractional derivative is defined as:

[0216]

[0217] in , It is the Gamma function;

[0218] The fractional derivative term of the time-fractional elastoviscous wave equation is solved using the short memory fast algorithm as follows:

[0219] Set memory length as The time step is The current time is Let the fractional derivative term at the current moment be expressed as follows:

[0220]

[0221] in, , ;

[0222] The coefficient of the weighting term is defined as:

[0223] , .

[0224] Step S4, based on seismic wave response signal Forward modeling response signal The full-space forward modeling wave field and the initial wave velocity field are used to perform full-waveform inversion to obtain the optimized wave velocity field;

[0225] The steps of full waveform inversion include:

[0226] S401. Using the inversion objective function Calculate the forward modeling response signal and seismic wave response signals Differences between them:

[0227]

[0228] in, To describe the propagation speed of P waves Spatial distribution function;

[0229] Wave field evolution time The total length;

[0230] S402. Calculate the objective function using the adjoint state method. right gradient ;

[0231] S403. Iteratively update using an optimization algorithm. Once the objective function converges, the optimized wave velocity field is obtained. Optimization algorithms include the conjugate gradient method or the L-BFGS method.

[0232] Step S5: Construct a three-parameter coupled model of wave velocity, quality factor, and stress by posterior experimental calibration. Based on the optimized wave velocity field and the initial quality factor field, use the three-parameter coupled model of wave velocity, quality factor, and stress to map the stress value at each location point and construct a three-dimensional stress field.

[0233] The specific steps for constructing the wave velocity-quality factor-stress three-parameter coupled model using the post-experimental calibration method are as follows:

[0234] Standard samples of various media are taken, and different static stresses are applied to each standard sample under laboratory conditions. At the same time, the wave velocity and quality factor of the samples are measured. The empirical curves of wave velocity and quality factor as a function of stress for each medium are fitted by multiple sets of wave velocity and quality factor data under static stress. The formula of the empirical curve is used as the formula of the three-parameter coupling model of wave velocity-quality factor-stress for that medium.

[0235] Step S6 involves identifying high stress concentration regions and anomalous structures in the three-dimensional stress field. The specific steps are as follows:

[0236] Extract the extreme regions in the stress field that exceed a preset multiple of the regional average stress and mark them as high stress concentration regions;

[0237] Analyze the spatial gradients of the wave velocity field and stress field, locate regions where gradient abrupt changes coincide with stress distortion locations, and mark them as anomalous structures.

[0238] Output risk identification results, including coordinates of high stress concentration areas and boundary markers of anomalous structures.

[0239] Step S7: Integrate the material's three-dimensional spatial distribution map, the coordinates of the artificial seismic source and detection points, the optimized wave velocity field, the three-dimensional stress field, and the risk identification results into the three-dimensional visualization platform.

[0240] The following are embodiments of the coal and rock stress inversion system based on multi-material level set topology optimization provided in this application. This coal and rock stress inversion system based on multi-material level set topology optimization belongs to the same inventive concept as the coal and rock stress inversion methods in the above embodiments. For details not described in the embodiments of the coal and rock stress inversion system, please refer to the embodiments of the coal and rock stress inversion methods based on multi-material level set topology optimization described above.

[0241] like Figure 2 As shown, the coal and rock stress inversion system based on multi-material level set topology optimization includes:

[0242] The data acquisition module is used to set up artificial seismic sources and multiple detection points in the area to be detected within the coal and rock mass. The artificial seismic source generates a source signal and simultaneously acquires the seismic wave response signals of each detection point.

[0243] The geological modeling and initial parameter field construction module is used to acquire empirical geological information of the area to be explored, and based on the empirical geological information, to model the geological structure using at least two level set functions, to obtain the level set function values ​​of each location point in the area to be explored, and to construct a three-dimensional spatial distribution map of materials based on the level set function values; and to map the level set function values ​​to wave velocity and quality factor, and to construct an initial multi-material physical property parameter field, including an initial three-dimensional wave velocity field and an initial three-dimensional quality factor field.

[0244] The forward modeling module is used to perform forward modeling based on the source signal and the initial three-dimensional wave velocity field, using the time fractional viscoelastic wave equation, and outputs the forward modeling wave field in the whole space and the forward modeling response signal of each detection point.

[0245] The full waveform inversion module is used to perform full waveform inversion based on seismic wave response signal, forward modeling response signal, full-space forward modeling wave field and initial wave velocity field to obtain optimized wave velocity field;

[0246] The stress field construction module is used to construct a three-parameter coupled model of wave velocity, quality factor and stress by posterior experimental calibration. Based on the optimized wave velocity field and the initial quality factor field, the three-parameter coupled model of wave velocity, quality factor and stress is used to map the stress value at each location point to construct a three-dimensional stress field.

[0247] The risk identification module is used to identify high stress concentration areas and anomalous structures in a three-dimensional stress field and output the risk identification results.

[0248] The 3D visualization integration module is used to integrate the material's 3D spatial distribution map, the coordinates of artificial seismic sources and detection points, the optimized wave velocity field, the 3D stress field, and the risk identification results into the 3D visualization platform.

[0249] The coal and rock stress inversion system in this embodiment is used to implement a coal and rock stress inversion method based on multi-material level set topology optimization.

[0250] This application also provides an electronic device for implementing the various embodiments of this application. Figure 3 To illustrate the hardware structure of an electronic device according to various embodiments of this application, as shown in the following diagram... Figure 3 As shown, the electronic device includes a memory, a processor, and a computer program stored in the memory and capable of running on the processor.

[0251] Those skilled in the art will understand that the electronic device structure involved in the embodiments of this application does not constitute a limitation on the electronic device. The electronic device may include more or fewer components than shown in the figure, or combine certain components, or have different component arrangements.

[0252] In embodiments of this application, electronic devices include, but are not limited to, laptop computers, desktop computers, workstations, personal digital assistants, servers, blade servers, mainframe computers, and other suitable computers. Electronic devices may also represent various forms of mobile devices and other similar computing devices. The components shown herein, their connections and relationships, and their functions are merely examples and are not intended to limit the implementation of the embodiments of this application described and / or claimed herein.

[0253] In this application embodiment, the processor can be implemented using at least one of an Application-Specific Integrated Circuit (ASIC), a Digital Signal Processor (DSP), a Digital Signal Processing Device (DSPD), a processor, a controller, a microcontroller, a microprocessor, or an electronic unit designed to perform the functions described herein. In some cases, such implementations can be implemented within a controller. For software implementations, implementations such as processes or functions can be implemented with separate software modules that allow the performance of at least one function or operation. The software code can be implemented by a software application (or program) written in any suitable programming language, and the software code can be stored in memory and executed by the controller.

[0254] In addition, the electronic device includes some functional modules not shown, which will not be described in detail here.

[0255] Those skilled in the art will understand that the various aspects of the electronic device provided in this application can be implemented as a system, method, or program product. Therefore, the various aspects of this application can be specifically implemented in the following forms: a completely hardware implementation, a completely software implementation (including firmware, microcode, etc.), or a combination of hardware and software aspects, collectively referred to herein as a "circuit," "module," or "system."

[0256] This application also provides a storage medium storing a program product capable of implementing a coal and rock stress inversion method based on multi-material level set topology optimization. In some possible embodiments, various aspects of this application may also be implemented as a program product comprising program code that, when run on a terminal device, causes the terminal device to perform the steps described in the foregoing "Exemplary Methods" section of this specification according to various exemplary embodiments of this application.

[0257] The storage medium may be any combination of one or more readable media. A readable medium may be a readable signal medium or a readable storage medium. A readable storage medium may be, for example,, but not limited to, an electrical, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any combination thereof. More specific examples (a non-exhaustive list) of readable storage media include: electrical connections having one or more wires, portable disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), optical fiber, portable compact disk read-only memory (CD-ROM), optical storage devices, magnetic storage devices, or any suitable combination thereof.

[0258] The above description of the disclosed embodiments enables those skilled in the art to make or use this application. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of this application. Therefore, this application is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for coal and rock stress inversion based on multi-material level set topology optimization, characterized in that, include: S1. Artificial seismic sources and multiple detection points are set up in the area to be explored within the coal and rock mass. Seismic source signals are excited at the artificial seismic sources, and seismic wave response signals at each detection point are collected simultaneously. ; in, Spatial location coordinates; The wavefield evolution time is calculated from the start of the excitation source signal. Number the testing points; S2. Obtain empirical geological information of the area to be explored, and based on this empirical geological information, utilize at least two level set functions. Geological structure modeling is performed to obtain the location points of each location in the area to be explored. The level set function values ​​are used to construct a three-dimensional spatial distribution map of the material; The level set function values ​​are mapped to wave velocity and quality factor to construct an initial multi-material property parameter field, including an initial three-dimensional wave velocity field and an initial three-dimensional quality factor field. in, , The total number of functions in the level set; Wave velocity is the propagation speed of P-waves in the medium, and the three-dimensional wave velocity field is the spatial distribution of the propagation speed of P-waves in the medium; S3. Based on the source signal and the initial three-dimensional wave velocity field, forward modeling is performed using the time-fractional viscoelastic wave equation, outputting the full-space forward modeling wave field and the forward modeling response signals at each detection point. In the forward modeling process, the short memory fast algorithm is used to solve the fractional derivative terms of the time fractional elastoviscous wave equation. The forward modeling wave field in the whole space is either a sound pressure field or a displacement field. S4. Based on seismic wave response signals Forward modeling response signal The full-space forward modeling wave field and the initial wave velocity field are used to perform full-waveform inversion to obtain the optimized wave velocity field; S5. A three-parameter coupled model of wave velocity, quality factor and stress is constructed by posterior experimental calibration. Based on the optimized wave velocity field and the initial quality factor field, the stress value of each location point is mapped using the three-parameter coupled model of wave velocity, quality factor and stress to construct a three-dimensional stress field. S6. Identify high stress concentration regions and anomalous structures in a three-dimensional stress field, and output risk identification results, including the coordinates of high stress concentration regions and boundary markers of anomalous structures; S7. Integrate the material's three-dimensional spatial distribution map, the coordinates of artificial seismic sources and detection points, the optimized wave velocity field, the three-dimensional stress field, and the risk identification results into a three-dimensional visualization platform.

2. The coal and rock stress inversion method as described in claim 1, characterized in that, Step S2, which involves modeling the geological structure using multi-material level set functions, includes: S201. Initialize the level set function ; S202. Iteratively update the level set function using a topological evolution method based on the reaction-diffusion equation. The iterative update formula is: in, Indicates dynamic weighting factor; This represents the sensitivity normalization factor; Represents the objective function value; This represents the diffusion term that controls the smoothness of the boundary. Represents the regularization coefficient; Indicates the location obtained based on empirical geological information. Place Known value.

3. The coal and rock stress inversion method as described in claim 1, characterized in that, In step S202, the iterative update terminates if at least one of the following conditions is met: The number of iterations has reached the preset maximum number of iterations. The difference between the objective function values ​​updated in two consecutive iterations is less than the preset convergence tolerance; In two consecutive iterations, the minimum rate of change of the level set function at the same position is lower than the set threshold.

4. The coal and rock stress inversion method as described in claim 1, characterized in that, In step S2, the level set values ​​are mapped to wave velocity and quality factor using the SIMP interpolation strategy.

5. The coal and rock stress inversion method as described in claim 1, characterized in that, In step S3, the time-fractional viscoelastic wave equation is expressed as: in, Spatial location coordinates; For position P-wave propagation speed at that location; Represents the forward-modeling wave field across the entire space; These are scale parameters related to the relaxation properties of the medium. Indicates the order of the Caputo type is The time fractional derivative is defined as: in , This is the Gamma function.

6. The coal and rock stress inversion method as described in claim 1, characterized in that, The fractional derivative term of the time-fractional elastoviscous wave equation is solved using the short memory fast algorithm as follows: Set memory length as The time step is The current time is Let the fractional derivative term at the current moment be expressed as follows: in, , ; Spatial location coordinates; Represents the forward-modeling wave field across the entire space; The coefficient of the weighting term is defined as: ； This is the Gamma function.

7. The coal and rock stress inversion method as described in claim 1, characterized in that, In step S4, the steps for full waveform inversion include: S401. Using the inversion objective function Calculate the forward modeling response signal and seismic wave response signals Differences between them: in, To describe the propagation speed of P waves Spatial distribution function; Wave field evolution time The total length; S402. Calculate the objective function using the adjoint state method. right gradient ; S403. Iterative Update The optimized wave velocity field is obtained when the objective function converges.

8. The coal and rock stress inversion method as described in claim 1, characterized in that, In step S5, the specific operation of constructing the three-parameter coupled model of wave velocity, quality factor, and stress through posterior experimental calibration is as follows: Standard samples of various media are taken, and different static stresses are applied to each standard sample under laboratory conditions. At the same time, the wave velocity and quality factor of the samples are measured. The empirical curves of wave velocity and quality factor as a function of stress for each medium are fitted by multiple sets of wave velocity and quality factor data under static stress. The formula of the empirical curve is used as the formula of the three-parameter coupling model of wave velocity-quality factor-stress for that medium.

9. The coal and rock stress inversion method as described in claim 1, characterized in that, In step S6, high stress concentration regions and anomalous structures are identified in the three-dimensional stress field by combining threshold criteria and spatial gradient analysis. The specific operation is as follows: Extract the extreme regions in the stress field that exceed a preset multiple of the regional average stress and mark them as high stress concentration regions; Analyze the spatial gradients of the wave velocity field and stress field, locate regions where gradient abrupt changes coincide with stress distortion locations, and mark them as anomalous structures.

10. A coal and rock stress inversion system based on multi-material level set topology optimization, characterized in that, To implement the coal and rock stress inversion method as described in any one of claims 1-9, the method includes: The data acquisition module is used to set up artificial seismic sources and multiple detection points in the area to be detected within the coal and rock mass. The artificial seismic source generates a source signal and simultaneously acquires the seismic wave response signals of each detection point. The geological modeling and initial parameter field construction module is used to acquire empirical geological information of the area to be explored, and based on the empirical geological information, to model the geological structure using at least two level set functions, to obtain the level set function values ​​of each location point in the area to be explored, and to construct a three-dimensional spatial distribution map of materials based on the level set function values; and to map the level set function values ​​to wave velocity and quality factor, and to construct an initial multi-material physical property parameter field, including an initial three-dimensional wave velocity field and an initial three-dimensional quality factor field. The forward modeling module is used to perform forward modeling based on the source signal and the initial three-dimensional wave velocity field, using the time fractional viscoelastic wave equation, and outputs the forward modeling wave field in the whole space and the forward modeling response signal of each detection point. The full waveform inversion module is used to perform full waveform inversion based on seismic wave response signal, forward modeling response signal, full-space forward modeling wave field and initial wave velocity field to obtain optimized wave velocity field; The stress field construction module is used to construct a three-parameter coupled model of wave velocity, quality factor and stress by posterior experimental calibration. Based on the optimized wave velocity field and the initial quality factor field, the three-parameter coupled model of wave velocity, quality factor and stress is used to map the stress value at each location point to construct a three-dimensional stress field. The risk identification module is used to identify high stress concentration areas and anomalous structures in a three-dimensional stress field and output the risk identification results. The 3D visualization integration module is used to integrate the material's 3D spatial distribution map, the coordinates of artificial seismic sources and detection points, the optimized wave velocity field, the 3D stress field, and the risk identification results into the 3D visualization platform.

Images