Coal rock stress inversion method and system based on multi-material level set topological optimization

Abstract

The invention relates to the technical field of coal mine safety, in particular to a coal rock stress inversion method and system based on multi-material level set topological optimization, and the method comprises the steps: setting an artificial seismic source and a plurality of detection points, exciting a seismic source signal, and collecting a seismic wave response signal; utilizing at least two level set functions to construct a material three-dimensional space distribution map and mapping wave velocity and quality factors; performing forward modeling by using a time fractional order viscoelastic wave equation, and then performing full-waveform inversion to obtain an optimized wave velocity field; constructing a wave velocity-quality factor-stress three-parameter coupling model, and constructing a three-dimensional stress field; identifying a high stress concentration area and an abnormal structure, and outputting a risk identification result; and integrating the three-dimensional space distribution map of the material, the coordinates of the artificial seismic source and the detection point, the optimized wave velocity field, the three-dimensional stress field and the risk identification result into a three-dimensional visual platform. According to the invention, high-precision wave velocity reconstruction, stress quantitative conversion and risk dynamic display in the coal rock mass can be realized.

Description

Technical Field

[0001] The present application relates to the field of coal mine safety technology, and in particular to a coal rock stress inversion method and system based on multi-material level set topology optimization. Background Art

[0002] Coal seam rock burst is a sudden, highly destructive dynamic disaster deep within the coal rock mass, posing a serious threat to the safety of underground personnel and mining equipment. As coal seam mining depth increases, surrounding rock stress levels continue to rise, significantly increasing the probability of rock burst accidents. Therefore, timely characterization of the three-dimensional stress distribution within the coal 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 levels. The three-dimensional P-wave velocity field within coal and rock can indirectly indicate areas of high stress concentration, providing a basis for identifying rock burst risks. Therefore, existing technologies primarily capture rock failure events through microseismic monitoring and acoustic emission, combining seismic exploration techniques to invert the three-dimensional P-wave velocity field within coal and rock, indirectly predicting the internal stress distribution characteristics of coal and rock. Furthermore, the quality factor (Q) of the medium, as a measure of its stiffness and energy dissipation characteristics, is used in some existing technologies to infer changes in stress states.

[0004] However, existing methods for predicting the internal stress distribution characteristics of coal and rock masses still have significant limitations: first, the inversion process of the three-dimensional wave velocity field in coal and rock masses generally ignores the viscoelastic dissipation characteristics of coal and rock masses, resulting in waveform amplitude attenuation and distortion of the dispersion behavior description, affecting the accuracy of wave velocity field reconstruction; second, existing technologies mostly use grid units to directly invert physical parameters, which has a high calculation dimension and makes it difficult to clearly characterize the boundaries of various media, and the interface fuzziness problem is prominent; third, there is a lack of a coordinated constraint mechanism between wave velocity, quality factor and static stress, and a quantitative conversion relationship between dynamic wave field parameters and static stress field has not been established. Summary of the Invention

[0005] In response to the technical problems of existing coal and rock internal stress prediction methods, such as waveform distortion, interface fuzziness, and lack of multi-parameter constraints, which lead to insufficient accuracy of wave velocity field reconstruction, unclear medium boundary characterization, and inability to quantitatively convert dynamic parameters into static stress fields, this application provides a coal and rock stress inversion method and system based on multi-material level set topology optimization, which can realize high-precision wave velocity reconstruction, quantitative stress conversion and dynamic risk display inside the coal and rock mass.

[0006] In a first aspect, the present application provides a coal rock stress inversion method based on multi-material level set topology optimization, comprising the following steps: S1. Set up an artificial seismic source and multiple detection points in the area to be detected within the coal rock body, stimulate the source signal at the artificial seismic source, and synchronously collect the seismic wave response signals at each detection point. ; in, is the spatial position coordinate; is the wave field evolution time, starting from the excitation source signal; Number the test points; S2. Obtain empirical geological information of the area to be explored, and based on this empirical geological information, use at least two level set functions Carry out geological structure modeling to obtain each location point in the area to be detected The level set function value is used to construct the three-dimensional spatial distribution map of the material based on the level set function value; Map the level set function value to wave velocity and quality factor to construct the initial multi-material physical property parameter field, including the initial three-dimensional wave velocity field and the initial three-dimensional quality factor field; in, , is the total number of level set functions; The wave velocity is the propagation velocity of the P wave in the medium, and the three-dimensional wave velocity field is the spatial distribution of the propagation velocity of the P wave in the medium; S3. Based on the source signal and the initial 3D wave velocity field, forward modeling is performed using the time fractional viscoelastic wave equation, outputting the full-space forward modeled wave field and the forward modeled response signals of each detection point. In the forward simulation process, a short memory fast algorithm is used to solve the fractional derivative terms of the time fractional elastic-viscous wave equation. The full-space forward wave field is the acoustic pressure field or displacement field. S4. Based on seismic wave response signal , forward simulation response signal , perform full waveform inversion on the full-space forward wave field and the initial velocity field to obtain the optimized velocity field; S5. A three-parameter coupled model of velocity, quality factor, and stress is constructed using a posteriori experimental calibration method. Based on the optimized velocity field and initial quality factor field, the three-parameter coupled model is used to map the stress values ​​at each location and construct a three-dimensional stress field. S6. Identify high stress concentration areas and abnormal structures in the three-dimensional stress field and output risk identification results, including the coordinates of the high stress concentration areas and the boundary markers of the abnormal structures; S7. Integrate the 3D spatial distribution of materials, the coordinates of artificial seismic sources and detection points, the optimized wave velocity field, the 3D stress field, and risk identification results into a 3D visualization platform.

[0007] It should be further explained that in step S2, the empirical geological information includes the distribution of three media: the main body of coal rock, the water-bearing weak zone, and the hard interbedded gangue area; Level set function include 、 ,in: and Indicates location The medium at the location is mainly coal rock; and Indicates location The medium at the location is a water-bearing weak zone; Indicates location The medium at the location is a hard gangue area.

[0008] It should be further explained that, in step S2, the steps of using the multi-material level set function to perform geological structure modeling include: S201. Initialize the level set function ; S202. Iteratively updating the level set function using a topological evolution method based on the reaction-diffusion equation , the iterative update formula is:

[0009] in, represents the dynamic weight factor; represents the sensitivity normalization factor; represents the objective function value; represents the diffusion term that controls the smoothness of the boundary; represents the regularization coefficient; Indicates the location based on empirical geological information Department Known value.

[0010] It should be further explained that in step S202, the iterative update is terminated if at least one of the following conditions is met: The number of iterative updates reaches the preset maximum number of iterative updates; The difference between the objective function values ​​updated in two consecutive iterations is less than the preset convergence tolerance; In two consecutive iterative updates, the minimum value of the rate of change of the level set function at the same position is lower than the set threshold.

[0011] It should be further explained that, in step S2, the SIMP interpolation strategy is used to map the level set value to the wave velocity and quality factor.

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

[0013] in, is the spatial position coordinate; For location P-wave propagation velocity at ; represents the full-space forward wave field; is a scale parameter related to the relaxation characteristics of the medium; The Caputo type order is The time fractional derivative of is defined as:

[0014] in , is the Gamma function.

[0015] It should be further explained that the fractional derivative term of the time fractional-order elastic-viscous wave equation solved using the short-memory fast algorithm is specifically: Set the memory length to , the time step is , the current time is , assume the fractional derivative term at the current moment to be the following expression:

[0016] in, , ; is the weight term coefficient, defined as: .

[0017] It should be further explained that .

[0018] It should be further explained that in step S4, the full waveform inversion step includes: S401. Using the Inversion Objective Function Calculate the forward simulation response signal and seismic wave response signals The differences between:

[0019] in, To describe the P-wave propagation speed The spatial distribution function of TRON evolution time Total length; S402. Calculate the objective function using the adjoint state method right Gradient ; S403. Iterative Update Until the objective function converges, the optimized wave velocity field is obtained.

[0020] It should be further explained that step S403 uses an optimization algorithm to iteratively update , the optimization algorithms include conjugate gradient method or L-BFGS method.

[0021] It should be further explained that, in step S5, the specific operation of constructing the three-parameter coupling model of wave velocity, quality factor and stress by the a posteriori experimental calibration method is as follows: Standard samples of various media are taken, and different static stresses are applied to each standard sample under laboratory conditions. The wave velocity and quality factor of the sample are measured at the same time. The empirical curve of the wave velocity and quality factor of each medium changing with stress is 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 wave velocity-quality factor-stress three-parameter coupling model of the medium.

[0022] It should be further explained that in step S6, high stress concentration areas and abnormal structures are identified in the three-dimensional stress field by combining the threshold criterion and the spatial gradient analysis method. The specific operations are as follows: Extract the extreme value areas in the stress field that exceed the preset multiples of the regional average stress and mark them as high stress concentration areas; The spatial gradients of the wave velocity field and stress field are analyzed, and the areas where the gradient changes suddenly and coincide with the stress distortion position are located and marked as abnormal structures.

[0023] In a second aspect, the present application provides a coal rock stress inversion system based on multi-material level set topology optimization, which is used to implement the above-mentioned coal rock stress inversion method, including: The data acquisition module is used to set an artificial seismic source and multiple detection points in the area to be detected in the coal rock body, stimulate the source signal at the artificial seismic source, and synchronously collect the seismic wave response signal of each detection point; A geological modeling and initial parameter field construction module is used to obtain empirical geological information of the area to be explored, perform geological structure modeling based on the empirical geological information using at least two level set functions, obtain the level set function value for each location point in the area to be explored, and construct a three-dimensional spatial distribution map of the materials based on the level set function values; and map the level set function values ​​to wave velocity and quality factor 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 output the full-space forward modeling wave field and the forward modeling response signals of each detection point; Full waveform inversion module, used to perform full waveform inversion based on seismic wave response signals, forward simulation response signals, full-space forward wave fields and initial velocity fields to obtain optimized velocity fields; The stress field construction module is used to construct a three-parameter coupling model of wave velocity, quality factor, and stress through a posteriori experimental calibration method. Based on the optimized wave velocity field and initial quality factor field, the three-parameter coupling model of wave velocity, quality factor, and stress is used to map the stress value of each position point and construct a three-dimensional stress field. Risk identification module, used to identify high stress concentration areas and abnormal structures in the three-dimensional stress field and output risk identification results; The three-dimensional visualization integration module is used to integrate the three-dimensional spatial distribution map of materials, artificial seismic sources and detection point coordinates, optimized wave velocity field, three-dimensional stress field and risk identification results into the three-dimensional visualization platform.

[0024] In a third aspect, the present application provides an electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor is configured to implement the steps of the above-mentioned coal rock stress inversion method when executing the computer program.

[0025] In a fourth aspect, the present application provides a storage medium having a computer program stored thereon, which implements the steps of the above-mentioned coal rock stress inversion method when executed by a processor.

[0026] It can be seen from the above technical solutions that this application has the following advantages: 1. This application solves the problems of high computational dimension and fuzzy interface of existing grid inversion methods by utilizing at least two level set functions to model geological structures, construct a three-dimensional spatial distribution map of materials, map wave velocities and quality factors, 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.

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

[0028] 3. This application constructs a three-parameter coupling model of wave velocity, quality factor and stress through the method of a posteriori experimental calibration, maps the stress value of each position point and constructs a three-dimensional stress field, which solves the problem of lack of multi-parameter collaborative constraints and realizes the quantitative conversion from dynamic parameters to three-dimensional static stress field.

[0029] 4. This application identifies high-stress areas and abnormal structures based on the three-dimensional stress field and integrates the results into a visualization platform, solving the problem of insufficient early warning information output and achieving accurate labeling of risk sources and three-dimensional dynamic display. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] In order to more clearly illustrate the technical solution of the present application, the following is a brief introduction to the drawings required for the description. Obviously, the drawings described below are only some embodiments of the present application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without any creative work.

[0031] Figure 1 This is a flow chart of a coal rock stress inversion method based on multi-material level set topology optimization in one embodiment of the present application.

[0032] Figure 2 It is a schematic block diagram of a coal rock stress inversion system based on multi-material level set topology optimization in one embodiment of the present application.

[0033] Figure 3 It is a schematic diagram of the hardware structure of an electronic device in one embodiment of the present application. DETAILED DESCRIPTION

[0034] In order to make the application objectives, features, and advantages of this application more obvious and easy to understand, the technical solutions protected by this application will be clearly and completely described below using specific embodiments and drawings. Obviously, the embodiments described below are only part of the embodiments of this application, not all of them. Based on the embodiments in this application, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of this application.

[0035] The following describes in detail the coal rock stress inversion method involved in this application. Specific details, such as specific system structures and techniques, are provided for illustrative purposes, not for limitation, to facilitate a thorough understanding of the embodiments of this application. However, it should be apparent to those skilled in the art that this application may also be implemented in other embodiments without these specific details.

[0036] In the coal rock stress inversion method involved in this application, the term "comprising" is used to indicate the presence of the described features, integral bodies, steps, operations, elements and / or components, but does not exclude the presence or addition of one or more other features, integral bodies, steps, operations, elements, components and / or their collections. The terms "including," "comprising," "having" and their variations all mean "including but not limited to," unless otherwise specifically emphasized.

[0037] To facilitate the clear description of the technical solutions of this application, the words "first" and "second" are used to distinguish between identical or similar items with substantially the same functions and effects. Those skilled in the art will understand that the words "first" and "second" do not limit the quantity or order of execution, and the words "first" and "second" do not necessarily mean different.

[0038] The phrases "one embodiment" or "some embodiments" described in this application mean that the specific features, structures, or characteristics described in the embodiment are included in one or more embodiments of the application. Therefore, the phrases "in one embodiment," "in some embodiments," "in other embodiments," "in other embodiments," etc. that appear in different places in this application do not necessarily refer to the same embodiment, but rather mean "one or more but not all embodiments," unless otherwise specifically emphasized.

[0039] The technical solutions in the embodiments of the present application will be described clearly and completely below in conjunction with the drawings in the embodiments of the present application.

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

[0041] Figure 1 This is a flow chart of a coal rock stress inversion method based on multi-material level set topology optimization according to an embodiment of the present application. Figure 1 The execution subject may be a coal rock stress inversion system. According to different requirements, the order of the steps in the flow chart may be changed, and some steps may be omitted.

[0042] like Figure 1 As shown in Figure 2, the coal rock stress inversion method based on multi-material level set topology optimization includes: Step S1: Set up an artificial seismic source and multiple detection points in the area to be detected in the coal rock body, stimulate the source signal at the artificial seismic source, and synchronously collect the seismic wave response signal of each detection point. ; in, is the spatial position coordinate; is the wave field evolution time, starting from the excitation source signal; Number the test point.

[0043] By deploying artificial seismic sources and multiple detection points in the area to be detected and synchronously collecting seismic wave response signals, comprehensive, high signal-to-noise ratio observation data of the dynamic response of underground structures is obtained, providing a real and reliable input information source for subsequent inversion.

[0044] Step S2: obtain empirical geological information of the area to be detected, and use at least two level set functions based on the empirical geological information. Carry out geological structure modeling to obtain each location point in the area to be detected The level set function value is used to construct the three-dimensional spatial distribution map of the material based on the level set function value; Map the level set function value to wave velocity and quality factor to construct the initial multi-material physical property parameter field, including the initial three-dimensional wave velocity field and the initial three-dimensional quality factor field; in, , is the total number of level set functions; The wave velocity is the P-wave propagation velocity in the medium, and the three-dimensional wave velocity field is the spatial distribution of the P-wave propagation velocity in the medium.

[0045] By using multi-material level set functions to model geological structures based on empirical geological information and mapping the level set values ​​into wave velocities and quality factors, the reasonable construction of the initial three-dimensional physical parameter field (wave velocity field, quality factor field) of complex coal and rock formations (containing multiple media) is achieved, providing an initial physical model that conforms to the geological background for the forward modeling of the wave equation.

[0046] In some specific embodiments, the empirical geological information includes the distribution of three media: coal rock body, water-bearing weak zone and hard interbedded gangue area; Level set function include 、 ,in: and Indicates location The medium at the location is mainly coal rock; and Indicates location The medium at the location is a water-bearing weak zone; Indicates location The medium at the location is a hard gangue area.

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

[0048] In some specific embodiments, the step of using a multi-material level set function to perform geological structure modeling includes: S201. Initialize the level set function ; S202. Iteratively updating the level set function using a topological evolution method based on the reaction-diffusion equation , the iterative update formula is:

[0049] in, represents the dynamic weight factor; represents the sensitivity normalization factor; represents the objective function value; represents the diffusion term that controls the smoothness of the boundary; represents the regularization coefficient; Indicates the location based on empirical geological information Department Known value.

[0050] By iteratively updating the level set function using a topological evolution method based on the reaction-diffusion equation, 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 a precise foundation for the subsequent initialization of the physical property parameter field.

[0051] In some specific embodiments, in step S202, the iterative update is terminated if at least one of the following conditions is met: The number of iterative updates reaches the preset maximum number of iterative updates; The difference between the objective function values ​​updated in two consecutive iterations is less than the preset convergence tolerance; In two consecutive iterative updates, the minimum value of the rate of change of the level set function at the same position is lower than the set threshold.

[0052] By setting clear iterative termination conditions, including the maximum number of iterations, the tolerance for objective function changes, and the threshold for the level set function change rate, intelligent control of the level set topology evolution process and optimization of computing resources are achieved. This can effectively prevent the algorithm from falling into invalid iterations or premature convergence, ensuring a balance between efficiency and accuracy in geological modeling.

[0053] In some embodiments, a SIMP interpolation strategy is used to map level set values ​​to wave speeds and quality factors.

[0054] By using an interpolation strategy to map the level set function values ​​into physical parameters, including wave velocity and quality factor, a smooth transition from a discrete geometric model to a continuous physical parameter field is achieved. This can avoid non-physical jumps in physical parameters at the material interface, ensure the physical consistency and numerical stability of the initial wave velocity field and quality factor field, and support subsequent forward and inversion processes.

[0055] Step S3: Based on the source signal and the initial three-dimensional wave velocity field, the time fractional viscoelastic wave equation is used to perform forward simulation, and the full-space forward wave field and the forward simulation response signal of each detection point are output. In the forward simulation process, a short memory fast algorithm is used to solve the fractional derivative terms of the time fractional elastic-viscous wave equation, and the full-space forward wave field is the acoustic pressure field or displacement field.

[0056] By using the time fractional-order viscoelastic wave equation for forward simulation based on the initial wave velocity field and using the short-memory fast algorithm to solve the fractional-order derivative terms, an accurate numerical simulation of the seismic wave field propagation process with attenuation and dispersion effects is achieved, and the theoretical simulation response signal and spatial wave field evolution at each detection point are obtained, providing forward prediction data for inversion.

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

[0058] in, is the spatial position coordinate; For location P-wave propagation velocity at ; represents the full-space forward wave field; is a scale parameter related to the relaxation characteristics of the medium; The Caputo type order is The time fractional derivative of is defined as:

[0059] in , is the Gamma function.

[0060] By using the viscoelastic wave equation containing time fractional derivative terms to describe wave propagation, an accurate mathematical representation of the dispersion and attenuation effects of seismic waves in coal and rock formations is achieved. This can significantly improve the fit of the forward simulation to the actual observed wave field and provide a high-fidelity forward modeling engine for full waveform inversion.

[0061] In some specific embodiments, the fractional-order derivative term of the time fractional-order elastic-viscous wave equation is solved using a short-memory fast algorithm as follows: Set the memory length to , the time step is , the current time is , assume the fractional derivative term at the current moment to be the following expression:

[0062] in, , ; is the weight term coefficient, defined as: .

[0063] By adopting a short-memory fast algorithm to approximate the calculation of fractional-order derivative terms (using only historical wave field values ​​of a finite time step and specific weight coefficients), the computational complexity of the solution process of the fractional-order viscoelastic wave equation is significantly reduced. This can greatly improve the computational efficiency of forward simulation while ensuring accuracy, making full waveform inversion feasible on an engineering scale.

[0064] In some specific embodiments, .

[0065] By limiting the upper limit of the memory length of the short-memory fast algorithm, further optimization control of calculation accuracy and efficiency is achieved, which can effectively suppress the cumulative error caused by excessively long memory and control memory consumption, thereby ensuring the numerical stability and efficiency of the forward simulation.

[0066] Step S4: Based on the seismic wave response signal , forward simulation response signal , perform full waveform inversion on the full-space forward wave field and the initial velocity field to obtain the optimized velocity field.

[0067] By performing full waveform inversion based on the measured response signal, forward simulation response signal, spatial wave field and initial wave velocity field, the iterative optimization and correction of the initial wave velocity model is achieved, and a high-precision three-dimensional optimized wave velocity field that is more consistent with the actual observation data is obtained.

[0068] In some specific embodiments, the step of full waveform inversion includes: S401. Using the Inversion Objective Function Calculate the forward simulation response signal and seismic wave response signals The differences between:

[0069] in, To describe the P-wave propagation speed The spatial distribution function of TRON evolution time Total length; S402. Calculate the objective function using the adjoint state method right Gradient ; S403. Iterative Update Until the objective function converges, the optimized wave velocity field is obtained.

[0070] By defining an inversion objective function based on the sum of squares of waveform residuals and applying a gradient optimization algorithm for iterative updating, the optimal estimation of the spatial distribution function of wave velocity is achieved in the least squares sense. The difference between the observed wave field and the simulated wave field can be minimized, thereby obtaining a high-resolution optimized wave velocity field that conforms to actual geophysical laws.

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

[0072] By limiting the use of a specific gradient optimization algorithm to update the velocity field, an efficient and stable decrease of the objective function is achieved during the iteration process, which can accelerate the convergence speed of the full waveform inversion and enhance the robustness of the algorithm under complex models, ensuring the acquisition of a reliable optimized velocity field.

[0073] Step S5: construct a three-parameter coupling model of wave velocity, quality factor, and stress through a posteriori experimental calibration method. Based on the optimized wave velocity field and the initial quality factor field, the three-parameter coupling model of wave velocity, quality factor, and stress is used to map the stress value of each position point to construct a three-dimensional stress field.

[0074] By utilizing the velocity-quality factor-stress coupling model calibrated by a posteriori experiments and combining the optimized velocity field and initial quality factor field to map stress values, the dynamic physical parameters obtained by seismic wave inversion are converted into a static three-dimensional absolute stress field, providing direct mechanical state information for risk assessment.

[0075] In some specific embodiments, the specific operations for constructing the three-parameter coupling model of wave velocity, quality factor, and stress by the a posteriori experimental calibration method are as follows: Standard samples of various media are taken, and different static stresses are applied to each standard sample under laboratory conditions. The wave velocity and quality factor of the sample are measured at the same time. The empirical curve of the wave velocity and quality factor of each medium changing with stress is 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 wave velocity-quality factor-stress three-parameter coupling model of the medium.

[0076] By adopting the laboratory calibration method, an empirical coupling relationship model between the wave velocity, quality factor and stress of different media under static stress is established. The wave velocity field and quality factor information obtained by field inversion are used to quantitatively map the absolute stress value inside the coal rock mass. The seismic wave inversion results can be converted into a three-dimensional stress field that can be directly used for engineering risk assessment.

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

[0078] By using quantitative methods to identify high stress concentration areas and abnormal structures in the three-dimensional stress field, the automatic positioning and boundary characterization of potential dangerous sources of coal and rock dynamic disasters (such as rock burst) are achieved, and specific risk identification results are output to guide safety prevention and control.

[0079] In some specific embodiments, a method combining threshold criteria and spatial gradient analysis is used to identify high stress concentration areas and abnormal structures in a three-dimensional stress field. The specific operations are as follows: Extract the extreme value areas in the stress field that exceed the preset multiples of the regional average stress and mark them as high stress concentration areas; The spatial gradients of the wave velocity field and stress field are analyzed, and the areas where the gradient changes suddenly and coincide with the stress distortion position are located and marked as abnormal structures.

[0080] By combining stress threshold criteria with physical field gradient analysis to identify risk areas, the automated, quantitative positioning and boundary extraction of high stress concentration areas (corresponding to abnormally high stress values) and potential abnormal structures (corresponding to stress / wave velocity gradient mutation areas) are achieved. This can accurately identify potential sources of danger for disasters such as rock burst and guide the formulation of prevention and control measures.

[0081] Step S7: Integrate 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 the three-dimensional visualization platform.

[0082] By integrating the geological model, observation system, velocity field, stress field and risk identification results obtained by inversion into a three-dimensional visualization platform, a comprehensive and intuitive display of the stress state of coal and rock and the spatial distribution of risks is achieved, providing an intuitive and efficient interactive analysis environment for geological interpretation, engineering decision-making and disaster warning.

[0083] In some specific embodiments, the three-dimensional visualization platform is based on the three-dimensional spatial coordinates inside the coal rock, and displays the coordinates of the artificial seismic source and detection point, the optimized wave velocity field, the three-dimensional stress field, the risk identification results, the mine tunnel layout, the drilling measurement points, etc. in a three-dimensional overlay display. The platform software supports multi-window display, real-time data input, background calculation and result query and other functions.

[0084] For example, on the main interface, the three-dimensional structural model of the coal seam, the location and frequency of real-time microseismic events, the readings of each monitoring borehole stress sensor, and the comprehensive risk assessment index can be displayed simultaneously. Users can rotate and zoom the three-dimensional view to intuitively view the relative position and range of the high-stress area in the tunnel.

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

[0086] The 3D visualization platform can provide timely forecasts and early warnings for potential rock bursts caused by high stress concentrations in actual production, significantly improving the scientific nature and effectiveness of mine safety management. With the 3D visualization platform, mine workers can monitor coal seam stress evolution trends and potential rock burst hazards around the clock, shifting from passive rescue to active early warning.

[0087] In a specific embodiment, the steps of the coal rock stress inversion method based on multi-material level set topology optimization include: Step S1: Set up an artificial seismic source and multiple detection points in the area to be detected in the coal rock body, stimulate the source signal at the artificial seismic source, and synchronously collect the seismic wave response signal of each detection point. ; in, is the spatial position coordinate; is the wave field evolution time, starting from the excitation source signal; Number the test point.

[0088] Step S2: obtain empirical geological information of the area to be detected, and use at least two level set functions based on the empirical geological information. Carry out geological structure modeling to obtain each location point in the area to be detected The level set function value is used to construct the three-dimensional spatial distribution map of the material based on the level set function value; Among them, the empirical geological information includes the distribution of three kinds of media: coal rock body, water-bearing weak zone and hard interbedded gangue area; Level set function include 、 ,in: and Indicates location The medium at the location is mainly coal rock; and Indicates location The medium at the location is a water-bearing weak zone; Indicates location The medium at the location is a hard gangue area; The steps for geological structure modeling using multi-material level set functions include: S201. Initialize the level set function ; S202. Iteratively updating the level set function using a topological evolution method based on the reaction-diffusion equation , the iterative update formula is:

[0089] in, represents the dynamic weight factor; represents the sensitivity normalization factor; represents the objective function value; represents the diffusion term that controls the smoothness of the boundary; represents the regularization coefficient; Indicates the location based on empirical geological information Department known value; The iterative update is terminated if at least one of the following conditions is met: The number of iterative updates reaches the preset maximum number of iterative updates; The difference between the objective function values ​​updated in two consecutive iterations is less than the preset convergence tolerance; In two consecutive iterative updates, the minimum value of the rate of change of the level set function at the same position is lower than the set threshold; The SIMP interpolation strategy is used to map the level set function values ​​into wave velocity and quality factor, and the initial multi-material physical property parameter field is constructed, including the initial three-dimensional wave velocity field and the initial three-dimensional quality factor field; The wave velocity is the P-wave propagation velocity in the medium, and the three-dimensional wave velocity field is the spatial distribution of the P-wave propagation velocity in the medium.

[0090] Step S3: Based on the source signal and the initial three-dimensional wave velocity field, the time fractional viscoelastic wave equation is used to perform forward simulation, and the full-space forward wave field and the forward simulation response signal of each detection point are output. In the forward simulation process, a short memory fast algorithm is used to solve the fractional derivative terms of the time fractional elastic-viscous wave equation. The full-space forward wave field is the acoustic pressure field or displacement field. The time fractional viscoelastic wave equation is expressed as:

[0091] in, is the spatial position coordinate; For location P-wave propagation velocity at ; represents the full-space forward wave field; is a scale parameter related to the relaxation characteristics of the medium; The Caputo type order is The time fractional derivative of is defined as:

[0092] in , is the Gamma function; The fractional derivative terms of the time fractional-order elastic-viscous wave equation are solved using the short-memory fast algorithm as follows: Set the memory length to , the time step is , the current time is , assume the fractional derivative term at the current moment to be the following expression:

[0093] in, , ; is the weight term coefficient, defined as: , .

[0094] Step S4: Based on the seismic wave response signal , forward simulation response signal , perform full waveform inversion on the full-space forward wave field and the initial velocity field to obtain the optimized velocity field; The steps of full waveform inversion include: S401. Using the Inversion Objective Function Calculate the forward simulation response signal and seismic wave response signals The differences between:

[0095] in, To describe the P-wave propagation speed The spatial distribution function of TRON evolution time Total length; S402. Calculate the objective function using the adjoint state method right Gradient ; S403. Iterative update using optimization algorithm When the objective function converges, the optimized wave velocity field is obtained. The optimization algorithm includes the conjugate gradient method or the L-BFGS method.

[0096] Step S5, constructing a three-parameter coupling model of wave velocity, quality factor, and stress through a method of a posteriori experimental calibration. Based on the optimized wave velocity field and the initial quality factor field, the three-parameter coupling model of wave velocity, quality factor, and stress is used to map the stress value of each position point to construct a three-dimensional stress field. The specific operation of constructing the three-parameter coupling model of wave velocity, quality factor and stress through the a posteriori experimental calibration method is as follows: Standard samples of various media are taken, and different static stresses are applied to each standard sample under laboratory conditions. The wave velocity and quality factor of the sample are measured at the same time. The empirical curve of the wave velocity and quality factor of each medium changing with stress is 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 wave velocity-quality factor-stress three-parameter coupling model of the medium.

[0097] Step S6, identifying high stress concentration areas and abnormal structures in the three-dimensional stress field, specifically the following operations: Extract the extreme value areas in the stress field that exceed the preset multiples of the regional average stress and mark them as high stress concentration areas; Analyze the spatial gradients of the wave velocity field and stress field, locate the areas where the gradient changes suddenly and coincide with the stress distortion position, and mark them as abnormal structures; Output risk identification results, including coordinates of high stress concentration areas and abnormal structural boundary markers.

[0098] Step S7: Integrate 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 the three-dimensional visualization platform.

[0099] The following is an embodiment of a coal rock stress inversion system based on multi-material level set topology optimization provided in an embodiment of the present application. The coal rock stress inversion system based on multi-material level set topology optimization and the coal rock stress inversion methods of the above-mentioned embodiments belong to the same inventive concept. For details not fully described in the embodiment of the coal rock stress inversion system, please refer to the embodiment of the coal rock stress inversion method based on multi-material level set topology optimization.

[0100] like Figure 2 As shown in the figure, the coal rock stress inversion system based on multi-material level set topology optimization includes: The data acquisition module is used to set an artificial seismic source and multiple detection points in the area to be detected in the coal rock body, stimulate the source signal at the artificial seismic source, and synchronously collect the seismic wave response signal of each detection point; A geological modeling and initial parameter field construction module is used to obtain empirical geological information of the area to be explored, perform geological structure modeling based on the empirical geological information using at least two level set functions, obtain the level set function value for each location point in the area to be explored, and construct a three-dimensional spatial distribution map of the materials based on the level set function values; and map the level set function values ​​to wave velocity and quality factor 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 output the full-space forward modeling wave field and the forward modeling response signals of each detection point; Full waveform inversion module, used to perform full waveform inversion based on seismic wave response signals, forward simulation response signals, full-space forward wave fields and initial velocity fields to obtain optimized velocity fields; The stress field construction module is used to construct a three-parameter coupling model of wave velocity, quality factor, and stress through a posteriori experimental calibration method. Based on the optimized wave velocity field and initial quality factor field, the three-parameter coupling model of wave velocity, quality factor, and stress is used to map the stress value of each position point and construct a three-dimensional stress field. Risk identification module, used to identify high stress concentration areas and abnormal structures in the three-dimensional stress field and output risk identification results; The three-dimensional visualization integration module is used to integrate the three-dimensional spatial distribution map of materials, artificial seismic sources and detection point coordinates, optimized wave velocity field, three-dimensional stress field and risk identification results into the three-dimensional visualization platform.

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

[0102] This application also provides an electronic device for implementing each embodiment of this application. Figure 3 A hardware structure diagram of an electronic device for implementing various embodiments of the present application is shown in FIG. Figure 3 As shown, the electronic device includes a memory, a processor, and a computer program stored in the memory and executable on the processor.

[0103] Those skilled in the art will understand that the electronic device structure involved in the embodiments of the present 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 arrange the components differently.

[0104] In the embodiments of the present 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 the present application described and / or claimed herein.

[0105] In the embodiment of the present application, the processor can be implemented by 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, and an electronic unit designed to perform the functions described herein. In some cases, such an embodiment can be implemented in a controller. For software implementation, an embodiment such as a process or function can be implemented with a separate software module that allows the execution of at least one function or operation. The software code can be implemented by a software application (or program) written in any appropriate programming language, and the software code can be stored in a memory and executed by a controller.

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

[0107] Those skilled in the art will appreciate that various aspects of the electronic device provided herein may be implemented as a system, method, or program product. Therefore, various aspects of the present application may be specifically implemented in the following forms: a complete hardware implementation, a complete software implementation (including firmware, microcode, etc.), or a combination of hardware and software implementations, which may be collectively referred to herein as "circuits," "modules," or "systems."

[0108] This application also provides a storage medium storing a program product capable of implementing a coal rock stress inversion method based on multi-material level set topology optimization. In some possible implementations, various aspects of this application may also be implemented in the form of a program product comprising program code. When the program product is executed on a terminal device, the program code is configured to cause the terminal device to execute the steps described in the "Exemplary Methods" section above according to various exemplary embodiments of this application.

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

[0110] The above description of the disclosed embodiments is intended to enable one skilled in the art to implement or use the present application. Various modifications to these embodiments will be readily apparent to one skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present application. Therefore, the present application is not limited to the embodiments shown herein, but is intended to conform to the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A coal rock stress inversion method based on multi-material level set topology optimization, characterized in that: include: S1. Set up an artificial seismic source and multiple detection points in the area to be detected within the coal rock body, stimulate the source signal at the artificial seismic source, and synchronously collect the seismic wave response signals at each detection point. ; in, is the spatial position coordinate; is the wave field evolution time, starting from the excitation source signal; Number the test points; S2. Obtain empirical geological information of the area to be explored, and based on this empirical geological information, use at least two level set functions Carry out geological structure modeling to obtain each location point in the area to be detected The level set function value is used to construct the three-dimensional spatial distribution map of the material based on the level set function value; Map the level set function value to wave velocity and quality factor to construct the initial multi-material physical property parameter field, including the initial three-dimensional wave velocity field and the initial three-dimensional quality factor field; in, , is the total number of level set functions; The wave velocity is the propagation velocity of the P wave in the medium, and the three-dimensional wave velocity field is the spatial distribution of the propagation velocity of the P wave in the medium; S3. Based on the source signal and the initial 3D wave velocity field, forward modeling is performed using the time fractional viscoelastic wave equation, outputting the full-space forward modeled wave field and the forward modeled response signals of each detection point. In the forward simulation process, a short memory fast algorithm is used to solve the fractional derivative terms of the time fractional elastic-viscous wave equation. The full-space forward wave field is the acoustic pressure field or displacement field. S4. Based on seismic wave response signal , forward simulation response signal , perform full waveform inversion on the full-space forward wave field and the initial velocity field to obtain the optimized velocity field; S5. A three-parameter coupled model of velocity, quality factor, and stress is constructed using a posteriori experimental calibration method. Based on the optimized velocity field and initial quality factor field, the three-parameter coupled model is used to map the stress values ​​at each location and construct a three-dimensional stress field. S6. Identify high stress concentration areas and abnormal structures in the three-dimensional stress field and output risk identification results, including the coordinates of the high stress concentration areas and the boundary markers of the abnormal structures; S7. Integrate the 3D spatial distribution of materials, the coordinates of artificial seismic sources and detection points, the optimized wave velocity field, the 3D stress field, and risk identification results into a 3D visualization platform.

2. The coal rock stress inversion method according to claim 1, characterized in that: In step S2, the steps of using the multi-material level set function to perform geological structure modeling include: S201. Initialize the level set function ; S202. Iteratively updating the level set function using a topological evolution method based on the reaction-diffusion equation , the iterative update formula is: in, represents the dynamic weight factor; represents the sensitivity normalization factor; represents the objective function value; represents the diffusion term that controls the smoothness of the boundary; represents the regularization coefficient; Indicates the location based on empirical geological information Department Known value.

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

4. The coal rock stress inversion method according to claim 1, characterized in that: In step S2, the SIMP interpolation strategy is used to map the level set value to wave velocity and quality factor.

5. The coal rock stress inversion method according to claim 1, characterized in that: In step S3, the time fractional-order viscoelastic wave equation is expressed as: in, is the spatial position coordinate; For location P-wave propagation velocity at ; represents the full-space forward wave field; is a scale parameter related to the relaxation characteristics of the medium; The Caputo type order is The time fractional derivative of is defined as: in , is the Gamma function.

6. The coal rock stress inversion method according to claim 1, characterized in that: The fractional derivative terms of the time fractional-order elastic-viscous wave equation are solved using the short-memory fast algorithm as follows: Set the memory length to , the time step is , the current time is , assume the fractional derivative term at the current moment to be the following expression: in, , ; is the weight term coefficient, defined as: 。 7. The coal rock stress inversion method according to claim 1, characterized in that: In step S4, the full waveform inversion step includes: S401. Using the Inversion Objective Function Calculate the forward simulation response signal and seismic wave response signals The differences between: in, To describe the P-wave propagation speed The spatial distribution function of TRON evolution time Total length; S402. Calculate the objective function using the adjoint state method right Gradient ; S403. Iterative Update Until the objective function converges, the optimized wave velocity field is obtained.

8. The coal rock stress inversion method according to claim 1, characterized in that: In step S5, the specific operation of constructing the three-parameter coupling model of wave velocity, quality factor and stress by the a posteriori experimental calibration method is as follows: Standard samples of various media are taken, and different static stresses are applied to each standard sample under laboratory conditions. The wave velocity and quality factor of the sample are measured at the same time. The empirical curve of the wave velocity and quality factor of each medium changing with stress is 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 wave velocity-quality factor-stress three-parameter coupling model of the medium.

9. The coal rock stress inversion method according to claim 1, characterized in that: In step S6, high stress concentration areas and abnormal structures are identified in the three-dimensional stress field by combining threshold criterion and spatial gradient analysis. The specific operations are as follows: Extract the extreme value areas in the stress field that exceed the preset multiples of the regional average stress and mark them as high stress concentration areas; The spatial gradients of the wave velocity field and stress field are analyzed, and the areas where the gradient changes suddenly and coincide with the stress distortion position are located and marked as abnormal structures.

10. A coal rock stress inversion system based on multi-material level set topology optimization, characterized in that: A method for implementing the coal rock stress inversion method according to any one of claims 1 to 9, comprising: The data acquisition module is used to set an artificial seismic source and multiple detection points in the area to be detected in the coal rock body, stimulate the source signal at the artificial seismic source, and synchronously collect the seismic wave response signal of each detection point; A geological modeling and initial parameter field construction module is used to obtain empirical geological information of the area to be explored, perform geological structure modeling based on the empirical geological information using at least two level set functions, obtain the level set function value for each location point in the area to be explored, and construct a three-dimensional spatial distribution map of the materials based on the level set function values; and map the level set function values ​​to wave velocity and quality factor 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 output the full-space forward modeling wave field and the forward modeling response signals of each detection point; Full waveform inversion module, used to perform full waveform inversion based on seismic wave response signals, forward simulation response signals, full-space forward wave fields and initial velocity fields to obtain optimized velocity fields; The stress field construction module is used to construct a three-parameter coupling model of wave velocity, quality factor, and stress through a posteriori experimental calibration method. Based on the optimized wave velocity field and initial quality factor field, the three-parameter coupling model of wave velocity, quality factor, and stress is used to map the stress value of each position point and construct a three-dimensional stress field. Risk identification module, used to identify high stress concentration areas and abnormal structures in the three-dimensional stress field and output risk identification results; The three-dimensional visualization integration module is used to integrate the three-dimensional spatial distribution map of materials, artificial seismic sources and detection point coordinates, optimized wave velocity field, three-dimensional stress field and risk identification results into the three-dimensional visualization platform.

Images