Goaf joint detection and identification method fusing gravity, earthquake and transient electromagnetism

Abstract

The invention relates to the technical field of coal mine goaf exploration, in particular to a goaf combined detection and identification method fusing gravity, earthquake and transient electromagnetism, which comprises the following steps of: (1) arranging an earthquake, transient electromagnetism and gravity method combined observation system according to a detection geological task and a field condition of a coal goaf; (2) seismic data, transient electromagnetic data and gravity data are collected through the combined observation system; (3) performing independent inversion on the gravity data, the transient electromagnetic data and the seismic data to obtain density distribution, resistivity distribution and wave velocity distribution of the goaf; (4) carrying out a joint inversion method based on improved cross gradient seismic data, transient electromagnetic data and gravity data; (5) obtaining a joint inversion resistivity profile map, a joint inversion wave velocity profile map and a joint inversion density profile map of the goaf through joint inversion of the seismic data, the transient electromagnetic data and the gravity data of the goaf; and identifying development characteristics of the goaf.

Description

Technical Field

[0001] This invention relates to the field of coal mine goaf exploration technology, and particularly to a method for precise detection of geological anomalies such as coal mine goaf and underground cavities, which can provide accurate geological data support for safe coal mine production. In particular, it relates to a method for joint detection and identification of goaf by integrating gravity, seismic and transient electromagnetic methods. Background Technology

[0002] During coal mining, the existence of goaf areas can seriously threaten the safety of mine production and easily lead to geological disasters such as roof collapse, gas outburst, and water damage. Therefore, accurate detection and identification of goaf areas is one of the core tasks of coal mine geological work.

[0003] Currently, the core technology for detecting goaf areas is geophysical exploration technology. With its advantages of being non-destructive, efficient, economical, and repeatable, it is widely used in practice. Commonly used methods include seismic methods, gravity methods, transient electromagnetic methods, high-density resistivity methods, magnetotelluric sounding, and ground-penetrating radar methods.

[0004] To improve detection results, the industry has gradually developed a comprehensive geophysical exploration method, which involves using multiple geophysical methods simultaneously and comparing and interpreting the results to reduce the limitations of a single method through mutual verification.

[0005] Traditional methods for detecting and interpreting goaf areas still mainly rely on single-method inversion and multi-result comparison. This involves separately inverting gravity, seismic, and transient electromagnetic data, and then manually comparing the results of different methods to infer the characteristics of the goaf area. Although some integrated detection methods attempt joint inversion, they mostly use data splicing or simple weighting methods and have not established an effective multi-physics coupling mechanism.

[0006] Taking the transient electromagnetic detection method for water-filled goaf in steeply inclined coal seams (patent application number CN202410259639.5) as an example, this patent achieves precise positioning of deep water-filled goaf areas by densely laying transmitting source wires on the ground and receiving measurements in underground roadways. However, this method only improves upon transient electromagnetic methods and does not involve the combined application of multiple geophysical methods. Therefore, it remains difficult to comprehensively and accurately detect various characteristics of goaf areas under complex geological conditions.

[0007] The aforementioned literature and existing technologies have the following main drawbacks:

[0008] First, single-method inversion suffers from significant ambiguity and insufficient detection accuracy. Seismic methods have strong vertical resolution but weak horizontal resolution and cannot distinguish whether mined-out areas contain water; transient electromagnetic methods are sensitive to low resistivity but have poor vertical resolution and high low-frequency noise; gravity methods have high lateral resolution but low depth resolution. Single-method inversion or simple comparison results cannot resolve ambiguity, resulting in low data utilization and difficulty in leveraging the advantages of multiple methods.

[0009] Second, traditional joint exploration lacks an effective physical field coupling mechanism. It often uses data stitching or simple weighting, ignoring the structural correlation of physical property parameters such as density, wave velocity, and resistivity, and fails to address the issue of differences in the order of magnitude of different physical property parameters. This leads to conflicting inversion results, poor convergence, and increases the ambiguity of geological interpretation and the risks of engineering decision-making.

[0010] Therefore, it is necessary to design a joint detection and identification method for goaf areas that integrates gravity, seismic and transient electromagnetic methods. Summary of the Invention

[0011] To solve one of the above-mentioned technical problems, the technical solution adopted by the present invention is: a method for joint detection and identification of goaf areas integrating gravity, earthquake and transient electromagnetic methods, including the following steps: (1) According to the geological task of coal goaf detection and the site conditions, a joint observation system of earthquake, transient electromagnetic and gravity methods is set up;

[0012] (2) Seismic data, transient electromagnetic data and gravity data are acquired through the joint observation system;

[0013] (3) Perform separate inversions on gravity data, transient electromagnetic data, and seismic data to obtain the density distribution, resistivity distribution, and wave velocity distribution of the goaf.

[0014] (4) A joint inversion method based on improved cross-gradient seismic data, transient electromagnetic data, and gravity data is implemented:

[0015] ① Establish the objective function for the joint inversion of gravity, seismicity, and electrical systems with the addition of a normalization operator for the cross gradient:

[0016] ② Construct cross-gradient terms for the preliminary inversion results of gravity data, seismic data, and transient electromagnetic data;

[0017] ③ Perform joint inversion of gravity data, electrical resistivity data, and seismic data based on cross-gradient structural consistency;

[0018] (5) By jointly inverting the seismic data, transient electromagnetic data and gravity data of the goaf, the joint inverted resistivity profile, joint inverted wave velocity profile and joint inverted density profile of the goaf are obtained, and the development characteristics of the goaf are identified.

[0019] Based on any of the above technical solutions, the following optimization is made: In step (1), according to the geological task of coal goaf detection and the site conditions, the following is specifically made: the estimated depth, scale and strata lithology (such as the coal-bearing strata of the Carboniferous system and the limestone of the Liangjiashan Formation of the Ordovician system) of the target goaf in the detection area are clearly defined, while avoiding areas that affect the quality of data acquisition, such as surface high-voltage lines and industrial interference areas; the joint observation system must ensure that the survey lines of the three methods of seismic, transient electromagnetic and gravity are consistent and the positions of the measuring points correspond one by one to ensure the spatial matching of the collected data.

[0020] Based on any of the above technical solutions, the following further optimizations are made: The specific parameters of the transient electromagnetic method observation system in step (1) are as follows: the measurement network density is set to 40m×10m, that is, the measurement line spacing is 40m and the measurement point spacing is 10m; a center loop device is adopted, the side length of the transmitting frame is 240m, the receiving coil is placed in the center 1 / 4 area of ​​the transmitting frame, and the number of turns of the receiving coil is adjusted according to the detection depth to ensure effective coverage of the target goaf area; the transmission frequency parameter is set to 25Hz, the transmission current stability is controlled within ±1%, and the turn-off time is ≤50μs to avoid the turn-off time being too long and affecting the late signal acquisition.

[0021] Based on any of the above technical solutions, the following further optimization is made: In step (1), the seismic observation system is a three-dimensional seismic observation system with the following specific parameters: a 12-line, 6-shot observation method is adopted, with the excitation point located at the midpoint of the survey line; the shot point spacing is 10m, the shot row spacing is 40m, the receiver line spacing is 40m, and the channel spacing is 10m; each receiver line is equipped with 64 detectors, and the detectors are piezoelectric detectors with a main frequency of 30Hz to ensure sensitive reception of seismic wave signals; the sampling rate is set to 1ms, and the recording length is ≥2s to ensure that the reflected wave signals corresponding to the target goaf area can be collected; the underground CDP (common depth point) element size is 5m×5m, and the full coverage number reaches 48 times, which meets the resolution requirements for fine detection of the goaf area.

[0022] Based on any of the above technical solutions, the following further optimizations are made: The specific parameters of the gravity method observation system in step (1) are as follows: The measurement network density is consistent with that of the transient electromagnetic method, set at 40m×10m; a high-precision gravimeter is used, with an instrument resolution ≤0.1μGal and an observation accuracy ≤5μGal; static zero-point drift correction is required before observation, with a correction time ≥24h; during the observation process, after completing the observation of 20 measurement points, the reference point must be returned for repeated observation, with the mean error of repeated observation ≤3μGal; the elevation of each measurement point is recorded synchronously during observation, and the elevation measurement is performed using a total station with a measurement accuracy ≤0.1m, providing basic data for the subsequent topographic correction of gravity data.

[0023] Based on any of the above technical solutions, the following further optimization is made: the specific operational requirements for data acquisition in step (2) are as follows:

[0024] Seismic data acquisition: The excitation method uses explosives. The amount of explosives is adjusted according to the lithology of the strata in the detection area. In soft strata, the amount of explosives is controlled at 0.5-1kg, and in hard strata, the amount of explosives is controlled at 1-2kg to ensure effective excitation of seismic waves and avoid excessive interference with the surrounding environment. The geophone needs to be tightly coupled to the ground with a coupling coefficient ≥0.9 to reduce signal attenuation caused by poor contact between the geophone and the ground.

[0025] Transient electromagnetic data acquisition: Each measuring point is observed three times consecutively, and the average value is taken as the final data of the measuring point after removing outliers; the signal-to-noise ratio is monitored in real time during the observation to ensure that the signal-to-noise ratio of the late signal (corresponding to the deep goaf) is ≥5. If the signal-to-noise ratio does not meet the requirements, the observation parameters need to be adjusted or the observation needs to be repeated.

[0026] Gravity data acquisition: Three sets of data were observed at each measuring point, with each set of observations lasting ≥60s. After removing outlier data sets through data stability analysis, the average was taken. During the observation period, vehicles were strictly prohibited from passing within 50m of the measuring point, and personnel were prevented from walking near the measuring point to reduce the impact of vibration interference on the observation data.

[0027] Based on any of the above technical solutions, the following optimization is made: The specific process of the individual inversion in step (3) is as follows:

[0028] Gravity data inversion: First, the original gravity data is sequentially corrected by topography, Bouguer, and normal field to eliminate the influence of factors such as topographic relief and crustal density differences. Then, a three-dimensional apparent density inversion algorithm is used, with the resolution of the inversion grid matched to the measurement grid density, i.e., the horizontal grid spacing is 40m×10m, and the vertical grid spacing is adjusted according to the detection depth, 5m for shallow layers and 20m for deep layers. The initial model is set based on regional geological data, such as the normal stratum density is set to 2.6g / cm³, and the density of the goaf is initially estimated to be 2.0-2.2g / cm³. The inversion iteration is ≥10 times until the data fitting difference is ≤5%, and the three-dimensional density distribution model of the goaf is output.

[0029] Transient electromagnetic data inversion: The raw data is first filtered by power frequency (50Hz) and trend surface filtering to eliminate industrial interference and regional trend interference caused by surface topography; then topographic correction and device correction are performed to correct the influence of observation device and topographic undulation on the data; a one-dimensional Occam inversion algorithm is used, with the number of inversion layers set to 15-20, and the layer thickness gradually increasing from shallow to deep, 5m for shallow layers and 20m for deep layers; during the inversion process, the least squares criterion is used to fit the observation data, and the vertical resistivity distribution at each measuring point and the planar resistivity distribution of the entire detection area are output.

[0030] Seismic data inversion: The raw seismic data is first preprocessed, including static correction (eliminating the influence of surface elevation differences), denoising (suppressing random noise and surface wave interference), and amplitude compensation (compensating for energy attenuation of seismic waves during propagation); then dynamic correction and horizontal stacking are performed to obtain a stacked profile; based on the stacked profile, the reflection phase axis of the target layer (such as coal seam No. 6) is picked, and a post-stack wave velocity inversion algorithm is used, with the inversion interval consistent with the CDP surface element, which is 5m×5m; a three-dimensional wave velocity distribution model of the entire detection area is output, in which the wave velocity of normal coal seams is set to 2800-3000m / s, and the wave velocity of goaf is initially estimated to be 2200-2500m / s.

[0031] Based on any of the above technical solutions, a further optimization is made in step 4①: establishing a joint inversion objective function for gravity-seismic-electrical gradients with the addition of a normalization operator.

[0032] ;

[0033] Constraints: ;

[0034] Where, d i Represents observational data; m i Represents model parameters; m i _ref represents the reference model; f i (m i ) represents the model parameter m i The forward response, α i Represents the regularization factor; C di Representative observation data d i Data covariance matrix; C mi Representative observation data m i The model covariance matrix; ▽ represents the gradient, and τ represents the cross gradient constraint term.

[0035] Based on any of the above technical solutions, a further optimization is made to the cross-gradient term in step 4② of constructing the preliminary inversion results of gravity data, seismic data, and transient electromagnetic data:

[0036] The product of the gradient vectors of two physical property parameters can reflect the structural similarity between the two models. The cross gradient function of the two physical property parameters is defined as follows: ;

[0037] Based on the inversion results of different physical properties, the cross gradient in the structurally consistent region is close to 0, while the cross gradient magnitude is large in the structurally inconsistent region. The cross gradient term is added to the objective function to constrain the inversion results to tend towards structural consistency, thereby outputting a structurally consistent physical property distribution.

[0038] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0039] 1. Currently, most comprehensive geophysical exploration data of mining subsidence areas are processed using single-method inversion or by comparing and interpreting the results of different exploration methods. There is a lack of true multi-method joint exploration and multi-field data joint interpretation. There is also a lack of research on joint exploration technology that integrates gravity, seismic, and transient electromagnetic data.

[0040] This invention integrates gravity, seismic, and transient electromagnetic data, and performs multi-property parameter joint inversion on the three types of data to analyze the detection effect. It leverages the advantages of comprehensive detection using multiple geophysical methods and provides data support for the selection of detection methods for coal mining subsidence areas in the future.

[0041] 2. The magnitudes and units of different physical properties of underground media vary greatly (for example, the magnitudes of velocity, resistivity, and density of underground media differ significantly). If different physical properties are directly coupled, the joint inversion results of individual properties may fail to obtain structural similarity constraints transferred from other properties, and the joint inversion results will not be significantly improved compared to the individual inversion results.

[0042] Based on the above ideas, in this invention we improve the traditional cross gradient function by adding a normalization operator to the cross gradient function to overcome the influence of differences between different physical property parameters.

[0043] The main difference between multi-parameter joint inversion and dual-parameter joint inversion is that the cross-gradient function constraint term increases from one to multiple terms. That is, structural similarity constraints are required between each physical property parameter. Therefore, the problem of multiple solutions in inversion can be reduced to a greater extent, and the obtained similar structural model is more conducive to being converted into a geological model, thus helping geologists in their interpretation work.

[0044] 3. Significantly reduces inversion ambiguity and improves the accuracy of goaf detection:

[0045] This invention addresses the issue of result biases that easily occur in single-method or dual-property joint inversion by employing a triangular verification mechanism (density-wave velocity, density-resistivity, and wave velocity-resistivity cross-gradient constraints) based on three physical property parameters: gravity, seismic, and transient electromagnetic. For example, joint inversion can correct the shortcomings of inaccurate time-depth transformation in seismic inversion, insufficient depth resolution in gravity inversion, and poor lateral continuity in transient electromagnetic inversion, thus accurately matching actual geological conditions.

[0046] 4. Achieve deep integration of multiple data sources to improve data utilization:

[0047] Compared to the traditional single-method inversion followed by result comparison, this invention eliminates the order-of-magnitude differences in density, wave velocity, and resistivity by incorporating a normalization operator into the joint inversion objective function. This allows the three types of data to mutually transmit structural constraints during the inversion process, rather than simply splicing or weighting them. Simultaneously, the linearization of nonlinear problems is achieved through Taylor series expansion and the Lagrange operator method, fully leveraging the geological information value of each data type and avoiding information waste caused by noise and interference from individual data. This results in improved data utilization compared to traditional methods.

[0048] 5. Optimize the ability to identify features of goaf areas, taking into account both morphology and water content:

[0049] The combined resistivity, wave velocity, and density profiles output by this invention can identify goaf areas using a quantitative standard of ≥80% overlap of the three low anomalies (low resistivity, low wave velocity, and low density). It can also distinguish whether the goaf contains water (resistivity in water-bearing goaf areas is ≥50% lower than in dry goaf areas). For example, at a horizontal range of 2150-2520 meters and a depth of 40-80 meters, the combined inversion can clearly depict the layered / capillary morphology of the goaf, and the continuity characteristics of the high-resistivity layer can reflect the stability of the upper strata of the goaf, providing a more comprehensive geological basis for subsequent engineering decisions.

[0050] 6. Enhance adaptability to complex geological conditions and broaden exploration scenarios:

[0051] For typical coal mine geological conditions such as Carboniferous coal-bearing strata and Ordovician limestone, this invention optimizes observation system parameters (such as a 40m×10m transient electromagnetic network, full coverage of 48 seismic events, and high-precision gravity correction) and inversion algorithms (3D apparent density inversion, 1D Occam inversion, and post-stack wave velocity inversion). This effectively avoids the influence of surface high-voltage lines, industrial interference, and strata fracture zones. Compared to single methods that can only adapt to simple geological scenarios, this invention can cover complex scenarios such as deep goaf areas (≥80m), waterlogged goaf areas, and multi-layered goaf areas, significantly broadening its applicability.

[0052] 7. Establish an evaluation-feedback-optimization closed loop to ensure the practicality of the project:

[0053] This invention adds an accuracy evaluation step, using known geological data of the goaf as a benchmark to quantify the errors in horizontal range, vertical depth, and physical parameters. If the results are unsatisfactory, the weighting factors can be adjusted to re-invert the data, forming a closed-loop optimization mechanism. This mechanism avoids the drawback of traditional methods that involve one-time inversion without verification, ensuring that the detection results meet the requirements of accurate goaf location for safe coal mine production. It can be directly used for goaf remediation scheme formulation and has strong engineering practicality.

[0054] 8. Balancing detection efficiency and cost, it has the potential for large-scale application:

[0055] In this invention, the measurement network density (e.g., 40m × 10m for both transient electromagnetic and gravity methods) and observation procedures are mutually matched, eliminating the need for separate and repeated measurement network deployments. Compared to conducting detection using the three methods separately, this can save more than 30% of on-site operation time and labor costs. Furthermore, the inversion algorithm reduces computational complexity through linearization and iterative convergence control, allowing data processing to be completed using ordinary computers without relying on high-end computing equipment, thus facilitating large-scale application in small and medium-sized coal mining enterprises. Attached Figure Description

[0056] To more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the accompanying drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. In all the drawings, similar elements or components are generally identified by similar reference numerals. In the drawings, the elements or components are not necessarily drawn to scale.

[0057] Figure 1 Flowchart of a method for identifying coal mine goaf based on joint detection of gravity, seismic and transient electromagnetic data and joint inversion of multiphysics field data.

[0058] Figure 2 : A diagram of the joint inversion system architecture that integrates gravity, seismic and transient electromagnetic multiphysics data.

[0059] Figure 3 Flowchart of the multiphysics joint inversion algorithm.

[0060] Figure 4 Resistivity profiles of goaf areas obtained by electrical resistivity methods alone and in combination.

[0061] Figure 5 Seismic inversion of goaf wave velocity profiles (individual and combined).

[0062] Figure 6 : Gravity-based individual and combined inversion density profiles of goaf areas.

[0063] In the figure, (a) is the result of individual inversion; (b) is the result of joint inversion. Detailed Implementation

[0064] The embodiments of the technical solution of the present invention will now be described in detail with reference to the accompanying drawings. The following embodiments are only used to more clearly illustrate the technical solution of the present invention, and are therefore merely examples and should not be used to limit the scope of protection of the present invention. The specific structure of the present invention is as follows: Figures 1-6 As shown in the image.

[0065] Example 1: A method for joint detection and identification of goaf areas integrating gravity, seismic and transient electromagnetic methods, including the following steps: (1) Based on the geological task of coal goaf detection and the site conditions, deploy a joint observation system of seismic, transient electromagnetic and gravity methods;

[0066] (2) Seismic data, transient electromagnetic data and gravity data are acquired through the joint observation system;

[0067] (3) Perform separate inversions on gravity data, transient electromagnetic data, and seismic data to obtain the density distribution, resistivity distribution, and wave velocity distribution of the goaf.

[0068] (4) A joint inversion method based on improved cross-gradient seismic data, transient electromagnetic data, and gravity data is implemented:

[0069] ① Establish the objective function for the joint inversion of gravity, seismicity, and electrical systems with the addition of a normalization operator for the cross gradient:

[0070] ② Construct cross-gradient terms for the preliminary inversion results of gravity data, seismic data, and transient electromagnetic data;

[0071] ③ Perform joint inversion of gravity data, electrical resistivity data, and seismic data based on cross-gradient structural consistency;

[0072] (5) By jointly inverting the seismic data, transient electromagnetic data and gravity data of the goaf, the joint inverted resistivity profile, joint inverted wave velocity profile and joint inverted density profile of the goaf are obtained, and the development characteristics of the goaf are identified.

[0073] Based on any of the above technical solutions, the following optimization is made: In step (1), according to the geological task of coal goaf detection and the site conditions, the following is specifically made: the estimated depth, scale and strata lithology (such as the coal-bearing strata of the Carboniferous system and the limestone of the Liangjiashan Formation of the Ordovician system) of the target goaf in the detection area are clearly defined, while avoiding areas that affect the quality of data acquisition, such as surface high-voltage lines and industrial interference areas; the joint observation system must ensure that the survey lines of the three methods of seismic, transient electromagnetic and gravity are consistent and the positions of the measuring points correspond one by one to ensure the spatial matching of the collected data.

[0074] It should be noted that by clearly identifying the geological features and sources of interference in the detection area in advance, external interference in data acquisition can be avoided in a targeted manner. At the same time, the survey lines and points of the three methods can be accurately matched, so that the subsequently acquired seismic, transient electromagnetic, and gravity data have the same spatial coordinate reference. This avoids deviations in the inversion results caused by spatial misalignment, lays a spatial consistency foundation for multi-data joint inversion, and improves the accuracy of subsequent data fusion.

[0075] Based on any of the above technical solutions, the following further optimizations are made: The specific parameters of the transient electromagnetic method observation system in step (1) are as follows: the measurement network density is set to 40m×10m, that is, the measurement line spacing is 40m and the measurement point spacing is 10m; a center loop device is adopted, the side length of the transmitting frame is 240m, the receiving coil is placed in the center 1 / 4 area of ​​the transmitting frame, and the number of turns of the receiving coil is adjusted according to the detection depth to ensure effective coverage of the target goaf area; the transmission frequency parameter is set to 25Hz, the transmission current stability is controlled within ±1%, and the turn-off time is ≤50μs to avoid the turn-off time being too long and affecting the late signal acquisition.

[0076] It should be noted that: the 40m×10m measurement network density can balance detection efficiency and resolution, covering a large area while capturing the fine features of the goaf; the combination of the central loop device and the 240m transmitting frame can enhance the propagation depth and intensity of transient electromagnetic signals, and the reception in the central 1 / 4 area can reduce edge effect interference; the combination of 25Hz transmission frequency and ≤50μs off time can effectively collect late signals corresponding to deep goaf areas (late signals reflect the characteristics of deep geological bodies), and the current stability control can reduce signal noise, ensuring the sensitivity and accuracy of transient electromagnetic data for goaf areas (especially waterlogged goaf areas).

[0077] Based on any of the above technical solutions, the following further optimization is made: In step (1), the seismic observation system is a three-dimensional seismic observation system with the following specific parameters: a 12-line, 6-shot observation method is adopted, with the excitation point located at the midpoint of the survey line; the shot point spacing is 10m, the shot row spacing is 40m, the receiver line spacing is 40m, and the channel spacing is 10m; each receiver line is equipped with 64 detectors, and the detectors are piezoelectric detectors with a main frequency of 30Hz to ensure sensitive reception of seismic wave signals; the sampling rate is set to 1ms, and the recording length is ≥2s to ensure that the reflected wave signals corresponding to the target goaf area can be collected; the underground CDP (common depth point) element size is 5m×5m, and the full coverage number reaches 48 times, which meets the resolution requirements for fine detection of the goaf area.

[0078] It should be noted that: the 12-line 6-shot pattern and 48 full-coverage times can suppress random noise in the propagation of seismic waves through multiple superpositions, thereby enhancing the effective reflection signal; the 5m×5m CDP surface cells can improve the lateral resolution and accurately locate the horizontal range of the goaf; the 30Hz main frequency detector paired with a 1ms sampling rate can sensitively capture the wave velocity difference signal between the goaf and normal strata (the wave velocity in the goaf is lower than that in the normal strata), and the ≥2s recording length can ensure the acquisition of reflected waves from deep goaf areas, providing high-quality data for subsequent wave velocity inversion and ensuring the accuracy of vertical depth identification of the goaf.

[0079] Based on any of the above technical solutions, the following further optimizations are made: The specific parameters of the gravity method observation system in step (1) are as follows: The measurement network density is consistent with that of the transient electromagnetic method, set at 40m×10m; a high-precision gravimeter is used, with an instrument resolution ≤0.1μGal and an observation accuracy ≤5μGal; static zero-point drift correction is required before observation, with a correction time ≥24h; during the observation process, after completing the observation of 20 measurement points, the reference point must be returned for repeated observation, with the mean error of repeated observation ≤3μGal; the elevation of each measurement point is recorded synchronously during observation, and the elevation measurement is performed using a total station with a measurement accuracy ≤0.1m, providing basic data for the subsequent topographic correction of gravity data.

[0080] It should be noted that: a measurement network density consistent with that of the transient electromagnetic method ensures spatial matching between gravity data and transient electromagnetic data, facilitating subsequent joint inversion; a high-precision gravimeter with a resolution of ≤0.1μGal and an observation accuracy of ≤5μGal can capture minute density differences between goaf areas and normal strata (the density of goaf areas is lower than that of normal strata); static zero-point drift correction and repeated observations can eliminate the influence of instrument errors and environmental changes (such as temperature and air pressure) on the observation data; and elevation data with an accuracy of ≤0.1m provides an accurate basis for terrain correction, avoiding misjudgments of gravity anomalies caused by terrain undulations, and ensuring that gravity data can accurately reflect the density anomaly characteristics of goaf areas.

[0081] Based on any of the above technical solutions, the following further optimization is made: the specific operational requirements for data acquisition in step (2) are as follows:

[0082] Seismic data acquisition: The excitation method uses explosives. The amount of explosives is adjusted according to the lithology of the strata in the detection area. In soft strata, the amount of explosives is controlled at 0.5-1kg, and in hard strata, the amount of explosives is controlled at 1-2kg to ensure effective excitation of seismic waves and avoid excessive interference with the surrounding environment. The geophone needs to be tightly coupled to the ground with a coupling coefficient ≥0.9 to reduce signal attenuation caused by poor contact between the geophone and the ground.

[0083] Transient electromagnetic data acquisition: Each measuring point is observed three times consecutively, and the average value is taken as the final data of the measuring point after removing outliers; the signal-to-noise ratio is monitored in real time during the observation to ensure that the signal-to-noise ratio of the late signal (corresponding to the deep goaf) is ≥5. If the signal-to-noise ratio does not meet the requirements, the observation parameters need to be adjusted or the observation needs to be repeated.

[0084] Gravity data acquisition: Three sets of data were observed at each measuring point, with each set of observations lasting ≥60s. After removing outlier data sets through data stability analysis, the average was taken. During the observation period, vehicles were strictly prohibited from passing within 50m of the measuring point, and personnel were prevented from walking near the measuring point to reduce the impact of vibration interference on the observation data.

[0085] It should be noted that: Seismic data acquisition, by adjusting the explosive charge to suit different lithologies, ensures the generation of sufficiently strong seismic waves (ensuring deep signal propagation) while avoiding excessive interference; a high coupling coefficient reduces signal loss between the geophone and the ground, ensuring complete reception of reflected wave signals. Averaging three transient electromagnetic observations and controlling the signal-to-noise ratio reduce the impact of random noise on later signals, ensuring the validity of signals from deep goaf areas. Multiple observations, long-term acquisition, and interference avoidance for gravity data reduce gravity value fluctuations caused by vibrations (such as vehicle and personnel activity); stability analysis eliminates outliers, ensuring the authenticity of gravity data and providing reliable raw data support for subsequent inversion.

[0086] Based on any of the above technical solutions, the following optimization is made: The specific process of the individual inversion in step (3) is as follows:

[0087] Gravity data inversion: First, the original gravity data is sequentially corrected by topography, Bouguer, and normal field to eliminate the influence of factors such as topographic relief and crustal density differences. Then, a three-dimensional apparent density inversion algorithm is used, with the resolution of the inversion grid matched to the measurement grid density, i.e., the horizontal grid spacing is 40m×10m, and the vertical grid spacing is adjusted according to the detection depth, 5m for shallow layers and 20m for deep layers. The initial model is set based on regional geological data, such as the normal stratum density is set to 2.6g / cm³, and the density of the goaf is initially estimated to be 2.0-2.2g / cm³. The inversion iteration is ≥10 times until the data fitting difference is ≤5%, and the three-dimensional density distribution model of the goaf is output.

[0088] Transient electromagnetic data inversion: The raw data is first filtered by power frequency (50Hz) and trend surface filtering to eliminate industrial interference and regional trend interference caused by surface topography; then topographic correction and device correction are performed to correct the influence of observation device and topographic undulation on the data; a one-dimensional Occam inversion algorithm is used, with the number of inversion layers set to 15-20, and the layer thickness gradually increasing from shallow to deep, 5m for shallow layers and 20m for deep layers; during the inversion process, the least squares criterion is used to fit the observation data, and the vertical resistivity distribution at each measuring point and the planar resistivity distribution of the entire detection area are output.

[0089] Seismic data inversion: The raw seismic data is first preprocessed, including static correction (eliminating the influence of surface elevation differences), denoising (suppressing random noise and surface wave interference), and amplitude compensation (compensating for energy attenuation of seismic waves during propagation); then dynamic correction and horizontal stacking are performed to obtain a stacked profile; based on the stacked profile, the reflection phase axis of the target layer (such as coal seam No. 6) is picked, and a post-stack wave velocity inversion algorithm is used, with the inversion interval consistent with the CDP surface element, which is 5m×5m; a three-dimensional wave velocity distribution model of the entire detection area is output, in which the wave velocity of normal coal seams is set to 2800-3000m / s, and the wave velocity of goaf is initially estimated to be 2200-2500m / s.

[0090] It should be noted that: gravity data undergoes multiple rounds of correction to eliminate interference from non-mining-out area factors (topography, normal crustal density); 3D apparent density inversion and matching the grid resolution of the survey network accurately characterize the density distribution of the mining-out area; initial model and iterative fitting difference control ensure that the inversion results closely match the actual geological conditions. Transient electromagnetic data undergoes filtering and correction to eliminate industrial and topographic interference; the layered design of the one-dimensional Occam inversion clearly reflects resistivity changes at different depths; the least squares criterion improves data fitting accuracy and accurately identifies low-resistivity anomalies (waterlogged mining-out areas). Seismic data preprocessing and overlay processing enhance effective signals; post-stack wave velocity inversion and CDP surface matching accurately capture the low wave velocity characteristics of the mining-out area; target layer picking anchors the layer where the mining-out area is located, providing an independent and reliable distribution model of various physical parameters for subsequent joint inversion.

[0091] Based on any of the above technical solutions, a further optimization is made in step 4①: establishing a joint inversion objective function for gravity-seismic-electrical gradients with the addition of a normalization operator.

[0092] ;

[0093] Constraints: ;

[0094] Where, d i Represents observational data; m i Represents model parameters; m i _ref represents the reference model; f i (m i ) represents the model parameter m i The forward response, α i Represents the regularization factor; C di Representative observation data d i Data covariance matrix; C mi Representative observation data m i The model covariance matrix; ▽ represents the gradient, and τ represents the cross gradient constraint term.

[0095] The magnitudes and units of different physical properties of underground media vary significantly (for example, the magnitudes of velocity, resistivity, and density differ considerably). Directly coupling different physical properties may result in the joint inversion results for individual properties failing to capture the structural similarity constraints passed down to other properties, ultimately failing to significantly improve the joint inversion results compared to individual inversion results. Based on this idea, we improve the traditional cross-gradient function by adding a normalization operator to overcome the influence of differences between different physical properties.

[0096] It should be noted that the joint inversion objective function, through the combination of data fitting terms (the first term) and model constraint terms (the second term), ensures the consistency between the inversion results and the observed data while preventing the model from deviating excessively from the reference model (a reasonable initial model based on geological data). The cross-gradient constraint term (the third term) promotes structural consistency among different physical property models. The core advantage of adding the normalization operator lies in solving the constraint imbalance problem caused by the difference in magnitude of density (g / cm³), wave velocity (m / s), and resistivity (Ω·m). This allows the three physical property parameters to mutually transfer structural constraints in the joint inversion, preventing a single physical property parameter from being ignored due to its small magnitude, ensuring the structural consistency and reliability of the joint inversion results, and improving the multi-data fusion effect.

[0097] Based on any of the above technical solutions, a further optimization is made by adding a normalization operator to the cross gradient function expression, as follows:

[0098] ;

[0099] Where, m i Let k be the vector of the i-th physical property parameter. i Let be the normalization operator for the i-th physical property parameter. For the model trials in this paper, the normalization operator can be expressed as:

[0100] ;

[0101] Max and Min represent the maximum and minimum values, respectively. This is the individual inversion model obtained by the i-th property method.

[0102] For convenience, we will integrate the model vectors. , .

[0103] Then, the objective function It can be rewritten as:

[0104] ;

[0105] Constraints: ;

[0106] in:

[0107] ;

[0108] , ,

[0109] , , .

[0110] It should be noted that the normalization operator, by mapping the physical property parameters of individual inversion models to the [0,1] interval, eliminates the order-of-magnitude differences between different physical property parameters (such as the difference between wave velocities of thousands of m / s and resistivity of hundreds of Ω·m), enabling the cross-gradient function to fairly measure the structural similarity of each physical property model. The integrated model vector and objective function simplify the calculation process for multiple physical property parameters. At the same time, through the weight factors and covariance matrix in diagonal matrix form, the weights can be adjusted according to the quality of different data (such as signal-to-noise ratio and observation accuracy), allowing high-quality data to play a greater role in the inversion and further improving the accuracy and stability of the joint inversion results.

[0111] Based on any of the above technical solutions, the following optimization is made: The main difference between multi-parameter joint inversion and dual-parameter joint inversion is that the cross gradient function constraint term is increased from one to multiple terms. That is, structural similarity constraints are required between each physical property parameter. Therefore, the problem of multiple solutions in inversion can be reduced to a greater extent, and the obtained similar structural model is more conducive to being converted into a geological model, thus helping geologists in their interpretation work.

[0112] It should be noted that dual-property joint inversion can only establish structural constraints (such as density and wave velocity) between two properties, and may still lead to deviations in the inversion results due to the lack of verification by a third property (such as resistivity). In contrast, multi-property (gravity-seismic-transient electromagnetic) joint inversion, through cross-gradient constraints of density and wave velocity, density and resistivity, and wave velocity and resistivity, forms a triangular verification mechanism. This allows the three property models to mutually constrain and verify each other, significantly reducing the problem of multiple solutions under a single constraint. This multi-constraint mechanism allows the inversion results to more closely resemble the actual geological structure. When subsequently converted into a geological model, it can more clearly reflect the boundaries, depths, and other characteristics of the goaf, reducing the difficulty and error of geological interpretation.

[0113] Building upon any of the aforementioned technical solutions, a further optimization involves solving the joint inversion objective function for multiple physical properties and constraints. First, a Taylor series expansion is used to linearize the objective function and the cross-gradient constraint function, transforming the nonlinear problem into a linear one. Then, the constrained problem is transformed into an unconstrained problem using the Lagrange operator method (Tarantola, 1987), allowing the cross-gradient constraint function to be incorporated into the objective function.

[0114] ;

[0115] in, , .

[0116] ;

[0117] ;

[0118] Where Λ represents the Lagrange operator; A represents the Jacobian matrix; and B represents the normalized cross-gradient constraint function. The partial derivatives of .

[0119] Find the extrema of the objective function (3.8). ,

[0120] ;

[0121] ;

[0122] Combining the above two equations, we can obtain the expression for the Lagrange operator:

[0123] ;

[0124] ;

[0125] ;

[0126] ;

[0127] The final model change is:

[0128] .

[0129] It should be noted that: Taylor series expansion transforms the complex nonlinear objective function into an easily solvable linear form, reducing computational difficulty; the Lagrange operator rule, by introducing an operator, incorporates the constraint that the inversion result must fit the observed data into the objective function, avoiding the computational complexity caused by handling constraint problems separately. The final derived formula for model change can be used to gradually adjust model parameters based on the differences between the current iterative model and the reference model, the observed data, and the influence of cross-gradient constraints. This ensures that each iteration moves closer to a direction with optimal data fit and consistent structure, improving the convergence speed and accuracy of the joint inversion.

[0130] Based on any of the above technical solutions, the following further optimization is made: in step 4②, the cross gradient term of the preliminary inversion results of gravity data, seismic data and transient electromagnetic data is constructed;

[0131] The product of the gradient vectors of two physical property parameters can reflect the structural similarity between the two models. The cross gradient function of the two physical property parameters is defined as follows: ;

[0132] Based on the inversion results of different physical properties, the cross gradient in the structurally consistent region is close to 0, while the cross gradient magnitude is large in the structurally inconsistent region. The cross gradient term is added to the objective function to constrain the inversion results to tend towards structural consistency, thereby outputting a structurally consistent physical property distribution.

[0133] It should be noted that the mathematical properties of the gradient vector cross product dictate that when two physical property models have the same structure in a certain area (e.g., both are goaf areas or both are normal strata), the gradient directions are similar, and the cross product result approaches 0. Conversely, when the structures are inconsistent (e.g., one physical property model shows a goaf area while the other shows a normal strata), the gradient directions differ significantly, and the amplitude of the cross product result increases substantially. Adding the cross gradient term to the objective function automatically minimizes the cross gradient value during the inversion process, gradually unifying the different physical property models structurally. Ultimately, the output density, wave velocity, and resistivity distributions all point to the same geological body (goaf area), improving the reliability of the detection results.

[0134] Based on any of the above technical solutions, a further optimization is made to the inverted resistivity profile obtained by the transient electromagnetic method, as shown in the figure. Figure 4 As shown, where, Figure 4The attached figure in section a) shows the resistivity profile obtained by a single transient electromagnetic method. It can be seen that, due to the use of a single method, the horizontal continuity of the inverted resistivity profile is poor. Furthermore, the vertical resistivity of the inverted profile exhibits a low-high-low variation. However, the strata in this area are typical Carboniferous coal-bearing strata. The target coal seam 6 is located in the Taiyuan Formation, beneath which lies the Ordovician Liangjiashan Formation, a coal-bearing strata basement primarily composed of limestone with relatively high resistivity. Therefore, in the absence of goafs or local geological anomalies, the overall resistivity of the strata in this area should exhibit a clear low-high variation from shallow to deep. Thus, the resistivity pattern obtained using a single method differs somewhat from the actual strata. The range of the significant low-resistivity anomaly in the resistivity profile obtained by the single method is inconsistent with the known goaf range. When using the joint inversion method, the lateral continuity of the inverted resistivity profile is significantly improved, while the vertical electrical variation pattern shows a low-to-high characteristic from shallow to deep, which is consistent with the actual strata. The joint inversion resistivity profile shows a significant low-resistivity anomaly around 2150-2520 meters, which is basically consistent with the actual goaf area. At the same time, its vertical resolution is significantly improved compared with the single-method inversion, and this high-resistivity layer is relatively continuous, indicating that the strata above the goaf are relatively stable and unbroken, which is consistent with the actual geological data.

[0135] It should be noted that the single transient electromagnetic method, due to the influence of noise and volume effects, cannot fully reflect the actual electrical properties of the strata, leading to misjudgments of the low-resistivity anomaly range. However, the joint inversion, through the structural constraints of seismic and gravity data, corrects the biases of the transient electromagnetic inversion. The high vertical resolution of seismic data helps anchor stratigraphic layers (such as the interface between the No. 6 coal seam and the Ordovician limestone), and the density anomalies of gravity data help delineate the horizontal range of the goaf. The combination of these three factors makes the lateral continuity and vertical regularity of the resistivity profile more closely match the actual geological conditions. The precise match between the low-resistivity anomaly at 2150-2520 meters and the actual goaf, as well as the continuity characteristics of the high-resistivity layer, proves that the joint inversion can effectively improve the interpretation accuracy of the resistivity profile and provide a reliable basis for goaf location.

[0136] Based on any of the above technical solutions, a further optimization is made by: Figure 5 The seismic results show that inaccurate estimation of the seismic time-depth transformation parameters in seismic inversion alone leads to discrepancies between the depth of the inverted wave velocity anomaly and the depth of the electrical resistivity inversion results. By incorporating electrical resistivity and gravity inversion results as structural constraints into the seismic inversion, the depth matching degree between the wave velocity results and the resistivity and gravity results can be effectively improved, thereby enhancing the inversion accuracy. (Comparison) Figure 4 and Figure 5It can be seen that the seismic inversion results have high accuracy in detecting the strata. There is a low wave velocity stratum at a depth of 60 meters, which corresponds well with the No. 6 coal seam. Furthermore, there is a significant low wave velocity anomaly with a length of 370 meters near 2150-2520 meters in the profile, which is within the goaf area.

[0137] It should be noted that the time-depth conversion in individual seismic inversion relies on velocity model assumptions; deviations in these assumptions can lead to inaccurate depth measurements of wave velocity anomalies. In contrast, in joint inversion, electrical resistivity (resistivity) and gravity (density) data provide independent depth benchmarks. For example, electrical resistivity data can determine the depth of coal seam No. 6 using known formation resistivity characteristics, while gravity data can verify the vertical extent of the goaf through density differences. Both data jointly constrain the seismic time-depth conversion parameters, ensuring that the depth of wave velocity anomalies remains consistent with the depths of other physical property anomalies. The precise correspondence between the low-velocity layer at 60 meters and coal seam No. 6, and the clear representation of the 370-meter low-velocity anomaly at 2150-2520 meters, demonstrate that joint inversion retains the high accuracy of seismic methods in identifying stratigraphic positions while correcting depth deviations, thus improving the accuracy of goaf extent and depth detection.

[0138] pass Figure 6 The gravity inversion results show that, due to insufficient gravity depth resolution, the gravity inversion results alone cannot be concentrated on the goaf area, making it impossible to determine the goaf boundary. By incorporating seismic and electrical resistivity inversion results as structural constraints into the gravity inversion, the boundary gradient change of the density results can be effectively improved, thereby increasing the inversion accuracy. The density profile in the joint inversion results indicates that the depth of the goaf is approximately 40-80 meters, located at a horizontal distance of 2100-2500 meters.

[0139] It should be noted that the physical nature of the gravity method dictates its low sensitivity to depth; when used alone, density anomalies tend to disperse, making it impossible to accurately focus on goaf areas. However, in the joint inversion, the clear stratigraphic boundaries provided by the seismic method (such as the top coal seam interface at 40 meters and the bottom coal seam interface at 80 meters) and the low resistivity anomaly boundaries provided by the electrical resistivity method serve as hard constraints for the gravity inversion, causing density anomalies to concentrate at these boundaries and enhancing boundary gradient changes. The density anomaly range at depths of 40-80 meters and horizontal distances of 2100-2500 meters highly overlaps with seismic wave velocity anomalies and electrical resistivity anomalies, demonstrating that the joint inversion effectively compensates for the insufficient depth resolution of the gravity method, achieving accurate delineation of goaf boundaries.

[0140] Based on any of the above technical solutions, a further optimization is made: In step (5), the joint inversion profile includes:

[0141] Jointly inverted resistivity profile: Low resistivity anomalies correspond to goaf areas (low resistivity is more pronounced in water-bearing goaf areas).

[0142] Jointly inverted wave velocity profile: Low wave velocity anomalies correspond to goaf areas, with an overlap of ≥80% with the boundary of low resistivity anomalies;

[0143] Joint inversion density profile: Low density anomaly corresponds to goaf area, which is consistent with the range of low resistivity and low wave velocity anomaly;

[0144] Furthermore, the joint inversion profile exhibits better lateral continuity and longitudinal resolution than the individual inversion profile.

[0145] It should be noted that the three joint inversion profiles characterize the goaf from three physical dimensions: electrical properties (resistivity), elasticity (wave velocity), and density. Due to pore development (or water accumulation), the goaf exhibits a triple-low anomaly of low resistivity, low wave velocity, and low density. Furthermore, the high overlap of these three anomalies (boundary overlap ≥80%) forms a triple verification, completely eliminating the possibility of misjudgment from a single anomaly. Simultaneously, the joint inversion, through multi-property constraints, corrects the problems of lateral discontinuity and low vertical resolution of individual inversion profiles, making the profiles more clearly reflect the morphology, extent, and depth of the goaf, providing an intuitive and reliable basis for geological interpretation.

[0146] Based on any of the above technical solutions, the following further optimization is made: the goaf identification standard is:

[0147] Resistivity: ≥30% lower than normal strata (≥50% lower in water-bearing goaf areas);

[0148] Wave velocity: ≥20% lower than normal coal seam;

[0149] Density: ≥10% lower than normal formations;

[0150] The overlap of the three anomalies is ≥80%. Record the horizontal range (±5m), vertical depth (±1m), morphology and water content of the goaf.

[0151] It should be noted that this identification standard is based on the difference in physical properties between goaf areas and normal strata—goaf areas have higher porosity than normal strata, resulting in reduced density (≥10%) and wave velocity (≥20%); if water accumulates, the pores are filled with low-resistivity water, further reducing resistivity (≥50%). The quantified percentage thresholds avoid subjective judgment errors, while the requirement of ≥80% overlap among the three anomalies ensures that the identified goaf area is a geological body pointed to by all three physical mechanisms, significantly reducing ambiguity. Simultaneously, the accuracy requirements of ±5m horizontal range and ±1m vertical depth clarify the error range of the detection results, meeting the needs of precise goaf location for safe coal mine production.

[0152] Example 2: Compared with Example 1, this example also includes the following steps:

[0153] Based on any of the above technical solutions, the following optimization is made: The method also includes an accuracy evaluation step: when comparing the joint inversion results with the known geological data of the goaf, the horizontal range error ≤10m, the vertical depth error ≤5m, and the physical property parameter error ≤15% are deemed qualified; if unqualified, return to step (4) to adjust the weight factor.

[0154] It should be noted that the accuracy evaluation step uses known geological data of the goaf (such as actual parameters of mined and verified areas) as the gold standard. It determines whether the inversion results meet engineering requirements by quantifying error indicators (horizontal range, vertical depth, and physical property parameters). This closed-loop mechanism of evaluation-feedback-optimization ensures that the final output detection results have engineering application value and provide reliable protection for safe coal mine production.

[0155] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features therein. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention. For those skilled in the art, any alternative improvements or transformations made to the implementation of the present invention fall within the protection scope of the present invention.

[0156] Any aspects of this invention not described in detail are well-known to those skilled in the art.

Claims

1. A method for detecting and identifying goaf by combining gravity, earthquake and transient electromagnetic methods, characterized in that, The method comprises the following steps: (1) according to the detection geological task and the field conditions of the coal goaf, a joint observation system of seismic, transient electromagnetic and gravity method is arranged; (2) collecting seismic data, transient electromagnetic data and gravity data through the joint observation system; (3) separately inverting the gravity data, the transient electromagnetic data and the seismic data to obtain the density distribution, the resistivity distribution and the wave velocity distribution of the goaf; (4) performing joint inversion of the seismic data, the transient electromagnetic data and the gravity data based on the improved cross gradient method: ① establishing a gravity-seismic-electric joint inversion objective function of the cross gradient with a normalization operator added: ② constructing a cross gradient term of the preliminary inversion results of the gravity data, the seismic data and the transient electromagnetic data; ③ performing joint inversion of the gravity data, the electric method data and the seismic data based on the cross gradient structure consistency; (5) obtaining the joint inversion resistivity profile, the joint inversion wave velocity profile and the joint inversion density profile of the goaf through the joint inversion of the seismic data, the transient electromagnetic data and the gravity data of the goaf, and identifying the development characteristics of the goaf.

2. The method according to claim 1, wherein, In step (1), according to the detection geological task and the field conditions of the coal goaf, the estimated depth, scale and stratum lithology (such as the carboniferous coal-bearing stratum and the ordovician liangjiashan group limestone) distribution of the target goaf in the detection area are determined, and the areas affecting the data collection quality such as the surface high-voltage line and the industrial interference area are avoided; the joint observation system needs to ensure that the measurement line directions of the three methods of seismic, transient electromagnetic and gravity are consistent, the measurement point positions are one-to-one corresponding, and the spatial matching of the collected data is ensured.

3. The method according to claim 2, wherein, In step (1), the specific parameters of the transient electromagnetic method observation system are as follows: the measurement network density is set to 40m*10m, that is, the measurement line spacing is 40m and the measurement point spacing is 10m; a central loop device is used, the outer frame length of the transmitter is 240m, the receiving coil is placed in the 1 / 4 area of the center of the transmitting outer frame, the number of turns of the receiving coil is adjusted according to the detection depth to ensure effective coverage of the target goaf; the transmitting frequency parameter is set to 25Hz, the transmitting current stability is controlled within ±1%, and the off time is ≤50us to avoid that the long off time affects the late signal collection.

4. The method according to claim 3, wherein, In step (1), the seismic method observation system is a three-dimensional seismic observation system, and the specific parameters are as follows: a 12-line 6-shot observation mode is used, and the shot point is located at the midpoint of the measurement line; the shot point spacing is 10m, the shot array spacing is 40m, the receiving line spacing is 40m, and the channel spacing is 10m; 64 channels of geophones are arranged on each receiving line, and the geophones are selected to be piezoelectric geophones with a main frequency of 30Hz to ensure sensitive reception of the seismic wave signal; the sampling rate is set to 1ms, and the recording length is ≥2s to ensure that the reflection wave signal corresponding to the target goaf can be collected; the underground CDP (common depth point) bin size is 5m*5m, and the full coverage number reaches 48 times to meet the resolution requirement of fine detection of the goaf.

5. The method according to claim 4, wherein, In step (1), the specific parameters of the gravity method observation system are as follows: the measurement network density is consistent with that of the transient electromagnetic method, which is set to 40m*10m; The high-precision gravimeter is used, the instrument resolution is less than or equal to 0.1 muGal, and the observation accuracy is less than or equal to 5 muGal; before observation, static zero drift correction is required, the correction time is greater than or equal to 24 hours; during observation, after completing the observation of 20 measuring points, the reference point needs to be returned for repeated observation, and the mean error of repeated observation is less than or equal to 3 muGal; during observation, the elevation of each measuring point is recorded synchronously, the elevation measurement is performed by using a total station, and the measurement accuracy is less than or equal to 0.1 m, which provides basic data for subsequent terrain correction of gravity data.

6. The method according to claim 5, wherein, The specific operation requirements of data collection in step (2) are as follows: Seismic data acquisition: the explosive excitation is used, the explosive amount is adjusted according to the lithology of the detection area, the explosive amount is controlled to be 0.5-1 kg in soft stratum and 1-2 kg in hard stratum, so as to ensure effective excitation of seismic waves and avoid excessive interference with the surrounding environment; the detector needs to be tightly coupled with the ground, the coupling coefficient is greater than or equal to 0.9, and the signal attenuation caused by poor contact between the detector and the ground is reduced; Transient electromagnetic data collection: each measuring point is observed continuously for 3 times, and the average value after removing abnormal values is taken as the final data of the measuring point; during observation, the signal-to-noise ratio of the signal is monitored in real time to ensure that the signal-to-noise ratio of the late signal (corresponding to the deep goaf) is greater than or equal to 5; if the signal-to-noise ratio does not meet the requirements, the observation parameters need to be adjusted or re-observed; Gravity data collection: 3 groups of data are observed for each measuring point, each group of observation time is greater than or equal to 60 s, and the average value is taken after removing abnormal data groups through data stability analysis; during observation, vehicles are strictly prohibited to pass within a range of 50 m around the measuring point, and people are prohibited to walk near the measuring point, so as to reduce the influence of vibration interference on the observation data.

7. The method according to claim 6, wherein, The specific process of separate inversion in step (3) is as follows: Gravity data separate inversion: first, the original gravity data is sequentially subjected to terrain correction, Bouguer correction and normal field correction to eliminate the influence of factors such as terrain undulation and crust density difference; then, a three-dimensional apparent density inversion algorithm is used, the resolution of the grid matches the measurement grid density, that is, the horizontal grid spacing is 40 m x 10 m, and the vertical grid spacing is adjusted according to the detection depth, which is 5 m for shallow layer and 20 m for deep layer; the initial model is set based on regional geological data, such as normal stratum density of 2.6 g / cm³ and preliminary estimation of goaf density of 2.0-2.2 g / cm³; the inversion iteration number is greater than or equal to 10 times, until the data fitting difference is less than or equal to 5%, and the three-dimensional density distribution model of the goaf is output; Transient electromagnetic data separate inversion: first, the original data is subjected to power frequency filtering (50 Hz) and trend surface filtering to eliminate industrial interference and regional trend interference caused by surface terrain; then, terrain correction and device correction are performed to modify the influence of the observation device and terrain undulation on the data; a one-dimensional Occam inversion algorithm is used, the number of layers is set to 15-20 layers, the layer thickness increases gradually from shallow to deep, which is 5 m for shallow layer and 20 m for deep layer; during inversion, the least square criterion is used to fit the observation data, and the vertical distribution of resistivity at each measuring point and the plane resistivity distribution of the entire detection area are output. Seismic data inversion alone: The original seismic data is pre-processed, including static correction (eliminate the influence of surface elevation differences), denoising (suppress random noise, surface wave interference), amplitude compensation (compensate for the energy attenuation of seismic waves during propagation); then dynamic correction, horizontal stacking processing is carried out to obtain the stacked section; based on the stacked section, the reflection events of the target layer (such as No. 6 coal seam) are picked up, and the post-stack wave velocity inversion algorithm is used, the inversion interval is consistent with the CDP bin, which is 5m x 5m; the three-dimensional wave velocity distribution model of the whole detection area is output, wherein the normal coal seam wave velocity is set to 2800-3000m / s, and the wave velocity of the goaf is preliminarily estimated to be 2200-2500m / s.

8. The method according to claim 7, wherein, The step 4① builds a re-inversion objective function of cross gradient with the addition of a normalization operator: ； Constraints: ; where d i represents the observed data; m i represents the model parameters; m i _ref represents the reference model; f i (m i ) represents the forward response of the model parameters m i , a i represents the regularization factor; C di represents the data covariance matrix of the observed data d i ; C mi represents the model covariance matrix of the observed data m i ; ▽ represents the gradient, and τ represents the cross-gradient constraint term.

9. The method according to claim 8, wherein, The step 4② builds a cross gradient term of the preliminary inversion results of gravity data, seismic data and transient electromagnetic data: The gradient vector of two physical property parameters multiplied can reflect the structural similarity degree of two models, and the cross gradient function of two physical property parameters is defined as follows: ; Through different physical property inversion results, the cross gradient of the structure consistent area is close to 0, and the cross gradient amplitude of the structure inconsistent area is large, the cross gradient term is added to the objective function to constrain the inversion results to be consistent with the structure, and then the consistent physical property distribution is output.

Images