Goaf detection method based on transient electromagnetic and seismic wave field joint inversion

Abstract

According to the goaf detection method based on transient electromagnetic and seismic wave field joint inversion, data acquisition and joint inversion are carried out on the goaf through the combination of the transient electromagnetic detection system and the reflection seismic detection system, the complementary advantages of the two detection methods are played, the bottlenecks of a single method in the aspects of resolution and signal-to-noise ratio are broken through, and the goaf detection method based on transient electromagnetic and seismic wave field joint inversion is realized. Fine characterization of the structure, the physical property and the stability of the goaf is realized; a machine learning classification algorithm is introduced to analyze a goaf multi-attribute feature data set, and intelligent discrimination of the type, scale and stability of the goaf under complex geological conditions is realized; establishing a quantitative response relationship between goaf physical property parameters and goaf electrical and elastic difference characteristic parameters, constructing a goaf transient electromagnetic and seismic wave field joint forward modeling numerical model, and constructing an electromagnetic-seismic joint inversion objective function; and calculating and solving to obtain goaf resistivity, seismic wave velocity three-dimensional distribution and uncertainty quantitative evaluation results, and reducing blind areas of goaf judgment.

Description

Technical Field

[0001] The invention belongs to the technical field of geological engineering and relates to a method for detecting goaf areas. Background Art

[0002] The goaf caused by coal mining refers to the hollow area formed after coal seam mining, which mainly includes three parts: coal seam goaf, tunnel engineering and chamber. The formation and evolution of goaf is directly related to the surface subsidence, karst collapse, water resource destruction and ecological environment deterioration in the mining area, and is the bottleneck restricting the sustainable development of the mining area. Therefore, timely detection and identification of the spatial distribution and stability of goaf is of great significance to ensure the safety of coal mining, guide the ecological restoration of mining areas, and promote the economic transformation of mining areas.

[0003] Conventional methods for goaf detection mainly include drilling, geophysical exploration and mine surveying, but the above methods have the problems of high cost, low resolution and poor real-time performance. In recent years, geophysical exploration technology has been widely used and rapidly developed in goaf detection. For example, transient electromagnetic method was used to conduct detection research on a coal mining subsidence area, revealing the three-dimensional electrical structure of the subsidence area; shallow seismic reflection method was used to determine the burial depth and range of the goaf roof damage zone; comprehensive gravity and apparent resistivity data were used to invert the density and resistivity distribution of the goaf. However, the current joint inversion methods for goaf detection are mostly limited to qualitative descriptions, lacking a unified mathematical model and inversion framework, and it is difficult to accurately characterize the fine structure of the goaf under complex conditions. Summary of the invention

[0004] In order to solve the problem that the joint inversion method for goaf detection described in the background technology is difficult to accurately depict the fine structure of the goaf under complex conditions, the present invention provides a goaf detection method based on the joint inversion of transient electromagnetic and seismic wave fields.

[0005] The method of the present invention comprises the following steps:

[0006] According to the burial conditions of the goaf, a transient electromagnetic detection system and a reflection seismic detection system are deployed in the tunnel crossing area where the goaf is located;

[0007] Collect transient electromagnetic apparent resistivity and reflection seismic waveform data, and perform preprocessing including data registration, noise suppression, and synchronous imaging to extract electrical and elastic difference characteristic parameters of the goaf area;

[0008] Based on the electrical and elastic difference characteristic parameters of the goaf, a multi-attribute feature data set is extracted. The multi-attribute feature data set is trained, verified and tested using a machine learning classification algorithm to establish a goaf classification model to identify the type, scale and stability of the goaf.

[0009] Based on Maxwell's equations and elastic wave equation, finite difference and finite element methods are used for discrete solution, and the quantitative response relationship between the physical parameters of the goaf and the electrical and elastic difference characteristic parameters of the goaf is established, and a joint forward numerical model of transient electromagnetic and seismic wave fields in the goaf is constructed.

[0010] Based on the forward numerical model of transient electromagnetic and seismic wave fields in goaf, combined with the prior information of goaf including geometric shape, surrounding rock stress state and geology, the electromagnetic-seismic joint inversion objective function is constructed.

[0011] According to the electromagnetic-seismic joint inversion objective function, the Bayesian-Markov Chain Monte Carlo method is used to perform inversion and solve the problem, and the three-dimensional distribution of resistivity and seismic wave velocity in the goaf and the quantitative evaluation results of their uncertainty are obtained.

[0012] According to the identification results of goaf type, scale and stability, combined with the three-dimensional distribution of goaf resistivity and seismic wave velocity and quantitative assessment results of their uncertainty, a three-dimensional fine structure and physical property model of goaf is established.

[0013] Furthermore, the theoretical basis of the transient electromagnetic is Maxwell's equations, which are expressed as follows:

[0014]

[0015] Where, E is the electric field intensity, H is the magnetic field intensity, Js is the current source density, D is the electric displacement, B is the magnetic induction intensity, σ is the conductivity, and ρ is the charge density;

[0016] Considering the quasi-static approximation conditions and ignoring the displacement current, the above equations can be simplified to:

[0017]

[0018] In the numerical simulation of transient electromagnetic in goaf, the staggered grid finite difference method is used to discretely solve the above equations. First, the study area is divided into Nx×Ny×Nz cubic grids, and then the electric field and magnetic field components are respectively configured at the edges and center points of the grids, and the parameters such as conductivity and permeability are assigned at the center of the grid; the central difference format with second-order accuracy is used in the time direction, and the staggered grid difference operator with second-order accuracy is used in the space direction. The three-dimensional transient electromagnetic forward numerical equation is established as follows:

[0019] En+1-En-1,

[0020]

[0021] Hn+1-Hn,

[0022]

[0023] The above equations are assembled on the entire space grid and time step to form a large sparse linear equation system, which can be solved by direct or iterative methods to obtain the three-dimensional electromagnetic field space-time distribution.

[0024] Furthermore, the numerical simulation of reflection earthquake adopts elastic wave equation, which is expressed as:

[0025]

[0026] Where u is the displacement vector, σ is the stress tensor, ρ is the density, and f is the source term;

[0027] The stress-strain relationship satisfies Hooke's law:

[0028] σ=λtr(ε)δ+2με,

[0029]

[0030] Where δ is the unit tensor, λ and μ are the first and second Lame constants respectively;

[0031] In the numerical solution, the staggered grid finite difference method is used, the stress component is configured at the center of the grid, and the velocity component is configured at the midpoint of the corresponding edge; the central difference format of second-order accuracy is used in both time and space directions to establish the elastic wave equation difference format in the velocity-stress form:

[0032]

[0033] Where, v is the velocity component;

[0034] A giant sparse matrix equation set is assembled for all grids and time steps, and is solved using an explicit time-stepping method to obtain the spatiotemporal evolution information of the three-dimensional elastic wave field.

[0035] Furthermore, the electrical and elastic difference characteristic parameters of the goaf include transient electromagnetic multi-attribute characteristic parameters and reflection seismic multi-attribute characteristic parameters; the transient electromagnetic multi-attribute characteristic parameters include apparent resistivity, apparent conductivity, induced polarization coefficient, and transient electromagnetic waveform attenuation coefficient; the reflection seismic multi-attribute characteristic parameters include seismic wave first arrival time, amplitude, frequency, phase, reflection coefficient, and attenuation coefficient.

[0036] Furthermore, the method for establishing the goaf classification model includes:

[0037] S1. Based on the electrical and elastic difference characteristic parameters of the goaf, a multi-attribute characteristic data set is extracted, wherein the electrical and elastic difference characteristic parameters of the goaf include transient electromagnetic multi-attribute characteristic parameters and reflection seismic multi-attribute characteristic parameters, the transient electromagnetic multi-attribute characteristic parameters include apparent resistivity, apparent conductivity, induced polarization coefficient, transient electromagnetic waveform attenuation coefficient, and the reflection seismic multi-attribute characteristic parameters include seismic wave first arrival travel time, amplitude, frequency, phase, reflection coefficient, and attenuation coefficient;

[0038] S2, using Pearson correlation coefficient and mutual information index to analyze the correlation between each feature in the multi-attribute feature data set and the goaf type, and using recursive feature elimination and sequential feature selection methods to select the best feature combination as the feature subset with the greatest contribution to goaf classification;

[0039] S3, using random stratified sampling to divide the feature subset into three parts: training set, validation set and test set;

[0040] S4. The machine learning classification algorithm uses a support vector machine to construct an optimal classification hyperplane in a high-dimensional feature space to maximize the sample intervals of different categories. The training set is input into the support vector machine for training. The mathematical model of the support vector machine is:

[0041] min12||w||2+C∑i=1Nξ i ,

[0042] styi(w·x i +b)≥1-ξ i , ξi≥0, i=1, 2,...,N;

[0043] Among them, is the classification hyperplane normal vector, b is the bias term, C is the penalty factor, ξ i is the slack variable;

[0044] By solving the above optimization problem, we get the classification decision function:

[0045] f(x)=sign(∑i=1Nα i y i (x·xi)+b),

[0046] In the formula, α i is the Lagrange multiplier, y i is the support vector, (x·xi) is the kernel function;

[0047] During the training process, grid search and cross-validation methods are used to optimize model hyperparameters including kernel function type, penalty factor, and kernel function parameters to obtain the trained model;

[0048] S5. Use the test set to evaluate the performance of the goaf classification model of the trained model, and use the accuracy, recall rate, F1 value, and ROC curve indicators to comprehensively evaluate the classification effect; at the same time, the learning curve analyzes the trend of the accuracy of the training set and the validation set as the sample size changes to determine whether the trained model has overfitting or underfitting; for the trained model with poor performance, return to steps S2-S4, adjust the feature selection, and optimize the model parameters until a model with satisfactory classification performance is obtained, which is used as the goaf classification model.

[0049] Furthermore, the method for constructing the goaf area transient electromagnetic and seismic wave field joint forward numerical model includes:

[0050] S1. Based on prior geological information, the modeling software is used to construct the initial model of the goaf, including the spatial distribution of the goaf and the resistivity, density, and wave velocity physical parameters of the surrounding rock;

[0051] S2. Setting the excitation parameters of the transient electromagnetic emission source including the emission current and waveform and the receiving system parameters including the number of receiving channels and sampling interval, and configuring the spatial positions of the emission source and the receiving point in the initial model of the goaf;

[0052] S3. Based on the spatial and temporal distribution of the three-dimensional electromagnetic field and the spatial positions of the transmitting source and the receiving point, the three-dimensional transient electromagnetic response of the initial model of the goaf is forward calculated to obtain the electric field and magnetic field components at different time and spatial positions, and synthesize the apparent resistivity characterization parameters;

[0053] S4, setting reflection seismic excitation source parameters including source type and main frequency and receiving system parameters including number of detector groups and track spacing, and configuring the spatial positions of shot points and detector points in the initial model;

[0054] S5. Based on the spatial and temporal evolution information of the three-dimensional elastic wave field and the spatial positions of the shot points and receiver points, forward model the three-dimensional seismic wave field snapshot of the initial model of the goaf area, extract the seismic records of different gathers, and pick up attribute parameters including first arrival travel time, frequency dispersion, and amplitude attenuation;

[0055] S6. By changing the geometric parameters of the initial model of the goaf, including the goaf range and thickness, and the physical parameters, including the resistivity and density of the goaf filling material, steps S2-S6 are repeated to carry out parameter sensitivity analysis of electromagnetic and seismic responses, and quantitatively characterize the quantitative response relationship between the physical parameters of the goaf and the electrical and elastic difference characteristic parameters of the goaf;

[0056] S7. Based on the electric field and magnetic field components at different time and space positions, the apparent resistivity characterization parameters, the seismic records of different gathers, the attribute parameters including the first arrival travel time, dispersion, and amplitude attenuation, and the quantitative response relationship between the physical properties of the goaf and the electrical and elastic difference characteristic parameters of the goaf, a joint forward numerical model of transient electromagnetic and seismic wave fields in the goaf is constructed.

[0057] Furthermore, the electromagnetic-seismic joint inversion objective function is:

[0058] Φ(m)=1||Wd(df(m))|| 2 +α||Wm(m-m0)|| 2 +β||L(m)|| 2 ,

[0059] Wherein, the first term is the data fitting term, d is the observed data, including apparent resistivity, seismic travel time, and amplitude information; f(m) is the theoretical response of the goaf model m, and Wd is the data weighting matrix, which is used to balance the contribution of different data types; the second term is the model constraint term, m0 is the reference model, which is usually set based on geological data and experience, Wm is the model weighting matrix, which constrains the spatial distribution of the inversion model, and α is the model constraint factor, which controls the influence of prior information on the inversion results; the third term is the model smoothness term, L is the spatial second-order difference operator, and β is the smoothness constraint factor, which controls the smoothness of the inversion model, where the model parameter m mainly includes the resistivity, density, and P- and S-wave velocities of each grid in the goaf.

[0060] Furthermore, the solution process of the electromagnetic-seismic joint inversion objective function is:

[0061] The electromagnetic-seismic joint inversion objective function is transformed into solving the following minimization problem:

[0062] minΦ(m),stml≤m≤mu,

[0063] Where ml and mu are the upper and lower limits of the model parameters, respectively, reflecting the known geological constraints;

[0064] The Bayesian-Markov Chain Monte Carlo method is used for inversion. Given the observed data d, the posterior probability p(m|d) of the goaf model m is proportional to the product of the likelihood function p(d|m) and the prior probability p(m), that is: p(m|d)∝p(d|m)p(m); where the likelihood function p(d|m) reflects the degree of fit between the model m and the observed data d, and the prior probability p(m) reflects the physical constraints and statistical characteristics of the model parameters.

[0065] By constructing a Markov chain to realize random sampling of the posterior probability distribution, the optimal estimate and confidence interval of the goaf model parameters are obtained. The process is as follows:

[0066] S1. Based on geological data and empirical judgment, an initial model of the goaf is constructed as the starting point of the Markov chain;

[0067] S2. Based on the Metropolis-Hastings criterion, the current model mi is perturbed to the next state mi+1 with probability p; the expression of the transition probability p is:

[0068] p=min{1, p(mi+1|d)q(mi|mi+1)p(mi|d)q(mi+1|mi)},

[0069] Where q(mi+1|mi) is the transfer kernel function of the Markov chain, and a symmetric distribution such as Gaussian distribution is usually selected;

[0070] S3, repeat step S2, continuously update and adjust the state of the Markov chain until the convergence condition is reached: the posterior probability does not change significantly;

[0071] S4. Based on the stationary distribution of Markov chain, a series of inversion model samples are randomly selected and statistically analyzed to obtain the optimal estimation and confidence interval of the model parameters of the goaf area, and the three-dimensional distribution of resistivity and seismic wave velocity in the goaf area and the quantitative evaluation results of their uncertainty.

[0072] The present invention uses a transient electromagnetic detection system and a reflection seismic detection system to jointly collect data and jointly invert the goaf, giving full play to the complementary advantages of the two detection methods, breaking through the bottlenecks of a single method in terms of resolution and signal-to-noise ratio, and realizing the fine characterization of the structure, physical properties, and stability of the goaf; a machine learning classification algorithm is introduced to analyze the multi-attribute feature data set of the goaf, realizing the intelligent discrimination of the type, scale, and stability of the goaf under complex geological conditions, and expanding the universality of the method; a quantitative response relationship between the physical property parameters of the goaf and the electrical and elastic difference characteristic parameters of the goaf is established, a forward numerical model of the transient electromagnetic and seismic wave field of the goaf is constructed, an electromagnetic-seismic joint inversion objective function is constructed, and a Bayesian-Markov chain Monte Carlo method is used for inversion solution, and the three-dimensional distribution of the resistivity and seismic wave velocity of the goaf and the quantitative evaluation results of their uncertainty are obtained by calculation, thereby reducing the blind spots in the judgment of the goaf.

[0073] Compared with the prior art, the present invention has the following beneficial effects: First, through the joint inversion of transient electromagnetic and seismic wave fields, the problem of low resolution and insufficient reliability of a single geophysical exploration method in goaf detection is effectively solved, and meter-level fine characterization of goaf is achieved; Second, the efficiency and accuracy of intelligent identification and classification of goaf are significantly improved through machine learning classification algorithms; Third, the inversion solution process effectively reduces the uncertainty of goaf inversion, and the inversion results are in good agreement with actual verification. The three-dimensional morphological characterization of the three-dimensional fine structure and physical property model of the goaf that is finally established is more realistic and reliable, and can accurately characterize the fine structure of the goaf under complex geological conditions; Fourth, it has versatility and scalability, and has broad application prospects in the fields of coal mine goaf investigation, karst collapse area detection, and underground engineering safety evaluation, providing key technical support for the safe construction of major projects. BRIEF DESCRIPTION OF THE DRAWINGS

[0074] Figure 1 The figure is a flow chart of the method of the present invention.

[0075] Figure 2 It is a mapping relationship diagram of the goaf development degree and physical property parameters of the embodiment.

[0076] Figure 3 It is a three-dimensional fine structure model diagram of the goaf area of ​​the embodiment.

[0077] Figure 4 This is a tunnel optimization line plan diagram for an embodiment. DETAILED DESCRIPTION

[0078] In order to make the technical problems, technical solutions and beneficial effects to be solved by the present application more clearly understood, the present application is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application and are not used to limit the present application.

[0079] The goaf detection method based on the joint inversion of transient electromagnetic and seismic wave fields is shown in the flowchart. Figure 1 As shown, the details are as follows.

[0080] Firstly, according to the burial conditions of the goaf, a transient electromagnetic detection system and a reflection seismic detection system are deployed in the tunnel crossing area where the goaf is located.

[0081] Specifically, the detection frequency of the transient electromagnetic detection system is 0.1-100Hz, and the transmission power is not less than 500W; the excitation frequency of the reflection seismic detection system is 20-200Hz, and a controllable source is used for excitation. According to the burial conditions of the goaf, the spatial layout and survey line parameters of the two detection systems are optimized.

[0082] Then, transient electromagnetic apparent resistivity and reflection seismic waveform data are collected, and preprocessing including data registration, noise suppression, and synchronous imaging is performed to extract the electrical and elastic difference characteristic parameters of the goaf area.

[0083] Specifically, the electrical and elastic difference characteristic parameters of the goaf are composed of transient electromagnetic multi-attribute characteristic parameters and reflection seismic multi-attribute characteristic parameters. The transient electromagnetic multi-attribute characteristic parameters are: apparent resistivity, apparent conductivity, induced polarization coefficient, transient electromagnetic waveform attenuation coefficient; the reflection seismic multi-attribute characteristic parameters are: seismic wave first arrival time, amplitude, frequency, phase, reflection coefficient, attenuation coefficient.

[0084] Then, based on the electrical and elastic difference characteristic parameters of the goaf, a multi-attribute feature data set was extracted. The machine learning classification algorithm was used to train, verify and test the multi-attribute feature data set, and a goaf classification model was established to identify the type, scale and stability of the goaf.

[0085] Specifically, the method for establishing the goaf classification model includes:

[0086] S1. Based on the electrical and elastic difference characteristic parameters of the goaf, a multi-attribute feature data set is extracted;

[0087] S2, using Pearson correlation coefficient and mutual information index to analyze the correlation between each feature in the multi-attribute feature data set and the goaf type, and using recursive feature elimination and sequential feature selection methods to select the best feature combination as the feature subset with the greatest contribution to goaf classification;

[0088] S3, using random stratified sampling to divide the feature subset into three parts: training set, validation set and test set;

[0089] S4. The machine learning classification algorithm uses a support vector machine to construct an optimal classification hyperplane in a high-dimensional feature space to maximize the sample intervals of different categories. The training set is input into the support vector machine for training. The mathematical model of the support vector machine is:

[0090] min12||w||2+C∑i=1Nξ i ,

[0091] styi(w·x i +b)≥1-ξ i , ξi≥0, i=1, 2,...,N;

[0092] Among them, is the classification hyperplane normal vector, b is the bias term, C is the penalty factor, ξ i is the slack variable;

[0093] By solving the above optimization problem, we get the classification decision function:

[0094] f(x)=sign(∑i=1Nα i y i (x·xi)+b),

[0095] In the formula, α i is the Lagrange multiplier, y i is the support vector, (x·xi) is the kernel function;

[0096] During the training process, grid search and cross-validation methods are used to optimize model hyperparameters including kernel function type, penalty factor, and kernel function parameters to obtain the trained model;

[0097] S5. Use the test set to evaluate the performance of the goaf classification model of the trained model, and use the accuracy, recall rate, F1 value, and ROC curve indicators to comprehensively evaluate the classification effect; at the same time, the learning curve analyzes the trend of the accuracy of the training set and the validation set as the sample size changes to determine whether the trained model has overfitting or underfitting; for the trained model with poor performance, return to steps S2-S4, adjust the feature selection, and optimize the model parameters until a model with satisfactory classification performance is obtained, which is used as the goaf classification model.

[0098] Then, based on Maxwell's equations and elastic wave equation, finite difference and finite element methods were used for discrete solution, and the quantitative response relationship between the physical properties of the goaf and the electrical and elastic difference characteristic parameters of the goaf was established, and a joint forward numerical model of transient electromagnetic and seismic wave fields in the goaf was constructed.

[0099] Specifically, the theoretical basis of transient electromagnetics is Maxwell's equations, which are expressed as:

[0100]

[0101] Where, E is the electric field intensity, H is the magnetic field intensity, Js is the current source density, D is the electric displacement, B is the magnetic induction intensity, σ is the conductivity, and ρ is the charge density;

[0102] Considering the quasi-static approximation conditions and ignoring the displacement current, the above equations can be simplified to:

[0103]

[0104] In the numerical simulation of transient electromagnetic in goaf, the staggered grid finite difference method is used to discretely solve the above equations. First, the study area is divided into Nx×Ny×Nz cubic grids, and then the electric field and magnetic field components are respectively configured at the edges and center points of the grids, and the parameters such as conductivity and permeability are assigned at the center of the grid; the central difference format with second-order accuracy is used in the time direction, and the staggered grid difference operator with second-order accuracy is used in the space direction. The three-dimensional transient electromagnetic forward numerical equation is established as follows:

[0105] En+1-En-1,

[0106]

[0107] Hn+1-Hn,

[0108]

[0109] The above equations are assembled on the entire space grid and time step to form a large sparse linear equation system, which can be solved by direct or iterative methods to obtain the three-dimensional electromagnetic field space-time distribution.

[0110] Specifically, the numerical simulation of reflection seismic adopts the elastic wave equation, which is expressed as:

[0111]

[0112] Where u is the displacement vector, σ is the stress tensor, ρ is the density, and f is the source term;

[0113] The stress-strain relationship satisfies Hooke's law:

[0114] σ=λtr(ε)δ+2με,

[0115]

[0116] Where δ is the unit tensor, λ and μ are the first and second Lame constants respectively;

[0117] In the numerical solution, the staggered grid finite difference method is used, the stress component is configured at the center of the grid, and the velocity component is configured at the midpoint of the corresponding edge; the central difference format of second-order accuracy is used in both time and space directions to establish the elastic wave equation difference format in the velocity-stress form:

[0118]

[0119] Where, v is the velocity component;

[0120] A giant sparse matrix equation set is assembled for all grids and time steps, and is solved using an explicit time-stepping method to obtain the spatiotemporal evolution information of the three-dimensional elastic wave field.

[0121] More specifically, the method for constructing the goaf area transient electromagnetic and seismic wave field joint forward numerical model includes:

[0122] S1. Based on prior geological information, the modeling software is used to construct the initial model of the goaf, including the spatial distribution of the goaf and the resistivity, density, and wave velocity physical parameters of the surrounding rock;

[0123] S2. Setting the excitation parameters of the transient electromagnetic emission source including the emission current and waveform and the receiving system parameters including the number of receiving channels and sampling interval, and configuring the spatial positions of the emission source and the receiving point in the initial model of the goaf;

[0124] S3. Based on the spatial and temporal distribution of the three-dimensional electromagnetic field and the spatial positions of the transmitting source and the receiving point, the three-dimensional transient electromagnetic response of the initial model of the goaf is forward calculated to obtain the electric field and magnetic field components at different time and spatial positions, and synthesize the apparent resistivity characterization parameters;

[0125] S4, setting reflection seismic excitation source parameters including source type and main frequency and receiving system parameters including number of detector groups and track spacing, and configuring the spatial positions of shot points and detector points in the initial model;

[0126] S5. Based on the spatial and temporal evolution information of the three-dimensional elastic wave field and the spatial positions of the shot points and receiver points, forward model the three-dimensional seismic wave field snapshot of the initial model of the goaf area, extract the seismic records of different gathers, and pick up attribute parameters including first arrival travel time, frequency dispersion, and amplitude attenuation;

[0127] S6. By changing the geometric parameters of the initial model of the goaf, including the goaf range and thickness, and the physical parameters, including the resistivity and density of the goaf filling material, steps S2-S6 are repeated to carry out parameter sensitivity analysis of electromagnetic and seismic responses, and quantitatively characterize the quantitative response relationship between the physical parameters of the goaf and the electrical and elastic difference characteristic parameters of the goaf;

[0128] S7. Based on the electric field and magnetic field components at different time and space positions, the apparent resistivity characterization parameters, the seismic records of different gathers, the attribute parameters including the first arrival travel time, dispersion, and amplitude attenuation, and the quantitative response relationship between the physical properties of the goaf and the electrical and elastic difference characteristic parameters of the goaf, a joint forward numerical model of transient electromagnetic and seismic wave fields in the goaf is constructed.

[0129] Then, based on the forward numerical model of transient electromagnetic and seismic wave fields in the goaf, combined with the prior information of the goaf including the geometric shape and surrounding rock stress state geology, the electromagnetic-seismic joint inversion objective function is constructed.

[0130] Specifically, the objective function of electromagnetic-seismic joint inversion is:

[0131] Φ(m)=1||Wd(df(m))|| 2 +α||Wm(m-m0)|| 2 +β||L(m)|| 2 ,

[0132] Wherein, the first term is the data fitting term, d is the observed data, including apparent resistivity, seismic travel time, and amplitude information; f(m) is the theoretical response of the goaf model m, and Wd is the data weighting matrix, which is used to balance the contribution of different data types; the second term is the model constraint term, m0 is the reference model, which is usually set based on geological data and experience, Wm is the model weighting matrix, which constrains the spatial distribution of the inversion model, and α is the model constraint factor, which controls the influence of prior information on the inversion results; the third term is the model smoothness term, L is the spatial second-order difference operator, and β is the smoothness constraint factor, which controls the smoothness of the inversion model, where the model parameter m mainly includes the resistivity, density, and P- and S-wave velocities of each grid in the goaf.

[0133] Then, according to the electromagnetic-seismic joint inversion objective function, the Bayesian-Markov Chain Monte Carlo method was used to perform inversion and solve the three-dimensional distribution of resistivity and seismic wave velocity in the goaf and the quantitative evaluation results of their uncertainty were obtained.

[0134] Specifically, the solution process of the electromagnetic-seismic joint inversion objective function is:

[0135] The electromagnetic-seismic joint inversion objective function is transformed into solving the following minimization problem:

[0136] minΦ(m),stml≤m≤mu,

[0137] Where ml and mu are the upper and lower limits of the model parameters, respectively, reflecting the known geological constraints;

[0138] The Bayesian-Markov Chain Monte Carlo method is used for inversion. Given the observed data d, the posterior probability p(m|d) of the goaf model m is proportional to the product of the likelihood function p(d|m) and the prior probability p(m), that is: p(m|d)∝p(d|m)p(m); where the likelihood function p(d|m) reflects the degree of fit between the model m and the observed data d, and the prior probability p(m) reflects the physical constraints and statistical characteristics of the model parameters.

[0139] By constructing a Markov chain to realize random sampling of the posterior probability distribution, the optimal estimate and confidence interval of the goaf model parameters are obtained. The process is as follows:

[0140] S1. Based on geological data and empirical judgment, an initial model of the goaf is constructed as the starting point of the Markov chain;

[0141] S2. Based on the Metropolis-Hastings criterion, the current model mi is perturbed to the next state mi+1 with probability p; the expression of the transition probability p is:

[0142] p=min{1, p(mi+1|d)q(mi|mi+1)p(mi|d)q(mi+1|mi)},

[0143] Where q(mi+1|mi) is the transfer kernel function of the Markov chain, and a symmetric distribution such as Gaussian distribution is usually selected;

[0144] S3, repeat step S2, continuously update and adjust the state of the Markov chain until the convergence condition is reached: the posterior probability does not change significantly;

[0145] S4. Based on the stationary distribution of Markov chain, a series of inversion model samples are randomly selected and statistically analyzed to obtain the optimal estimate and confidence interval of the goaf model parameters, and the three-dimensional distribution of resistivity and seismic wave velocity in the goaf and the quantitative evaluation results of their uncertainty.

[0146] Finally, based on the identification results of goaf type, scale and stability, combined with the three-dimensional distribution of goaf resistivity and seismic wave velocity and quantitative assessment results of their uncertainty, a three-dimensional fine structure and physical property model of the goaf was established.

[0147] Example

[0148] In view of a water conservancy project in which a water delivery tunnel passes through an old goaf section, the present invention is used to detect the goaf section.

[0149] According to the known distribution of goaf, 8 transient electromagnetic survey lines were designed, with a line length of 5km and a measurement point spacing of 20m; 4 seismic survey lines were laid out, with a line length of 5km and a detection point spacing of 5m. The transient electromagnetic emission source uses the EMP-400T transmitter of the American EMIT company, with a maximum transmission power of 800A·m 2 , frequency 0.3Hz~300Hz; the seismic excitation source uses the BIS-SHD50 electric spark source of IDE Company of the United States, and the single-point excitation energy is 2000J. The transient electromagnetic receiving equipment uses the V8 multi-channel receiver of Phoenix Company of Canada, with a total of 150 receiving channels; the seismic acquisition uses the ARIES-Ⅱ system of Inova Company, with a total of 1000 detection points.

[0150] After the acquisition of transient electromagnetic apparent resistivity and seismic three-component waveform data of the entire survey line was completed, the electromagnetic data was preprocessed by instrument response correction, synchronous stacking, and three-dimensional spatial interpolation, and the signal-to-noise ratio was increased by more than 5 times; static correction, amplitude recovery, and surface wave attenuation preprocessing were performed on the seismic data to effectively suppress high-frequency random noise, optimize the event axis response, and obtain high-quality post-stack migration profiles. 12 attributes such as seismic reflection wave dispersion, amplitude, and apparent velocity were extracted to collaboratively characterize the electrical and elastic differences in the goaf.

[0151] A typical goaf section was selected to construct a seismic-electromagnetic wave field coupling model of three layers of media, namely, cover layer, bedrock and goaf. According to geological data, the filling types of goaf were set: fissure type, collapse type and water-filled type, with a thickness of 1m-20m and a burial depth of 10m-100m. Full response numerical simulation was carried out to obtain the response characteristics of apparent resistivity of 5Ω·m-500Ω·m and seismic wave velocity of 300m / s-2500m / s for different types of goaf, and a mapping relationship of "goaf development degree-physical property parameters" was established, such as Figure 2 As shown, constraints are provided for the joint inversion.

[0152] Using forward simulation data as training samples, the goaf classification model was constructed by integrating 8-dimensional multi-source attributes such as measured apparent resistivity and seismic wave dispersion. Under the constraints of well data, intelligent calibration of goaf type, scale, and stability classification labels was performed, and field collected data was classified and identified. The accuracy of goaf identification reached 85%, which is more than 20% higher than the traditional empirical interpretation method.

[0153] Taking transient electromagnetic apparent resistivity and seismic wave reflection velocity as constraints, and adding prior geological information such as the burial law of the goaf and the stress state of the surrounding rock, the electromagnetic-seismic joint inversion objective function is constructed:

[0154] min{||W d (d obs -F(m))|| 2 +α||W m (m-m0)|| 2 +β||Lm|| 2},

[0155] Among them, m is the spatial distribution and physical property parameters of the goaf; dobs is the measured geophysical data; F is the forward operator; Wd and Wm are the data and model weighting matrices respectively; m0 is the prior model; L is the smooth constraint operator; α and β are the prior information and smooth constraint weights. The Markov chain Monte Carlo method is used to randomly perturb the goaf model, obtain the posterior distribution of the inversion results, generate the confidence probability map, and quantitatively evaluate the uncertainty of the inversion parameters.

[0156] Based on the results of transient electromagnetic and seismic joint inversion and uncertainty analysis, combined with drilling exposure, a three-dimensional fine structure model of the goaf along the XX water diversion tunnel was established. Figure 3 As shown in the figure, 7 high-risk goaf sections were identified, with an overall positioning accuracy of better than 5m and a scale estimation error of less than 10%. The maximum burial depth of the goaf is more than 50 meters, and the volume is more than 50,000 cubic meters. Based on the stability prediction of the goaf, a tunnel optimization route plan was proposed, such as Figure 4As shown in the figure, two strongly coupled goaf areas were partially avoided, and the risk factor of goaf instability was reduced by more than 80%, saving 30 million yuan in project investment. The site was guided to carry out targeted goaf area grouting reinforcement treatment, achieving safe and efficient crossing construction of the XX water diversion tunnel, creating significant engineering and social benefits.

[0157] Those skilled in the art will appreciate that the embodiments of the present application can be provided as methods, systems, or computer program products. Therefore, the application can adopt the form of complete hardware embodiment, complete software embodiment, or the embodiment in combination with software and hardware. Moreover, the application can adopt the form of the computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) that contain computer-usable program code. The scheme in the embodiment of the present application can be implemented in various computer languages, for example, object-oriented programming language Java, C++, Python and literal scripting language JavaScript, etc.

[0158] The present application is described with reference to the flowcharts and / or block diagrams of the methods, devices (systems), and computer program products according to the embodiments of the present application. It should be understood that each process and / or box in the flowchart and / or block diagram, as well as the combination of the processes and / or boxes in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to generate a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 A process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0159] These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing device to operate in a specific manner, so that the instructions stored in the computer-readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 A process or multiple processes and / or boxes Figure 1 A function specified in one or more boxes.

[0160] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing instructions for implementing the process. Figure 1 A process or multiple processes and / or boxes Figure 1 The steps for the functions specified in one or more boxes.

[0161] Although the preferred embodiments of the present application have been described, those skilled in the art may make other changes and modifications to these embodiments once they have learned the basic creative concept. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments and all changes and modifications falling within the scope of the present application.

[0162] Obviously, those skilled in the art can make various changes and modifications to the present application without departing from the spirit and scope of the present application. Thus, if these modifications and variations of the present application fall within the scope of the claims of the present application and their equivalents, the present application is also intended to include these modifications and variations.

Claims

1. A method for detecting goaf areas based on the joint inversion of transient electromagnetic and seismic wave fields, characterized in that: The following steps are involved: According to the burial conditions of the goaf, a transient electromagnetic detection system and a reflection seismic detection system are deployed in the tunnel crossing area where the goaf is located; Collect transient electromagnetic apparent resistivity and reflection seismic waveform data, and perform preprocessing including data registration, noise suppression, and synchronous imaging to extract electrical and elastic difference characteristic parameters of the goaf area; Based on the electrical and elastic difference characteristic parameters of the goaf, a multi-attribute feature data set is extracted. The multi-attribute feature data set is trained, verified and tested using a machine learning classification algorithm to establish a goaf classification model to identify the type, scale and stability of the goaf. Based on Maxwell's equations and elastic wave equation, finite difference and finite element methods are used for discrete solution, and the quantitative response relationship between the physical parameters of the goaf and the electrical and elastic difference characteristic parameters of the goaf is established, and a joint forward numerical model of transient electromagnetic and seismic wave fields in the goaf is constructed. Based on the forward numerical model of transient electromagnetic and seismic wave fields in goaf, combined with the prior information of goaf including geometric shape, surrounding rock stress state and geology, the electromagnetic-seismic joint inversion objective function is constructed. According to the electromagnetic-seismic joint inversion objective function, the Bayesian-Markov Chain Monte Carlo method is used to perform inversion and solve the problem, and the three-dimensional distribution of resistivity and seismic wave velocity in the goaf and the quantitative evaluation results of their uncertainty are obtained. According to the identification results of goaf type, scale and stability, combined with the three-dimensional distribution of goaf resistivity and seismic wave velocity and quantitative assessment results of their uncertainty, a three-dimensional fine structure and physical property model of goaf is established.

2. The method for goaf detection based on transient electromagnetic and seismic wave field joint inversion according to claim 1 is characterized in that: The theoretical basis of the transient electromagnetic is Maxwell's equations, which are expressed as follows: Where, E is the electric field intensity, H is the magnetic field intensity, Js is the current source density, D is the electric displacement, B is the magnetic induction intensity, σ is the conductivity, and ρ is the charge density; Considering the quasi-static approximation conditions and ignoring the displacement current, the above equations can be simplified to: In the numerical simulation of transient electromagnetic in goaf, the staggered grid finite difference method is used to discretely solve the above equations. First, the study area is divided into Nx×Ny×Nz cubic grids, and then the electric field and magnetic field components are respectively configured at the edges and center points of the grids, and the parameters such as conductivity and permeability are assigned at the center of the grid; the central difference format with second-order accuracy is used in the time direction, and the staggered grid difference operator with second-order accuracy is used in the space direction. The three-dimensional transient electromagnetic forward numerical equation is established as follows: En+1-En-1, Hn+1-Hn, The above equations are assembled on the entire space grid and time step to form a large sparse linear equation system, which can be solved by direct or iterative methods to obtain the three-dimensional electromagnetic field space-time distribution.

3. The method for goaf detection based on transient electromagnetic and seismic wave field joint inversion according to claim 2 is characterized in that: The numerical simulation of reflection earthquake adopts elastic wave equation, which is expressed as: Where u is the displacement vector, σ is the stress tensor, ρ is the density, and f is the source term; The stress-strain relationship satisfies Hooke's law: σ=λtr(ε)δ+2με, Where δ is the unit tensor, λ and μ are the first and second Lame constants respectively; In the numerical solution, the staggered grid finite difference method is used, the stress component is configured at the center of the grid, and the velocity component is configured at the midpoint of the corresponding edge; the central difference format of second-order accuracy is used in both time and space directions to establish the elastic wave equation difference format in the velocity-stress form: Where, v is the velocity component; A giant sparse matrix equation set is assembled for all grids and time steps, and is solved using an explicit time-stepping method to obtain the spatiotemporal evolution information of the three-dimensional elastic wave field.

4. The method for goaf detection based on transient electromagnetic and seismic wave field joint inversion according to claim 3 is characterized in that: The electrical and elastic difference characteristic parameters of the goaf include transient electromagnetic multi-attribute characteristic parameters and reflected seismic multi-attribute characteristic parameters; the transient electromagnetic multi-attribute characteristic parameters include apparent resistivity, apparent conductivity, induced polarization coefficient, and transient electromagnetic waveform attenuation coefficient; the reflected seismic multi-attribute characteristic parameters include seismic wave first arrival time, amplitude, frequency, phase, reflection coefficient, and attenuation coefficient.

5. The method for goaf detection based on transient electromagnetic and seismic wave field joint inversion according to claim 4 is characterized in that: The method for establishing the goaf classification model comprises: S1. Based on the electrical and elastic difference characteristic parameters of the goaf, a multi-attribute feature data set is extracted; S2, using Pearson correlation coefficient and mutual information index to analyze the correlation between each feature in the multi-attribute feature data set and the goaf type, and using recursive feature elimination and sequential feature selection methods to select the best feature combination as the feature subset with the greatest contribution to goaf classification; S3, using random stratified sampling to divide the feature subset into three parts: training set, validation set and test set; S4. The machine learning classification algorithm uses a support vector machine to construct an optimal classification hyperplane in a high-dimensional feature space to maximize the sample intervals of different categories. The training set is input into the support vector machine for training. The mathematical model of the support vector machine is: min12||w||2+C∑i=1Nξ i , s.t.yi(w·x i +b)≥1-ξ i ,ξi≥0,i=1,2,…,N; Among them, is the classification hyperplane normal vector, b is the bias term, C is the penalty factor, ξ i is the slack variable; By solving the above optimization problem, we get the classification decision function: f(x)=sign(∑i=1Nα i y i (x·xi)+b), In the formula, α i is the Lagrange multiplier, y i is the support vector, (x·xi) is the kernel function; During the training process, grid search and cross-validation methods are used to optimize model hyperparameters including kernel function type, penalty factor, and kernel function parameters to obtain the trained model; S5. Use the test set to evaluate the performance of the goaf classification model of the trained model, and use the accuracy, recall rate, F1 value, and ROC curve indicators to comprehensively evaluate the classification effect; at the same time, the learning curve analyzes the trend of the accuracy of the training set and the validation set as the sample size changes to determine whether the trained model has overfitting or underfitting; for the trained model with poor performance, return to steps S2-S4, adjust the feature selection, and optimize the model parameters until a model with satisfactory classification performance is obtained, which is used as the goaf classification model.

6. The method for goaf detection based on transient electromagnetic and seismic wave field joint inversion according to claim 5 is characterized in that: The method for constructing the goaf area transient electromagnetic and seismic wave field joint forward numerical model includes: S1. Based on prior geological information, the modeling software is used to construct the initial model of the goaf, including the spatial distribution of the goaf and the resistivity, density, and wave velocity physical parameters of the surrounding rock; S2. Setting the excitation parameters of the transient electromagnetic emission source including the emission current and waveform and the receiving system parameters including the number of receiving channels and sampling interval, and configuring the spatial positions of the emission source and the receiving point in the initial model of the goaf; S3. Based on the spatial and temporal distribution of the three-dimensional electromagnetic field and the spatial positions of the transmitting source and the receiving point, the three-dimensional transient electromagnetic response of the initial model of the goaf is forward calculated to obtain the electric field and magnetic field components at different time and spatial positions, and synthesize the apparent resistivity characterization parameters; S4, setting reflection seismic excitation source parameters including source type and main frequency and receiving system parameters including number of detector groups and track spacing, and configuring the spatial positions of shot points and detector points in the initial model; S5. Based on the spatial and temporal evolution information of the three-dimensional elastic wave field and the spatial positions of the shot points and receiver points, forward model the three-dimensional seismic wave field snapshot of the initial model of the goaf area, extract the seismic records of different gathers, and pick up attribute parameters including first arrival travel time, frequency dispersion, and amplitude attenuation; S6. By changing the geometric parameters of the initial model of the goaf, including the goaf range and thickness, and the physical parameters, including the resistivity and density of the goaf filling material, steps S2-S6 are repeated to carry out parameter sensitivity analysis of electromagnetic and seismic responses, and quantitatively characterize the quantitative response relationship between the physical parameters of the goaf and the electrical and elastic difference characteristic parameters of the goaf; S7. Based on the electric field and magnetic field components at different time and space positions, the apparent resistivity characterization parameters, the seismic records of different gathers, the attribute parameters including the first arrival travel time, dispersion, and amplitude attenuation, and the quantitative response relationship between the physical properties of the goaf and the electrical and elastic difference characteristic parameters of the goaf, a joint forward numerical model of transient electromagnetic and seismic wave fields in the goaf is constructed.

7. The method for goaf detection based on transient electromagnetic and seismic wave field joint inversion according to claim 6 is characterized in that: The electromagnetic-seismic joint inversion objective function is: Φ(m)=1||Wd(d-f(m))|| 2 +α||Wm(m-m0)|| 2 +β||L(m)|| 2 , Wherein, the first term is the data fitting term, d is the observed data, including apparent resistivity, seismic travel time, and amplitude information; f(m) is the theoretical response of the goaf model m, and Wd is the data weighting matrix, which is used to balance the contribution of different data types; the second term is the model constraint term, m0 is the reference model, which is usually set based on geological data and experience, Wm is the model weighting matrix, which constrains the spatial distribution of the inversion model, and α is the model constraint factor, which controls the influence of prior information on the inversion results; the third term is the model smoothness term, L is the spatial second-order difference operator, and β is the smoothness constraint factor, which controls the smoothness of the inversion model, where the model parameter m mainly includes the resistivity, density, and P- and S-wave velocities of each grid in the goaf.

8. The method for goaf detection based on transient electromagnetic and seismic wave field joint inversion according to claim 7 is characterized in that: The solution process of the electromagnetic-seismic joint inversion objective function is: The electromagnetic-seismic joint inversion objective function is transformed into solving the following minimization problem: minΦ(m),stml≤m≤mu, Where ml and mu are the upper and lower limits of the model parameters, respectively, reflecting the known geological constraints; The Bayesian-Markov Chain Monte Carlo method is used for inversion. Given the observed data d, the posterior probability p(m|d) of the goaf model m is proportional to the product of the likelihood function p(d|m) and the prior probability p(m), that is: p(m|d)∝p(d|m)p(m); where the likelihood function p(d|m) reflects the degree of fit between the model m and the observed data d, and the prior probability p(m) reflects the physical constraints and statistical characteristics of the model parameters. By constructing a Markov chain to realize random sampling of the posterior probability distribution, the optimal estimate and confidence interval of the goaf model parameters are obtained. The process is as follows: S1. Based on geological data and empirical judgment, an initial model of the goaf is constructed as the starting point of the Markov chain; S2. Based on the Metropolis-Hastings criterion, the current model mi is perturbed to the next state mi+1 with probability p; the expression of the transition probability p is: p=min{1, p(mi+1|d)q(mi|mi+1)p(mi|d)q(mi+1|mi)}, Where q(mi+1|mi) is the transfer kernel function of the Markov chain, and a symmetric distribution such as Gaussian distribution is usually selected; S3, repeat step S2, continuously update and adjust the state of the Markov chain until the convergence condition is reached: the posterior probability does not change significantly; S4. Based on the stationary distribution of Markov chain, a series of inversion model samples are randomly selected and statistically analyzed to obtain the optimal estimate and confidence interval of the goaf model parameters, and the three-dimensional distribution of resistivity and seismic wave velocity in the goaf and the quantitative evaluation results of their uncertainty.

Images