Method for detecting three-dimensional density structure of stratum

Abstract

The invention relates to a method for detecting a stratum three-dimensional density structure, and the method comprises the steps: obtaining seismic background noise data of a stratum to be detected, extracting Rayleigh surface wave frequency dispersion data based on the seismic background noise data, calculating Rayleigh surface wave travel time, and determining a first relation between the Rayleigh surface wave travel time and the stratum density; the method comprises the following steps: obtaining muon data of a to-be-detected stratum, calculating the density length of the to-be-detected stratum based on the muon data, and determining a second relationship between the density length of the to-be-detected stratum and the stratum density; and determining the three-dimensional density structure of the stratum according to the first relationship and the second relationship. The three-dimensional density structure of the to-be-detected stratum is subjected to joint inversion based on the seismic surface wave data and the muon data, the internal structure of the shallow stratum can be detected through a lossless and non-invasive device and remote sensing without any interference on the structure of the shallow stratum, and the method has the advantages of being high in precision, high in resolution ratio, large in observable depth and the like.

Description

Technical Field

[0001] The present application relates to the field of geophysical exploration technology, and in particular to a method for detecting the three-dimensional density structure of a stratum. Background Art

[0002] Geological exploration is the exploration and detection of geology through various technical means and methods. There are many commonly used geophysical exploration methods. Exploration methods such as direct current exploration, alternating current exploration, ultrasonic method, infrared temperature recording method, radar method and X-ray method have always been widely known to mankind. However, these methods significantly limit the inspection objects, and there are limits on the resolution of internal states and the depth of the surface layer that can be detected. For example, direct current and alternating current exploration are easily interfered by low resistivity mineralized layers and steel bars, resulting in inaccurate detection results. In nuclear radiation imaging detection technology, X-ray imaging requires that light particles must have sufficient energy. Muon imaging observation technology has the disadvantages of long observation time and shallow observation depth. Therefore, a new method for surveying deep structures of strata is needed. Summary of the invention

[0003] In order to overcome the problems existing in the related art, the present application provides a method for detecting the three-dimensional density structure of a stratum.

[0004] According to a first aspect of an embodiment of the present application, a method for detecting a three-dimensional density structure of a formation is provided, comprising: acquiring seismic background noise data and muon data of a formation to be detected; extracting Rayleigh surface wave dispersion data based on the seismic background noise data, and calculating Rayleigh surface wave travel time; determining a first relationship between the Rayleigh surface wave travel time and formation density; calculating a density length of the formation to be detected based on the muon data; determining a second relationship between the density length of the formation to be detected and formation density; and determining the three-dimensional density structure of the formation based on the first relationship and the second relationship.

[0005] In some exemplary embodiments of the present application, extracting Rayleigh surface wave dispersion data based on the seismic background noise data and calculating the Rayleigh surface wave travel time include: extracting Rayleigh surface wave signals from the seismic background noise data based on a waveform cross-correlation method; acquiring the Rayleigh surface wave dispersion data using a multiple filtering method; and calculating the Rayleigh surface wave travel time based on the Rayleigh surface wave dispersion data.

[0006] In some exemplary embodiments of the present application, determining the first relationship between the Rayleigh surface wave travel time and the formation density includes: establishing a linear equation group of the Rayleigh surface wave travel time and the formation density:

[0007] (w s A)δρ=(w s Δt)

[0008] Where A represents the partial derivative matrix of the Rayleigh surface wave travel time for the formation density anomaly, δρ represents the density disturbance value of the three-dimensional grid of the formation density structure to be detected, Δt represents the Rayleigh surface wave travel time residual, w s Represents the normalized weight factor for seismic data.

[0009] In some exemplary embodiments of the present application, the muon data includes first muon data and second muon data; the first muon data includes: the distribution of the number of muons on the surface of the formation to be detected with respect to azimuth; and the second muon data includes: the distribution of the number of muons at each muon sampling point under the formation to be detected with respect to azimuth.

[0010] In some exemplary embodiments of the present application, the density length of the formation to be detected is calculated based on the muon data, including: calculating the muon attenuation of each muon sampling point based on the first muon data and the second muon data; calculating the minimum energy required for the muon to penetrate the formation to be detected based on the relationship between the muon flux and the muon attenuation of each muon sampling point; and calculating the density length of the formation to be detected passed by the muon based on the relationship between the minimum energy required for the muon to penetrate the formation to be detected and the density length of the formation to be detected passed by the muon.

[0011] In some exemplary embodiments of the present application, determining the second relationship between the density length of the to-be-detected formation and the formation density includes: establishing a linear equation group of the density length of the to-be-detected formation and the formation density:

[0012] (w μ L)δρ=(w μ ΔDL)

[0013] Wherein, L represents the partial derivative matrix of the density length of the formation to be detected that the muon passes through for the formation density anomaly, δρ represents the density perturbation value of the three-dimensional grid of the density structure of the formation to be detected, ΔDL represents the residual of the density length of the formation to be detected that the muon passes through, and w μ represents the normalized weight factor of the muon data.

[0014] In some exemplary embodiments of the present application, determining the three-dimensional density structure of the formation according to the first relationship and the second relationship includes: combining the first relationship and the second relationship to obtain a linear equation group for inverting the three-dimensional density structure of the formation based on the Rayleigh surface wave travel time and the density length of the formation to be detected:

[0015]

[0016] Wherein, δρ represents the density disturbance value of the three-dimensional grid of the density structure of the formation to be detected, L represents the partial derivative matrix of the density length of the formation to be detected through which the muon passes for the formation density anomaly, A represents the partial derivative matrix of the Rayleigh surface wave travel time for the formation density anomaly, ΔDL represents the residual of the density length of the formation to be detected through which the muon passes, Δt represents the residual of the Rayleigh surface wave travel time, w μ represents the normalized weight factor of the muon data, w s represents the normalized weight factor of seismic data, S represents the regularized smoothing constraint, I represents the positive definite diagonal matrix, σ represents the normalization coefficient, and γ represents the damping factor;

[0017] Determining the density disturbance value of the three-dimensional grid of the stratum density structure to be detected;

[0018] The three-dimensional density structure of the formation is determined according to the density disturbance value and a preset one-dimensional density model of the formation.

[0019] In some exemplary embodiments of the present application, the obtaining of seismic background noise data of the formation to be detected includes: arranging a dense seismic array in a first selected area of ​​the formation to be detected; and obtaining the seismic background noise data observed by the dense seismic array.

[0020] In some exemplary embodiments of the present application, the obtaining of muon data of the formation to be detected includes: arranging a muon observation system in a second selected area of ​​the formation to be detected; and obtaining the muon data observed by the muon observation system.

[0021] In some exemplary embodiments of the present application, the arranging a muon observation system in the second selected area of ​​the formation to be detected includes: arranging muon imagers at equal intervals below the formation to be detected, wherein the interval Gap is:

[0022] Gap=(n+depth)*tanβ

[0023] Wherein, n represents the distance between the muon imager and the lower surface of the formation to be detected, depth represents the thickness of the formation to be detected, and β represents the detectable angle range of the muon imager.

[0024] The technical solution provided by the embodiments of the present application may have the following beneficial effects:

[0025] The present application provides a method for detecting the three-dimensional density structure of a stratum, which is based on the establishment of the stratum wave velocity structure by seismic surface waves and the establishment of the density distribution by muon imaging for joint inversion and solution, and obtains the accurate three-dimensional density structure of the stratum to be detected, which effectively improves the resolution and reliability of density structure imaging. This method combines the advantages of joint imaging of seismic surface wave data and muon data, breaks free from the constraints of traditional physical exploration methods, and gets rid of the limitations of traditional single detection methods. It can detect the internal structure of shallow strata by non-destructive, non-invasive devices and remote sensing without any interference to the shallow stratum structure, and has the characteristics of high accuracy, high resolution and large observable depth.

[0026] It should be understood that the foregoing general description and the following detailed description are exemplary and explanatory only and are not restrictive of the present application. BRIEF DESCRIPTION OF THE DRAWINGS

[0027] The accompanying drawings, which are incorporated in and constitute a part of this specification, illustrate embodiments consistent with the invention and, together with the description, serve to explain the principles of the invention.

[0028] Figure 1 The present invention is a flow chart showing a method for detecting three-dimensional density structure of a stratum according to an exemplary embodiment.

[0029] Figure 2 It is a schematic diagram showing the arrangement of instruments in a dense seismic array and muon observation system according to an exemplary embodiment.

[0030] Figure 3 The present invention is a flow chart showing a method for detecting three-dimensional density structure of a stratum according to an exemplary embodiment.

[0031] Figure 4 The present invention is a flow chart showing a method for detecting three-dimensional density structure of a stratum according to an exemplary embodiment.

[0032] Figure 5 It is a schematic diagram showing the initial response relationship between the density length of the formation to be detected and the muon attenuation at several zenith angles according to an exemplary embodiment.

[0033] Figure 6 The present invention is a flow chart showing a method for detecting three-dimensional density structure of a stratum according to an exemplary embodiment. DETAILED DESCRIPTION

[0034] Exemplary embodiments will be described in detail herein, examples of which are shown in the accompanying drawings. When the following description refers to the drawings, the same numbers in different drawings represent the same or similar elements unless otherwise indicated. The embodiments described in the following exemplary embodiments do not represent all embodiments consistent with the present invention. Instead, they are merely examples of devices and methods consistent with some aspects of the present invention as detailed in the appended claims.

[0035] Geological exploration is the exploration and detection of geology through various technical means and methods. There are many commonly used geophysical exploration methods. Exploration methods such as direct current exploration, alternating current exploration, ultrasonic method, infrared temperature recording method, radar method and X-ray method have always been widely known to mankind. However, these methods significantly limit the inspection objects, and there are limits on the resolution of internal states and the depth of the surface layer that can be detected. For example, direct current and alternating current exploration are easily interfered by low resistivity mineralized layers and steel bars, resulting in inaccurate detection results. In nuclear radiation imaging detection technology, X-ray imaging requires that light particles must have sufficient energy. Muon imaging observation technology has the disadvantages of long observation time and shallow observation depth.

[0036] In order to solve the above problems, the present application provides a method for detecting the three-dimensional density structure of a stratum, by obtaining the seismic background noise data of the stratum to be detected, extracting the Rayleigh surface wave dispersion data based on the seismic background noise data, and calculating the Rayleigh surface wave travel time, and determining the first relationship between the Rayleigh surface wave travel time and the stratum density; by obtaining the muon data of the stratum to be detected, calculating the density length of the stratum to be detected based on the muon data, and determining the second relationship between the density length of the stratum to be detected and the stratum density; then determining the three-dimensional density structure of the stratum according to the first relationship and the second relationship. The present application provides a method for detecting the three-dimensional density structure of a stratum, which establishes the stratum wave velocity structure based on seismic surface waves and the density distribution based on muon imaging for joint inversion solution, and obtains an accurate three-dimensional density structure of the stratum to be detected, which effectively improves the resolution and reliability of density structure imaging. This method combines the advantages of joint imaging of seismic surface wave data and muon data, breaks free from the constraints of traditional geophysical exploration methods, and gets rid of the limitations of traditional single detection methods. It can detect the internal structure of shallow strata using non-destructive, non-invasive devices and remote sensing without any interference to the shallow strata structure. It has the characteristics of high accuracy, high resolution and large observable depth.

[0037] The exemplary embodiment of the present application provides a method for detecting the three-dimensional density structure of a stratum, such as Figure 1 As shown, the method for detecting the three-dimensional density structure of a stratum shown in this embodiment includes:

[0038] S100, obtaining seismic background noise data and muon data of a formation to be detected.

[0039] S200, extracting Rayleigh surface wave dispersion data based on seismic background noise data, and calculating Rayleigh surface wave travel time.

[0040] S300. Determine a first relationship between Rayleigh surface wave travel time and formation density.

[0041] S400, calculating the density length of the formation to be detected based on the muon data.

[0042] S500: Determine a second relationship between the density length of the formation to be detected and the formation density.

[0043] S600: Determine a three-dimensional density structure of the stratum according to the first relationship and the second relationship.

[0044] The method for detecting the three-dimensional density structure of the stratum in this embodiment breaks free from the constraints of traditional physical exploration methods and the limitations of traditional single detection methods by jointly inverting the three-dimensional density structure of the stratum to be detected based on seismic surface wave data and muon data. It can detect the internal structure of the shallow stratum by using non-destructive, non-invasive devices and remote sensing without any interference to the shallow stratum structure, and has the characteristics of high accuracy, high resolution and large observable depth.

[0045] In step S100, seismic background noise data of the to-be-detected stratum is obtained, including: arranging a dense seismic array in a first selected area of ​​the to-be-detected stratum, and obtaining seismic background noise data observed by the dense seismic array.

[0046] A dense seismic array is a seismic observation system used for seismic observation and research, which is formed by arranging seismic detection equipment in a first selected area according to certain rules and density. In an embodiment of the present application, the first selected area can be selected according to actual detection needs, and the seismic detection equipment in the dense seismic array can be arranged and adjusted according to actual detection needs. For example, a group of seismic detection equipment can be arranged at equal intervals in an area on the surface of the stratum to be detected according to the positioning of the positioning system, and the spacing can be 100 meters, 150 meters, etc., and the present application does not make specific restrictions.

[0047] Earthquake detection equipment can be equipment that records ground vibrations, such as a seismograph. Figure 2 The schematic diagram of the arrangement of seismographs in a dense seismic array is shown in FIG. The seismic detection equipment mainly includes: seismic detectors, power supply systems, data transmission systems, etc. The positioning system mainly includes: a total station type electronic rangefinder, a prism matched with the total station type electronic rangefinder, etc.

[0048] By arranging dense seismic arrays, the monitoring capability of seismic waves can be improved and the propagation characteristics of seismic waves can be studied.

[0049] Continuous seismic background noise data can be observed in the dense seismic array, which includes seismic wave signals and noise signals. Seismic wave signals include body waves and surface waves. Body waves include primary waves and secondary waves, and surface waves include Love waves and Rayleigh waves.

[0050] like Figure 3 As shown, in step S200, Rayleigh surface wave dispersion data is extracted based on seismic background noise data, and Rayleigh surface wave travel time is calculated, including:

[0051] S310. Extracting Rayleigh surface wave signals from seismic background noise data based on waveform cross-correlation method.

[0052] S320. Use a multiple filtering method to obtain Rayleigh surface wave dispersion data.

[0053] S330. Calculate the Rayleigh surface wave travel time based on the Rayleigh surface wave dispersion data.

[0054] In this embodiment, the Rayleigh surface wave signal in the seismic surface wave data is extracted by using the waveform cross-correlation method. In view of the small aperture of the dense array, generally several meters to hundreds of meters, waveform segmentation with different time lengths of 2 minutes to 30 minutes is adopted, and the linear superposition, phase weighted superposition and other methods are compared and analyzed to improve the surface wave signal-to-noise ratio and obtain accurate Rayleigh surface wave signals.

[0055] A multiple filtering method is used for the Rayleigh surface wave signal to extract the Rayleigh surface wave dispersion data. The Rayleigh surface wave dispersion data can be a Rayleigh surface wave dispersion curve.

[0056] Seismic background noise data is correlated between different seismic detection devices, so Rayleigh surface wave signals can be extracted by performing cross-correlation calculations on seismic background noise records received by different seismic detection devices. The multiple filtering method is a signal processing technology based on the frequency domain. By weighting and screening different frequency components, the frequency components of Rayleigh surface wave signals are highlighted, thereby effectively extracting Rayleigh surface wave dispersion information.

[0057] In step S330, the ray path and travel time of the Rayleigh surface wave from the source to the receiving point can be calculated based on the surface wave one-step imaging method of the fast marching method ray tracing. The travel time is the time from when the Rayleigh surface wave is emitted from the source to when it is received by the seismic detection equipment. For the Rayleigh surface wave travel time t with a frequency of ω, it can be expressed as:

[0058] t(ω)=∫ l S(l,ω)dl (1)

[0059] In formula (1), l represents the propagation path of the Rayleigh surface wave, and S(l,ω) represents the slowness of the Rayleigh surface wave with propagation path l and frequency ω. The surface wave slowness is the inverse of the surface wave propagation velocity.

[0060] In step S300, determining a first relationship between Rayleigh surface wave travel time and formation density includes:

[0061] The linear equations of Rayleigh surface wave travel time and formation density are established:

[0062] (w s A)δρ=(w s Δt) (2)

[0063] In formula (2), A represents the partial derivative matrix of Rayleigh surface wave travel time for formation density anomaly, δρ represents the density disturbance value of the three-dimensional grid of the formation density structure to be detected, Δt represents the residual of Rayleigh surface wave travel time, and w s Represents the normalized weight factor for seismic data.

[0064] The following describes the steps to establish the linear equations of Rayleigh surface wave travel time and formation density.

[0065] Firstly, the formation wave velocity structure model is established and the formation to be tested is parameterized.

[0066] In this embodiment, an initial model of the formation wave velocity structure is constructed based on the known geological information of the formation to be detected, including the initial distribution of physical parameters of the formation structure to be detected, including the thickness of each layer, the longitudinal wave velocity, the shear wave velocity, etc. The ray path and travel time information of the Rayleigh surface wave obtained by fast marching method ray tracing are combined with the dispersion characteristics of the surface wave to perform one-step imaging of the surface wave, and the velocity structure distribution in the formation to be detected can be obtained.

[0067] The structure of the stratum to be detected is parameterized and discretized into several grid cells. Assuming that the medium properties in each grid cell are uniform, solving the density of the stratum to be detected is to solve the density value of each grid cell.

[0068] Secondly, using the existing appropriate empirical relationship between wave velocity and formation density, a linear equation group of Rayleigh surface wave dispersion data and the formation density to be detected is established.

[0069] Using the empirical velocity-density relationship in physical experiments of different rock types and depths, the functional relationship between the P-wave velocity Vp, S-wave velocity Vs and formation density is established. The polynomial quantitative relationship between the P-wave velocity Vp and the formation density at depth z can be expressed as the following formula (3), and the polynomial quantitative relationship between the S-wave velocity Vs and the formation density at depth z can be expressed as the following formula (4):

[0070]

[0071] In equations (3) and (4), α represents the P-wave velocity Vp, β represents the S-wave velocity Vs, ρ represents the formation density, and z represents the formation depth. The nth-order parameter representing the polynomial relationship between the P-wave velocity Vp and the formation density ρ, The nth-order parameter representing the polynomial relationship between shear wave velocity Vs and formation density ρ.

[0072] According to the empirical relationship between wave velocity and formation density (Formula (3) and Formula (4)) and the surface wave travel time information obtained by ray tracing, for the ith Rayleigh surface wave (i is a positive integer, i represents the Rayleigh surface wave number), the partial derivative relationship between the Rayleigh surface wave travel time and the wave velocity (P-wave velocity Vp and S-wave velocity Vs) can be converted into the partial derivative relationship between the Rayleigh surface wave travel time and the formation density to be detected. It can be obtained that the partial derivative relationship between the travel time data of the ith Rayleigh surface wave with a frequency of ω and the formation density structure to be detected is:

[0073]

[0074] In formula (5), δt i represents the seismic surface wave propagation time of the i-th seismic ray (the i-th seismic ray represents the i-th Rayleigh surface wave propagation path), K represents the number of grid nodes along the ray path, and J represents the number of grid nodes along the depth direction. K and J can be set according to actual needs, and the specific values ​​are not limited. ik represents the propagation speed of the i-th seismic ray in the k-th grid, C k (ω) represents the surface wave phase velocity of the surface wave with frequency ω in the kth grid, z j represents the jth grid node in the depth direction, α k (z j ), β k (z j ),ρ k (z j ) represent the P-wave velocity Vp, S-wave velocity Vs and formation density parameters at the j-th depth grid node and the k-th ray of the seismic surface wave ray path, respectively. A il The parameter of the linear equation of the travel time of the ith Rayleigh wave for the lth grid parameter, m l represents the density parameter of the lth grid, R α (z j ) and R β (z i ) represent the partial derivative of the longitudinal wave velocity Vp with respect to density and the partial derivative of the shear wave velocity Vs with respect to density, which can be expressed as:

[0075]

[0076] and

[0077]

[0078] In this embodiment, formula (5) establishes a linear equation between the Rayleigh surface wave travel time and the density of the formation to be detected. By using the above formula (5), a group of linear equations of all Rayleigh surface wave travel times and the density of the formation to be detected is established, and formula (2) of the first relationship between the Rayleigh surface wave travel time and the formation density can be obtained.

[0079] In this embodiment, by establishing the formation wave velocity structure based on seismic surface waves, inverting the density structure of the formation to be detected based on surface wave dispersion data, and establishing the first relationship between the Rayleigh surface wave travel time and the formation density, the internal structure of the shallow formation can be detected by remote sensing using non-destructive, non-invasive devices without any interference to the shallow formation structure, with the characteristics of high accuracy, high resolution and large observable depth.

[0080] In step S100, the muon data of the formation to be detected is obtained, including: arranging a muon observation system in a second selected area of ​​the formation to be detected, and obtaining the muon data observed by the muon observation system.

[0081] The muon observation system consists of muon imagers arranged according to certain rules in the second selected area. The second selected area can be selected according to actual detection needs. The muon imager includes: a muon position sensitive detector, an electronics system, and a data acquisition and transmission system.

[0082] The muon imager in the muon observation system can be set below or to the side of the object to be observed, and is used to receive the muon information remaining after the muons in the universe pass through the object to be observed. If the muon imager is placed to the side of the object to be observed, 360 / γ muon imagers can be arranged in a circle around the object to be observed at equal angles γ. The distance D between the muon position sensitive detector and the object to be observed is calculated by the detection angle resolution α and the required imaging position resolution Y:

[0083] D=Y / tanα (8)

[0084] The object to be observed in the present application is a stratum, therefore, a group of muon imagers can be arranged at equal intervals below the stratum to be detected. For example, if the object to be observed is a stratum overlying a subway tunnel, a group of muon imagers can be arranged at equal intervals in the subway tunnel. If the muon imager is placed below the stratum to be detected, after selecting a starting position according to the detection needs, the muon imagers can be arranged at equal intervals according to the positioning of the positioning system, that is, the muon observation system is arranged in the second selected area of ​​the stratum to be detected, and the interval Gap of the arranged muon imagers can be expressed as:

[0085] Gap=(n+depth)*tanβ (9)

[0086] In formula (9), n represents the distance between the muon imager and the lower surface of the formation to be detected, depth represents the thickness of the formation to be detected, and β represents the detectable angle range of the muon imager.

[0087] Figure 2 Schematically shows a layout diagram of a muon imager in a muon observation system. Muon data includes first muon data and second muon data, and the muon imager arranged below the formation to be detected is used to observe the second muon data after passing through the formation to be detected. It can be understood that for the muon observation system arranged below the formation to be detected, it is necessary to additionally arrange a muon imager above the formation to be detected to receive muon information before passing through the formation to be detected. Therefore, a muon imager can be arranged in the same arrangement at the surface of the formation to be detected to observe the first muon data in the sky that has not passed through the formation.

[0088] The first muon data includes the distribution of the number of muons on the surface of the formation to be detected with respect to the azimuth, and the second muon data includes the distribution of the number of muons at each muon sampling point under the formation to be detected with respect to the azimuth. The muon sampling point is the location of each muon imager in the muon observation system under the formation.

[0089] like Figure 4 As shown, in step S400, the density length of the formation to be detected is calculated based on the muon data, including:

[0090] S410, calculating the muon attenuation of each muon sampling point based on the first muon data and the second muon data.

[0091] S420, calculating the minimum energy required for muons to penetrate the formation to be detected based on the relationship between the muon flux and the muon attenuation at each muon sampling point.

[0092] S430, based on the relationship between the minimum energy required for the muon to penetrate the formation to be detected and the density length of the formation to be detected that the muon passes through, calculate the density length of the formation to be detected that the muon passes through.

[0093] In step S410, by extracting the position coordinate information of each muon track, the number of muons on the surface of the formation to be detected can be obtained as a function of the zenith angle θ and the azimuth angle. Distribution of directions: And the number of muons at the qth (q is a positive integer) muon sampling point in the muon observation system below the formation to be detected changes with the zenith angle θ and the azimuth angle Distribution of directions:

[0094] According to the distribution information of the number of muons with the first muon data and the second muon data, the observed value of the muon attenuation at each observation point can be obtained as follows:

[0095]

[0096] In step S420, the ratio of the muon flux Φ(θ, E) before and after passing through the formation to be detected is:

[0097]

[0098] In formula (11), E min represents the minimum energy required for the muon to penetrate the formation to be detected, Φ(θ, E) represents the muon flux, E represents the muon energy, and θ represents the muon zenith angle.

[0099] Among them, the formula for the muon flux Φ(θ, E) is:

[0100]

[0101]

[0102] In equations (12), (13) and (14), Φ(θ, E) represents the muon flux, in units of cm -2 s -1 sr -1 , sr is the unit of solid angle, E represents the muon energy in GeV, θ represents the muon zenith angle, p1, p2, p3, p4, and p5 are all parameters, and their parameter values ​​are: p1 = 0.102573, p2 = -0.068287, p3 = 0.958633, p4 = 0.0407253, and p5 = 0.817285.

[0103] Therefore, according to the above formulas (10)-(14), the minimum energy E required for muons to penetrate the formation to be detected can be obtained: min .

[0104] In step S430, according to the minimum energy E required for muons to penetrate the formation to be detected, min Relationship with the density length a of the formation region that the muon passes through:

[0105]

[0106] In formula (15), E min It represents the minimum energy required for a muon to penetrate the formation to be detected, E μ represents the rest energy of the muon, i.e. 105.66 MeV, and a represents the density length of the muon penetrating the formation to be detected.

[0107] The relationship between the muon cutoff energy and the range in the formation medium in the existing data and the muon observation data in the muon observation system can be used to set the preset average density and length of the formation to be detected. According to the preset average density and length of the formation to be detected, the initial response relationship between the density length of the formation to be detected and the muon attenuation ratio and the muon incident direction θ can be established (the number of muons on the surface does not change with the azimuth angle). change), Figure 5 The relationship between the density length of the formation to be detected and the muon attenuation ratio at several zenith angles is shown in FIG. Based on the known topographic data, the average density length a of the muon penetrating the formation to be detected at a specific azimuth and zenith angle can be obtained. av for:

[0108] a av =∫ Ω ∫ r ρdrdΩ (16)

[0109] In formula (16), r is the muon penetration path, which can be obtained based on the geometric structure of the formation area to be detected, Ω represents the circumferential area of ​​the formation area to be detected, and ρ represents the average density value in the direction of a muon path.

[0110] In step S500, determining a second relationship between the density length of the formation to be detected and the formation density includes:

[0111] Establish a linear equation group of density length and formation density of the formation to be detected:

[0112] (w μ L)δρ=(w μ ΔDL) (17)

[0113] In formula (17), L represents the partial derivative matrix of the density length of the formation to be detected that the muon passes through with respect to the formation density anomaly, δρ represents the density perturbation value of the three-dimensional grid of the formation density structure to be detected, ΔDL represents the residual of the density length of the formation to be detected that the muon passes through, and w μ represents the normalization weight factor for the muon data.

[0114] The following describes the steps of establishing a linear equation group of density length and formation density of the formation to be detected.

[0115] In this embodiment, according to the average density length a of the muon penetrating the formation area to be detected at a specific azimuth and zenith angle, av , it can be obtained that the density of the area of ​​the formation to be detected where the muons penetrate at this specific azimuth and zenith angle is:

[0116]

[0117] In formula (18), ρ av is the average density value in the formation area to be detected covered by the muon path and solid angle range, r is the muon penetration path, which can be obtained based on the geometric structure of the formation area to be detected, Ω represents the circumferential area of ​​the formation area to be detected, and ρ represents the average density value in a certain muon path direction.

[0118] The basic principle of muon imaging is similar to the eikonal equation of seismic Rayleigh surface wave travel time tomography. Therefore, according to formula (16) and formula (18), the partial derivative relationship between muon data and the three-dimensional density data of the formation to be detected can be established by referring to the seismic Rayleigh surface wave travel time imaging method. The structure of the formation to be detected is parameterized and discretized into several grid units. Based on formula (16) and formula (18), the density length of different zenith angles and azimuth angles is given by the muon detector at the qth (q is a positive integer) muon sampling point, and the linear equation group of the density length of the formation to be detected and the formation density is constructed:

[0119] DL q =∑ j L qj ρ j (19)

[0120] In formula (19), DL q represents the density length of the formation area to be detected that the muon at the qth muon sampling point passes through, L qj represents the length of the jth grid that the muon at the qth muon sampling point passes through, ρ j Represents the average density value of the j-th grid.

[0121] Formula (19) represents the linear equation of the density length and formation density of the formation to be detected obtained based on the muon observation data of the muon detector at a muon sampling point. Therefore, based on the above formula (19), the linear equation group of the density length and formation density of the formation to be detected can be established according to the muon observation data of the muon detector at all muon sampling points. Set the average density value ρ of the jth grid j The density perturbation value δρ of the jth grid j and the background density of the jth grid Composition, namely:

[0122]

[0123] Substituting formula (20) into formula (19) yields:

[0124]

[0125] Defines the background density length of the jth grid for:

[0126]

[0127] Therefore, calculate the density length residual of the formation to be detected through which the muon passes:

[0128]

[0129] The density length residual ΔDL of the formation to be detected that the muon passes through at all muon sampling points and the density perturbation value δρ of the three-dimensional grid of the density structure of the formation to be detected are expressed in the form of vectors and matrices, and the normalized weight factor w of the muon data is introduced considering the weights of different muon data. μ , we can get the formula (17) for the second relationship between the density length of the formation to be detected and the formation density.

[0130] In this embodiment, density distribution is established through muon imaging, the density structure of the formation to be detected is inverted based on muon observation data, and a second relationship between the density length of the formation to be detected and the formation density is established. Without any interference with the shallow formation structure, the internal structure of the shallow formation can be detected by remote sensing using non-destructive, non-invasive devices, and has the characteristics of high accuracy, high resolution and large observable depth.

[0131] Step S600, determining the three-dimensional density structure of the formation according to the first relationship and the second relationship, comprising:

[0132] Based on the linear equations of Rayleigh surface wave travel time and formation density (Formula 2), and the linear equations of density length of the formation to be detected and formation density (Formula 17), a linear equation group for inverting the three-dimensional density data of the formation based on Rayleigh surface wave travel time and density length of the formation to be detected is established:

[0133]

[0134] In formula (24), δρ represents the density perturbation value of the three-dimensional grid of the density structure of the formation to be detected, L represents the partial derivative matrix of the density length of the formation to be detected that the muon passes through with respect to the formation density anomaly, A represents the partial derivative matrix of the Rayleigh surface wave travel time with respect to the formation density anomaly, ΔDL represents the residual of the density length of the formation to be detected that the muon passes through, Δt represents the residual of the Rayleigh surface wave travel time, and w μ represents the normalized weight factor of the muon data, w s represents the normalized weight factor of seismic data, S represents the regularized smooth constraint, I represents the positive definite diagonal matrix, σ represents the normalization coefficient, and γ represents the damping factor. The regularized smooth constraint S and the positive definite diagonal matrix I are used to ensure the mathematical stability of the solution of the linear equation system, suppress local singular values, ensure the smoothness and rationality of the obtained model, and ensure that a solution with clear physical meaning and a reliable underground density structure model are obtained.

[0135] The density disturbance value of the three-dimensional grid of the stratum density structure to be detected is determined.

[0136] The three-dimensional density structure of the formation is determined according to the density disturbance value and a preset one-dimensional density model of the formation.

[0137] Since the partial derivative matrix L of the density length of the formation to be detected that the muon passes through and the partial derivative matrix A of the Rayleigh surface wave travel time for the formation density anomaly are both sparse matrices, formula (24) is proposed to be iteratively solved using the damped least squares method (LSQR) suitable for sparse matrices. The solution δρ obtained is the density perturbation value of the three-dimensional grid of the density structure of the formation to be detected. By superimposing δρ with the initial one-dimensional density model of the formation to be detected, the three-dimensional density structure model of the underground formation can be obtained. The initial one-dimensional density model is a density model constructed based on prior information such as local geological conditions or known data.

[0138] The present application provides a method for detecting the three-dimensional density structure of a stratum, which is based on the establishment of the stratum wave velocity structure by seismic surface waves and the establishment of the density distribution by muon imaging for joint inversion and solution, and obtains the accurate three-dimensional density structure of the stratum to be detected, which effectively improves the resolution and reliability of density structure imaging. This method combines the advantages of joint imaging of seismic surface wave data and muon data, breaks free from the constraints of traditional geophysical methods, and gets rid of the limitations of traditional single detection methods. It can detect the internal structure of shallow strata by non-destructive, non-invasive devices and remote sensing without any interference to the shallow stratum structure, and has the characteristics of high accuracy, high resolution and large observable depth.

[0139] For ease of understanding, a specific embodiment is given below to describe the method for detecting the three-dimensional density structure of the stratum of the present application. Figure 6 In this embodiment, the method for detecting the three-dimensional density structure of the stratum in this application is used to observe the overburden of the subway tunnel.

[0140] S601. Arrange a dense seismic array in a first selected area of ​​a to-be-detected stratum and obtain seismic background noise data observed by the dense seismic array.

[0141] After arriving at the subway tunnel to be observed, the observation starting position is first selected, and the position of the seismometer is determined at intervals of 100 meters using a total station electronic rangefinder. Then the seismometer is arranged according to the determined position, powered on, and the seismometer begins to collect and record continuous seismic background noise data.

[0142] S602. Extracting Rayleigh surface wave signals from seismic background noise data based on waveform cross-correlation method.

[0143] S603. Use a multiple filtering method to obtain Rayleigh surface wave dispersion data.

[0144] S604. Calculate the Rayleigh surface wave travel time based on the Rayleigh surface wave dispersion data.

[0145] The travel time t of a Rayleigh surface wave with a frequency of ω can be expressed as:

[0146]

[0147] Where l represents the propagation path of the Rayleigh surface wave, and S(l,ω) represents the slowness of the Rayleigh surface wave with the propagation path l and frequency ω.

[0148] S605. Establish a linear equation system of Rayleigh surface wave travel time and formation density:

[0149] (w s A)δρ=(w s Δt)

[0150] Where A represents the partial derivative matrix of Rayleigh surface wave travel time for formation density anomaly, δρ represents the density perturbation value of the three-dimensional grid of the formation density structure to be detected, Δt represents the residual of Rayleigh surface wave travel time, and w s Represents the normalized weight factor for seismic data.

[0151] S606. Arrange a muon observation system in a second selected area of ​​the formation to be detected, and obtain muon data observed by the muon observation system.

[0152] According to the selected observation starting position in the subway tunnel, the placement of the muon detector is determined at 20-meter intervals using a total station electronic rangefinder. Then the muon detector is arranged according to the determined position, powered on, and the muon detector begins to collect and record the track position data of muons. A muon detector is arranged on the surface in the same way to collect muon data in the sky that does not pass through the tunnel overburden.

[0153] S607, calculating the muon attenuation of each muon sampling point based on the first muon data and the second muon data.

[0154] S608. Calculate the minimum energy required for muons to penetrate the formation to be detected based on the relationship between the muon flux and the muon attenuation at each muon sampling point.

[0155] S609, based on the relationship between the minimum energy required for the muon to penetrate the formation to be detected and the density length of the formation to be detected that the muon passes through, calculate the density length of the formation to be detected that the muon passes through.

[0156] S610, establish a linear equation group of density length and formation density of the formation to be detected:

[0157] (w μ L)δρ=(w μ ΔDL)

[0158] Where L represents the partial derivative matrix of the density length of the formation to be detected that the muon passes through for the formation density anomaly, δρ represents the density perturbation value of the three-dimensional grid of the formation density structure to be detected, ΔDL represents the residual of the density length of the formation to be detected that the muon passes through, and w μ represents the normalization weight factor for the muon data.

[0159] S611, according to the linear equation group of Rayleigh surface wave travel time and formation density and the linear equation group of density length and formation density of the formation to be detected, a linear equation group for inverting the three-dimensional density data of the formation based on the Rayleigh surface wave travel time and the density length of the formation to be detected is established:

[0160]

[0161] Among them, δρ represents the density perturbation value of the three-dimensional grid of the density structure of the formation to be detected, L represents the partial derivative matrix of the density length of the formation to be detected that the muon passes through for the formation density anomaly, A represents the partial derivative matrix of the Rayleigh surface wave travel time for the formation density anomaly, ΔDL represents the residual of the density length of the formation to be detected that the muon passes through, Δt represents the residual of the Rayleigh surface wave travel time, w μ represents the normalized weight factor of the muon data, w s represents the normalized weight factor of seismic data, S represents the regularized smoothness constraint, I represents the positive definite diagonal matrix, σ represents the normalization coefficient, and γ represents the damping factor.

[0162] S612. Use the damped least squares method suitable for sparse matrices to iteratively solve the problem. The solution δρ is the density disturbance value of the three-dimensional grid of the density structure of the formation to be detected. Superimposing δρ with the initial one-dimensional density model of the formation to be detected can obtain the three-dimensional density structure model of the underground formation.

[0163] The present application provides a method for detecting the three-dimensional density structure of a stratum, which is based on the establishment of the stratum wave velocity structure by seismic surface waves and the establishment of the density distribution by muon imaging for joint inversion and solution, and obtains the accurate three-dimensional density structure of the stratum to be detected, which effectively improves the resolution and reliability of density structure imaging. This method combines the advantages of joint imaging of seismic surface wave data and muon data, breaks free from the constraints of traditional geophysical methods, and gets rid of the limitations of traditional single detection methods. It can detect the internal structure of shallow strata by non-destructive, non-invasive devices and remote sensing without any interference to the shallow stratum structure, and has the characteristics of high accuracy, high resolution and large observable depth.

[0164] Those skilled in the art will readily appreciate other embodiments of the present invention after considering the specification and practicing the invention disclosed herein. This application is intended to cover any variations, uses or adaptations of the present invention that follow the general principles of the present invention and include common knowledge or customary techniques in the art that are not disclosed in this application. The specification and examples are to be considered as exemplary only, and the true scope and spirit of the present invention are indicated by the following claims.

[0165] It should be understood that the present invention is not limited to the exact construction that has been described above and shown in the drawings and that various modifications and changes may be made without departing from the scope thereof. The scope of the present invention is limited only by the appended claims.

Claims

1. A method for detecting the three-dimensional density structure of a stratum, characterized in that: include: Obtain seismic background noise data and muon data of the formation to be detected; Based on the seismic background noise data, Rayleigh surface wave dispersion data is extracted, and Rayleigh surface wave travel time is calculated; Determining a first relationship between the Rayleigh surface wave travel time and formation density; Based on the muon data, calculating the density length of the formation to be detected; Determine a second relationship between the density length of the to-be-detected formation and the formation density; The three-dimensional density structure of the formation is determined according to the first relationship and the second relationship.

2. The method for detecting the three-dimensional density structure of a stratum according to claim 1, characterized in that: The extracting Rayleigh surface wave dispersion data based on the seismic background noise data and calculating the Rayleigh surface wave travel time comprises: Extracting Rayleigh surface wave signals from the seismic background noise data based on a waveform cross-correlation method; A multiple filtering method is used to obtain the Rayleigh surface wave dispersion data; The Rayleigh surface wave travel time is calculated based on the Rayleigh surface wave dispersion data.

3. The method for detecting the three-dimensional density structure of a stratum according to claim 2, characterized in that: Determining the first relationship between the Rayleigh surface wave travel time and the formation density includes: The linear equations of the Rayleigh surface wave travel time and the formation density are established: (w s A)δρ=(w s Δt) Where A represents the partial derivative matrix of the Rayleigh surface wave travel time for the formation density anomaly, δρ represents the density disturbance value of the three-dimensional grid of the formation density structure to be detected, Δt represents the Rayleigh surface wave travel time residual, w s Represents the normalized weight factor for seismic data.

4. The method for detecting the three-dimensional density structure of a stratum according to claim 1, characterized in that: The muon data includes first muon data and second muon data; The first muon data includes: the distribution of the number of muons on the surface of the formation to be detected with respect to the direction; and The second muon data includes: the distribution of the number of muons at each muon sampling point under the formation to be detected with respect to the direction.

5. The method for detecting the three-dimensional density structure of a stratum according to claim 4, characterized in that: The step of calculating the density length of the formation to be detected based on the muon data includes: Calculating the muon attenuation of each muon sampling point based on the first muon data and the second muon data; Calculating the minimum energy required for the muon to penetrate the formation to be detected according to the relationship between the muon flux and the muon attenuation at each muon sampling point; Based on the relationship between the minimum energy required for the muon to penetrate the formation to be detected and the density length of the formation to be detected that the muon passes through, the density length of the formation to be detected that the muon passes through is calculated.

6. The method for detecting the three-dimensional density structure of a stratum according to claim 5, characterized in that: The determining of the second relationship between the density length of the to-be-detected formation and the formation density comprises: A linear equation group of the density length of the formation to be detected and the formation density is established: (w μ L)δρ=(w μ (DL) Wherein, L represents the partial derivative matrix of the density length of the formation to be detected that the muon passes through for the formation density anomaly, δρ represents the density perturbation value of the three-dimensional grid of the density structure of the formation to be detected, ΔDL represents the residual of the density length of the formation to be detected that the muon passes through, and w μ represents the normalized weight factor of the muon data.

7. The method for detecting the three-dimensional density structure of a stratum according to claim 3 or 6, characterized in that: Determining the three-dimensional density structure of the formation according to the first relationship and the second relationship includes: Combining the first relationship and the second relationship, a linear equation group is obtained for inverting the three-dimensional density structure of the formation based on the Rayleigh surface wave travel time and the density length of the formation to be detected: Wherein, δρ represents the density disturbance value of the three-dimensional grid of the density structure of the formation to be detected, L represents the partial derivative matrix of the density length of the formation to be detected through which the muon passes for the formation density anomaly, A represents the partial derivative matrix of the Rayleigh surface wave travel time for the formation density anomaly, ΔDL represents the residual of the density length of the formation to be detected through which the muon passes, Δt represents the residual of the Rayleigh surface wave travel time, w μ represents the normalized weight factor of the muon data, w s represents the normalized weight factor of seismic data, S represents the regularized smoothing constraint, I represents the positive definite diagonal matrix, σ represents the normalization coefficient, and γ represents the damping factor; Determining the density disturbance value of the three-dimensional grid of the stratum density structure to be detected; The three-dimensional density structure of the formation is determined according to the density disturbance value and a preset one-dimensional density model of the formation.

8. The method for detecting the three-dimensional density structure of a stratum according to claim 1, characterized in that: The step of obtaining seismic background noise data of the formation to be detected comprises: Arrange a dense seismic array in a first selected area of ​​the stratum to be detected; The seismic background noise data observed by the seismic dense array is obtained.

9. The method for detecting the three-dimensional density structure of a stratum according to claim 1, characterized in that: The step of obtaining muon data of the formation to be detected includes: Arranging a muon observation system in a second selected area of ​​the formation to be detected; The muon data observed by the muon observation system is obtained.

10. The method for detecting the three-dimensional density structure of a stratum according to claim 9, characterized in that: The arranging of the muon observation system in the second selected area of ​​the to-be-detected formation comprises: Muon imagers are arranged at equal intervals below the formation to be detected, wherein the interval Gap is: Gap=(n+depth)*tanβ Wherein, n represents the distance between the muon imager and the lower surface of the formation to be detected, depth represents the thickness of the formation to be detected, and β represents the detectable angle range of the muon imager.

Images