Method for detecting underground cavity by using irregular three-dimensional array

Abstract

A method for detecting an underground cavity by using an irregular three-dimensional array belongs to the technical field of underground detection, and is characterized in that a plurality of high-sensitivity digital seismographs are utilized, a detection area is well covered with a ray path, and each digital seismograph collects irregular micro-amplitude vibration information of the surface of the detection area all the time; a surface wave signal is extracted from micro-amplitude vibration information, then a cross-correlation function between digital seismometers is obtained through background noise cross-correlation calculation and superposition, and a Rayleigh wave phase velocity frequency dispersion curve is obtained based on a multiple filtering and extended space autocorrelation (ESPAC) method. And finally, inverting a three-dimensional phase velocity structure below the detection area based on the frequency dispersion curve of each grid node, and performing two-dimensional or three-dimensional mapping display on an inversion calculation result through mapping software. The method is simple in field work operation, flexible in laying according to field conditions, simple and rapid in data processing flow and mature in algorithm, and can rapidly and accurately position the underground cavity position.

Description

Technical Field

[0001] This invention belongs to the field of underground detection technology, specifically relating to a method for detecting underground cavities using an irregular three-dimensional array. Background Technology

[0002] Traditional methods for detecting underground cavities include high-density electrical resistivity tomography, seismic refraction, seismic reflection, surface wave, ground-penetrating radar, transient electromagnetic methods, cross-hole CT, and microgravity methods.

[0003] However, the size, depth, and filling material (air / water / silt) of cavities vary, and each detection method has its own limitations. For example, some methods are greatly affected by the site topography, grounding conditions, and surface uniformity; some methods have shallow detection depth and low resolution; some methods are not sensitive to small cavities or cavities filled with air; and some methods require drilling to assist in testing. It is difficult to accurately locate underground cavities using a single method and achieve rapid and accurate detection. In practice, multiple geophysical methods are often combined to accurately determine the location of underground cavities.

[0004] Therefore, there is a need to design a geophysical method that can quickly and accurately detect the location and size of underground cavities. Summary of the Invention

[0005] The purpose of this invention is to solve the problems existing in the background technology and to provide a method for detecting underground cavities using an irregular three-dimensional array.

[0006] This method utilizes several highly sensitive digital seismographs that provide good ray path coverage of the detection area. Each digital seismograph collects information on the irregular micro-amplitude vibrations occurring on the surface of the detection area every moment. Surface wave signals are extracted from the micro-amplitude vibration information. Then, the cross-correlation function between the digital seismographs is obtained through background noise cross-correlation calculation and superposition. The Rayleigh wave phase velocity dispersion curve is obtained using a multi-filter and extended spatial autocorrelation (ESPAC) method. Finally, the three-dimensional phase velocity structure below the detection area is inverted based on the dispersion curve of each grid node. The inversion calculation results are then displayed in two-dimensional or three-dimensional using mapping software.

[0007] The invention is simple in principle and easy to set up. It can accurately and quickly detect the location and size of underground cavities, which is conducive to its widespread application.

[0008] A method for detecting underground cavities using an irregular three-dimensional array is presented. The array consists of several digital seismographs, referred to as a three-dimensional array. This method involves arranging several irregularly shaped digital seismographs to record continuous seismic background noise waveform data. Cross-correlation and superposition calculations of the background noise are used to obtain the cross-correlation function between pairs of digital seismographs, fully utilizing the rich information from multiple digital seismographs in the array and maximizing the imaging potential of the array data. Rayleigh wave phase velocity dispersion curves in the frequency range of 0.2–200 Hz are extracted using a multiple filtering method. A low-squares iterative linear method is then used to invert and obtain a high-resolution three-dimensional phase velocity structure of the underground cavity. This allows for joint direct imaging of the three-dimensional phase velocity structure, avoiding the problem of stitching together images using one-dimensional models. This significantly improves the accuracy and reliability of underground structure detection, enhances imaging effects, and improves the scientific rigor and rationality of geological interpretation. This invention is not limited to a regular arrangement of measuring points; digital seismographs can be arranged more flexibly and conveniently according to site conditions, reducing the limitations and influences of site conditions.

[0009] A method for detecting underground cavities using an irregular three-dimensional array, the specific steps of which are as follows: Step 1: Preparation Stage Prepare several digital seismographs, each equipped with a three-component seismic sensor, a high-sensitivity BeiDou + GPS module, an electronic compass, an attitude sensor, and a rechargeable lithium battery. Place all the digital seismographs together, turn them on, and record seismic signals for a period of time. Compare the consistency of the digital seismographs to ensure that they are working properly and that the data acquisition is accurate.

[0010] Step Two, Data Collection Phase: Based on the planned exploration depth, design the boundary positions of the acquisition system. The distance between digital seismographs should be no less than twice the exploration depth. Deploy several densely packed observation arrays across the entire survey area as simultaneously as possible. If a limited number of digital seismographs prevent a single observation of the entire area, observations can be conducted in different phases. Adjust the positions of the digital seismographs for different phases to ensure the arrays are evenly distributed across the entire survey area. Arrange the positions of the digital seismographs according to site conditions, ensuring their azimuth and level are properly adjusted and they can receive good BeiDou and GPS signals. Accurate BeiDou + GPS timing is essential for the digital seismographs to function properly and ensure simultaneous data acquisition. Record the coordinates (X, Y, Z) of each digital seismograph using measuring instruments and calculate the distances between different digital seismographs.

[0011] Step 3: Data Processing Stage like Figure 2 As shown, the data processing flow of the three-dimensional array includes: 1. Data preprocessing: The data from each digital seismograph were resampled, de-instrument response was removed, mean was removed, and bandpass filtering was performed to eliminate bad data. The data were then normalized.

[0012] 2. Spatial autocorrelation waveform calculation and post-processing: After completing the preprocessing of a single digital seismograph, the noise cross-correlation function between different digital seismographs is calculated and superimposed in the time domain to improve the signal-to-noise ratio.

[0013] 3. Surface wave dispersion curve extraction: The fundamental group velocity and phase velocity dispersion curves are manually extracted from the cross-correlation function using a multiple filtering method. The group velocity dispersion curve is extracted based on the principle of maximum energy, while the phase velocity dispersion curve is extracted using a reference dispersion curve and in combination with the principle that the phase velocity value is greater than the group velocity value under the same period condition.

[0014] 4. Discrete frequency point fitting and verification: By comparing the phase velocity dispersion of Rayleigh waves and Love waves, the average Rayleigh wave dispersion is used as the reference dispersion. In order to ensure the reliability of the dispersion curve, the dispersion curve is manually screened to remove obvious errors or non-smooth dispersion curves. For the same period, dispersion data beyond twice the mean error are deleted.

[0015] 5. Structural Modeling: Based on the actual layout of the field array, different frequency intervals, imaging grids, and imaging control parameters (damping factor, smoothing factor, etc.) can be established; the initial model is a uniform model, and the velocity value is given by averaging the selected dispersion data.

[0016] 6. Inversion Calculation: Rayleigh wave phase velocity tomography is performed using the generalized least squares inversion method of continuous functions. Using this inversion method requires selecting appropriate correlation length L and prior model error σ. ｍ0 The correlation length L is equivalent to the spatial smoothing factor, corresponding to a spatial resolution of 2L. The choice of L needs to consider ray coverage and wavelength. From an optical point of view, even if the rays are very dense, it is impossible to accurately invert anomalies smaller than one wavelength. Therefore, the minimum value of L is half a wavelength (0.5λ). To prevent small-scale spurious anomalies from appearing in short periods, we set L = max [50; 0.5λ]; the prior model error σ ｍ0To control the model's resolution, choosing a larger prior model error can improve the resolution of the inverted model, but the posterior model error also increases accordingly. Therefore, a trade-off needs to be struck between model resolution and error. Considering that the phase velocity dispersion measurement error is typically 3% to 10%, the acceptable posterior model error for phase velocity tomography is 1%. The residual perturbation between the travel time calculated by forward modeling and the observed travel time is added to the model. When the perturbation tends to zero and stabilizes, the inversion can be stopped. Imaging parameters can be adjusted and the imaging results optimized based on Rayleigh wave data fitting and imaging results.

[0017] 7. Calculation results are plotted and displayed: Based on the inversion calculation of the dispersion curve at each grid node, the Rayleigh wave phase velocity structure at different burial depths below all grid nodes is obtained. The inversion calculation results are plotted and displayed in two or three dimensions using mapping software.

[0018] The drawing software mentioned is CAD, EVS, or Voxler4, etc.

[0019] Working principle and usage of this invention: This method is based on two assumptions: the theory of spatially stationary random fluctuations and the fact that the fundamental surface wave is the main energy in the micro-motion signal.

[0020] In a regular array observation system, let the vertical component wavefield at the center be... The vertical component wavefield of any digital seismograph on the circumference is, Then the spatial autocorrelation function It can be represented as: , In the formula, x、y Represents spatial location coordinates, t For the time of dissemination, and It is the change in spatial location. The overline represents the time average. Assuming the wave field propagates within the circle at a single velocity c, the azimuth average of the spatial autocorrelation function is: , The above formula is in polar coordinate form. , is the azimuth average value of the spatial autocorrelation function. The spatial autocorrelation function is in polar coordinate form. The distance between any digital seismograph and its center point. For radius The azimuth angle between the circular digital seismograph and its center point. Represent the independent variable The differential.

[0021] The relationship between wavefield power spectral density and the azimuth average of spatial autocorrelation function can be expressed as:

[0022] In the formula, Let be the power spectral density function of the wave field. The spatial autocorrelation function represents the azimuth average value. oh Angular frequency, It is a zeroth-order Bessel function of the first kind. The distance between any digital seismograph and its center point. dr Represent the independent variable The differential, dω Represent the independent variable oh The differential, c Given a single velocity in the wave field, applying narrowband filtering to the above signal, the power spectral density can be expressed as: ,

[0023] In the formula, P ( oh 0) represents the frequency at which oh Power spectral density at 0 oh Angular frequency, Let be the Dirac function, Substituting into formulas (5) and (4), the azimuth average of the spatial autocorrelation function can be expressed as: ,

[0024] In the formula, for At that time, the azimuth average value of the spatial autocorrelation function, It is a zeroth-order Bessel function of the first kind.

[0025] The autocorrelation coefficient is defined as: , In the formula, Indicates radius as Azimuth angle of the circular digital seismograph relative to the center point The azimuth average of the spatial autocorrelation function. Indicates the azimuth angle of the center point digital seismograph The spatial autocorrelation function of the location is the azimuth average.

[0026] Since the power spectral density is independent of the digital seismograph's location, the azimuth-averaged autocorrelation coefficient can be written as: , In the formula, The distance between any digital seismograph and its center point. oh 0 is the angular frequency, i.e. oh 0=2 πf , c ( oh 0) represents the phase velocity of the Rayleigh wave. It is a zeroth-order Bessel function of the first kind.

[0027] Studies have shown that long-term earthquake time series can replace spatial azimuth averaging of digital seismographs; that is, for two digital seismographs, by extending the observation and recording time, the noise source can be made to tend towards a spatially uniform distribution; for long-term micromotion records between digital seismographs, the autocorrelation coefficient is calculated according to the following formula: , In the formula, It is the autocorrelation coefficient, which is the frequency. f The function; U i ( f ) and U j ( f () indicates a digital seismograph i With digital seismograph j Frequency domain micro-motion data, * denotes conjugate. oh 0 is the angular frequency, i.e. oh 0=2 πf , For digital seismographs i The distance between the digital seismograph j and the digital seismograph j. c ( oh 0) represents the phase velocity of the Rayleigh wave.

[0028] When a digital seismograph generates a sufficiently large or isotropic wavefield on a circular path, or when digital seismographs record micro-motion signals for a sufficiently long time, the imaginary part of the autocorrelation coefficient spectrum is zero, and the real part is a zero-order Bessel function of the first kind. In actual calculations, the imaginary part is generally not zero, indicating that the spatial distribution of the noise source is not absolutely uniform, or that the seismic noise wavefield is not a stable plane wave. This introduces a certain error into the calculation of the autocorrelation coefficient spectrum. Generally, the real part of the autocorrelation coefficient spectrum is approximated as a zero-order Bessel function of the first kind, i.e. , In the formula, Re This means taking the real part of a complex number. f For frequency, U i ( f ) and U j ( f () indicates a digital seismograph i With digital seismograph jFrequency domain micro-motion data, * denotes conjugate. For digital seismographs i With digital seismograph j The distance between, c ( f () represents the phase velocity of the Rayleigh wave. It is a zeroth-order Bessel function of the first kind.

[0029] Conventional SPAC methods typically deploy the observation system in concentric circles to reduce the impact of uneven noise source azimuth distribution on the reliability of Rayleigh wave dispersion curves. To enable spatial autocorrelation methods to be used over larger areas, the Extended Spatial Autocorrelation Method (ESPAC) with irregular arrays is employed. The ESPAC method offers greater flexibility in deploying the observation system, allowing for linear, cross-shaped, circular, or mesh-like configurations. Irregular arrangements are simpler and more flexible, depending on site conditions. Rayleigh wave dispersion information is obtained through data processing between two digital seismographs. Rayleigh wave phase velocity dispersion curves are acquired using a multiple filtering method. Rayleigh wave phase velocity tomography is performed using a generalized least squares inversion method based on continuous functions to obtain the Rayleigh wave phase velocity structure at different burial depths below all grid nodes.

[0030] The beneficial effects of this invention are: This invention is simple to operate in the field, can be arranged in irregular arrays, and can be flexibly deployed according to site conditions. The data processing flow is streamlined and fast, and the algorithm is mature, which can quickly and accurately locate the underground cavity. By changing the spacing density and size of the digital seismographs, the exploration depth and accuracy can be controlled, overcoming the drawbacks of traditional detection methods. It is not affected by on-site grounding conditions, electromagnetic interference, etc., and effectively solves the problem of underground cavity detection. Attached Figure Description

[0031] Figure 1 This is a layout diagram of the digital seismograph used in this invention; Figure 2 This is a flowchart of the data processing process of the present invention. Detailed Implementation

[0032] See Figure 1As shown, an irregular three-dimensional array method for detecting underground cavities involves deploying digital seismographs (DSS) according to site conditions, ensuring good ray path coverage of the detection area, and using a geological compass to adjust the DSS's azimuth at locations with good BeiDou and GPS signal reception. The DSS is then leveled using a tail cone or foot screw. The process includes powering on, positioning, timing, and data acquisition, recording the time period for simultaneous acquisition. Depending on the exploration depth, the acquisition time period can be selected between 30 minutes and 24 hours. After completing the acquisition, the DSS is turned off, moved to the next location, and the above operation is repeated. The coordinates (X, Y, Z) of all deployed DSS locations are recorded using a measuring instrument.

[0033] The data processing workflow of the three-dimensional array includes preprocessing, Rayleigh wave dispersion curves, and Rayleigh wave phase velocity structure inversion imaging.

[0034] In the data preprocessing of a single digital seismograph, the vertical components of all digital seismographs are first resampled, instrument response is removed, mean and linear trend are removed, and a bandpass filter within the range of 0.2~200Hz is set. To suppress the influence of interference signals and distortion signals caused by instrument failure, a moving absolute average method is used for normalization, and its calculation formula is as follows: , In the formula, d n These are the original data points. It is a normalized time series. oh n This is the weight of that point, that is, the weight of the center point of the time window. It is obtained by calculating the average of the absolute values ​​of the amplitude within a certain time window, and its formula is expressed as: , In the formula, DJ It is the first j Waveform data at 2 time points, 2 N The next step is to perform spectral whitening on the normalized data, with the aim of widening the frequency band of the background noise signal, suppressing single-frequency signal interference, and thus obtaining a more continuous dispersion curve.

[0035] The vertical component micromotion data of each digital seismograph is segmented according to the simultaneous acquisition time. Each segment of data is subjected to Fourier transform. The autocorrelation coefficient is calculated by the conjugate of the Fourier spectrum of one digital seismograph and the spectrum of another digital seismograph according to formula (9). The real part is then taken according to formula (10) and superimposed. Then, the superimposed autocorrelation coefficient waveform is smoothed by the moving absolute average method to obtain the autocorrelation coefficient curve of each digital seismograph pair. The autocorrelation coefficient curves with low signal-to-noise ratio and those that fail to show the shape of the first type zero-order Bezier curve are removed. It can be seen from the autocorrelation coefficient waveforms of different digital seismograph distances that the autocorrelation coefficient curves are similar to the first type zero-order Bezier curves. As the distance between digital seismographs increases, the frequency corresponding to the first zero point shows a decreasing trend.

[0036] By analyzing the zeros and poles of the autocorrelation coefficient waveform and fitting it with the zeros and poles of the first-order zero Bessel function, the Rayleigh wave phase velocity dispersion curve was obtained, i.e., the zero-order Bessel function of the first order was known. The values ​​corresponding to the zeros and poles (x is the independent variable, y is the dependent variable) are obtained from the correlation coefficient curve. The zeros and poles can be used to obtain the corresponding frequencies. f From the formula The Rayleigh wave phase velocity dispersion curve for each digital seismograph can then be calculated. c ( f After quality control and selection, high-quality Rayleigh wave phase velocity dispersion curves were obtained. The Rayleigh wave tomography method was used to perform tomographic inversion on the phase velocity dispersion data from 0.2 to 200 Hz. After multiple imaging tests, appropriate regularization parameters and damping coefficients were determined to ensure that the imaging results were relatively smooth and could fit the observation data well.

[0037] The inversion calculation results can be displayed using 3D viewing software to show the 3D Rayleigh wave phase velocity structure below the survey area.

Claims

1. A method for detecting underground cavities using an irregular three-dimensional array, characterized in that: This method utilizes several highly sensitive digital seismometers that provide good ray path coverage of the detection area. Each digital seismometer collects information on the irregular micro-amplitude vibrations of the surface of the detection area at all times. Surface wave signals are extracted from the micro-amplitude vibration information. Then, the cross-correlation function between the digital seismometers is obtained through background noise cross-correlation calculation and superposition. Rayleigh wave phase velocity dispersion curves are obtained using a method based on multiple filtering and extended spatial autocorrelation. Finally, the three-dimensional phase velocity structure below the detection area is inverted based on the dispersion curves of each grid node. The inversion calculation results are then displayed in two-dimensional or three-dimensional using mapping software.

2. The method for detecting underground cavities using an irregular three-dimensional array according to claim 1, characterized in that: Includes the following steps: Step 1: Preparation Stage Prepare several digital seismographs, each equipped with a three-component seismic sensor, a high-sensitivity BeiDou + GPS module, an electronic compass, an attitude sensor, and a rechargeable lithium battery. Place all the digital seismographs together, turn them on, and record seismic signals for a period of time. Compare the consistency of the digital seismographs to ensure that they are working properly and that the data acquisition is accurate. Step Two, Data Collection Phase: Based on the planned exploration depth, design the boundary positions of the acquisition system. The distance between digital seismographs should be no less than twice the exploration depth. Deploy several densely packed observation arrays across the entire survey area at once. If a limited number of digital seismographs prevent a single observation of the entire area, observations can be conducted in different phases. Adjust the positions of the digital seismographs for different phases to ensure the arrays are evenly distributed across the survey area. Arrange the digital seismograph positions according to site conditions, ensuring their azimuth and level are properly adjusted and they can receive good BeiDou and GPS signals. Precise BeiDou + GPS timing is essential for the digital seismographs to function properly and ensure simultaneous data acquisition. Record the coordinates of each digital seismograph using measuring instruments and calculate the distances between different digital seismographs. Step 3: Data Processing Stage The data processing flow of a 3D array includes: (1) Data preprocessing: The data of each digital seismograph were resampled, de-instrument response was removed, mean was removed and bandpass filtering was performed to remove bad data and normalize the data. (2) Spatial autocorrelation waveform calculation and post-processing: After completing the preprocessing of a single digital seismograph, the noise cross-correlation function between different digital seismographs is calculated and superimposed in the time domain to improve the signal-to-noise ratio; (3) Surface wave dispersion curve extraction: The basic group velocity and phase velocity dispersion curves are manually extracted from the cross-correlation function using the multiple filtering method. The group velocity dispersion curve is extracted based on the principle of maximum energy, while the phase velocity dispersion curve is extracted with the help of the reference dispersion curve and in combination with the principle that the phase velocity value is greater than the group velocity value under the same period condition. (4) Fitting and verification of discrete frequency points: By comparing the phase velocity dispersion of Rayleigh wave and Love wave, the average Rayleigh wave dispersion is used as the reference dispersion; in order to ensure the reliability of the dispersion curve, the dispersion curve is manually screened to remove obvious errors or non-smooth dispersion curves; for the same period, dispersion data beyond twice the mean error are deleted. (5) Structural modeling: Based on the actual layout of the field array, different frequency intervals, imaging grids, and imaging control parameters can be established; the initial model is a uniform model, and the velocity value is given by averaging the selected dispersion data. (6) Inversion calculation: Rayleigh wave phase velocity tomography is performed using the generalized least squares inversion method of continuous functions. When using this inversion method, it is necessary to select an appropriate correlation length L and prior model error σ. ｍ0 The correlation length L is equivalent to the spatial smoothing factor, corresponding to a spatial resolution of 2L. The choice of L needs to consider ray coverage and wavelength. From an optical point of view, the minimum value of L is half a wavelength, i.e., 0.5λ. To prevent small-scale spurious anomalies from appearing in short periods, we set L = max [50; 0.5λ]; the prior model error σ ｍ0 To improve the resolution of the inverted model, a larger prior model error can be chosen, but the posterior model error also increases accordingly. Therefore, a trade-off needs to be struck between model resolution and error. Considering that the phase velocity dispersion measurement error is typically 3% to 10%, the acceptable posterior model error for phase velocity tomography is 1%. The residual perturbation between the travel time calculated from forward modeling and the observed travel time is added to the model. When the perturbation tends to zero and stabilizes, the inversion can be stopped. Imaging parameters can be adjusted and the imaging results optimized based on Rayleigh wave data fitting and imaging results. (7) Calculation results are plotted and displayed: Based on the inversion calculation of the dispersion curve of each grid node, the Rayleigh wave phase velocity structure at different burial depths below all grid nodes is obtained, and the inversion calculation results are plotted and displayed in two or three dimensions using plotting software.

3. The method for detecting underground cavities using an irregular three-dimensional array according to claim 2, characterized in that: The drawing software mentioned is CAD, EVS, or Voxler4.

Images