Method for Identifying Deep Geological Anomalies Based on Transient Electromagnetic and Muon Detection

Abstract

The present invention discloses a method for identifying deep geological anomalies based on transient electromagnetic and muon detection. Transient electromagnetic detection and muon detection are respectively carried out on the area to be detected, and the resistivity distribution data of the area to be detected are obtained after processing the data obtained by transient electromagnetic detection; at the same time, the density distribution data of the area to be detected are obtained after processing the data obtained by muon detection; then, a cross-gradient algorithm based on discrete smooth interpolation constraint is used to jointly invert the data of the two physical properties. The advantage of the cross-gradient function lies in that it does not rely on the relationship between physical properties, but performs inversion through the similarity of the underground structure; the discrete smooth interpolation algorithm can suppress local oscillations and noise effects to a certain extent by performing smooth constraints on the inversion results of the cross-gradient function, making the final inversion result more accurate and realizing the accurate identification of density anomalies and resistivity anomalies in the area to be detected.

Description

Technical Field

[0001] The present invention relates to a method for identifying deep geological anomalies based on transient electromagnetic and muon detection, belonging to the technical field of geological exploration. Background Art

[0002] Joint inversion has always been a new idea in geophysical inversion. Its principle is to use different types of geophysical parameters for joint constrained inversion to improve the reliability and accuracy of the inversion results. Although the transient electromagnetic technique has been widely applied in the fields of in-tunnel detection, mineral and groundwater resource exploration, etc., due to the receiving points being located inside the tunnel, the limitations of spatial and temporal observations lead to inaccurate positioning of the internal structure of the tunnel, and inevitable observation errors, so there will be a problem of non-uniqueness in the inversion process.

[0003] Cosmic ray muons (μ-mesons) are a kind of natural rays. Based on their characteristics such as wide energy range and strong penetration ability, muon detection technology has been widely applied in multiple fields. Among them, muon imaging inversion through plastic scintillator detectors can detect the density distribution inside the tunnel. Its principle is to construct a muon scattering point cloud through the PoCA algorithm, determine the μ-meson trajectory through the information of the incident ray, scattering point, and outgoing ray, calculate the scattering density using the parameter unbiased estimation method, and obtain the muon imaging result through image reconstruction, so as to realize the density inversion of the imaging area inside the tunnel and obtain the density distribution map inside the tunnel. The traditional density distribution map is usually obtained by gravity detection means. Due to the relatively sparse points of gravity detection, the inverted density distribution map will have problems such as low resolution and non-uniqueness, resulting in large density errors. The sensitivity of muon detection technology to the density of rock mass can reduce the density data error and present a better inversion effect. However, at present, only the single method of muon detection cannot effectively improve the inversion accuracy, and there will be noise interference when muon detection is jointly inverted with other detection methods, ultimately resulting in low joint inversion accuracy.

[0004] Therefore, how to jointly invert the transient electromagnetic detection data and muon detection data, combine the advantages of both and reduce the noise interference in the joint inversion process, and ultimately effectively improve the geological exploration accuracy is the research direction required by the present invention. Summary of the Invention

[0005] Aiming at the problems existing in the above-mentioned prior art, the present invention provides a method for identifying deep geological anomalies based on transient electromagnetic and muon detection, which jointly inverts the transient electromagnetic detection data and muon detection data, combines the advantages of both and can reduce the noise interference in the joint inversion process, and ultimately effectively improves the geological exploration accuracy.

[0006] To achieve the above object, the technical solution adopted by the present invention is: a method for identifying deep geological anomalies based on transient electromagnetic and muon detection, and the specific steps are as follows:

[0007] Step 1: Deploy muon detectors in the roadway to conduct muon detection on the geological anomalies in the required detection area around the roadway; drill holes from the roadway to the required detection area, connect the transient electromagnetic probe to the transient electromagnetic host, and then deploy the transient electromagnetic probe in the drill holes to conduct transient electromagnetic detection on the geological anomalies in the required detection area;

[0008] Step 2: Stop the transient electromagnetic probe after it penetrates a certain distance in the drill hole. At this time, start the transient electromagnetic host and rotate the transient electromagnetic probe to conduct an in-situ transient electromagnetic detection once, and return the detected data to the transient electromagnetic host until multiple in-situ transient electromagnetic detections of the entire drill hole depth are completed. Then, the transient electromagnetic host preprocesses the collected data, and performs transient electromagnetic resistivity inversion on the processed induced electromotive force data through the Occam algorithm to obtain the resistivity distribution data of the required detection area;

[0009] Step 3: The muon detector collects data on the required detection area based on the muon angle scattering imaging principle to obtain the distribution of muon scattering points inside the roadway; then, construct a muon scattering point cloud from the data collected by the muon detector through the PoCA algorithm, determine the muon trajectory based on the incident ray, scattering point, and outgoing ray information, calculate the scattering density using the parameter unbiased estimation method, and obtain the muon imaging result through image reconstruction to achieve density inversion of the required detection imaging area and obtain the density distribution data of the required detection area;

[0010] Step 4: Use the cross-gradient algorithm based on discrete smooth interpolation constraints to jointly invert the resistivity distribution data obtained in Step 2 and the density distribution data obtained in Step 3. By coupling the two physical properties based on discrete smooth interpolation constraints and cross-gradient, the correlation and structure between different inversion results are enhanced, and the discrete smooth interpolation algorithm can reduce the local oscillation and noise effects in the cross-gradient algorithm processing, realize the identification of density anomalies and resistivity anomalies in the required detection area, and finally perform three-dimensional imaging on the required detection area using the digital twin system platform according to the joint inversion result.

[0011] Furthermore, the data collection in Step 3 includes obtaining the incident track data of muons in the muon incident imaging area, the outgoing track data of muons exiting the imaging area, and the coordinate information of muons at different positions using the muon detector; where the muon scattering angle is constrained by the Gaussian distribution of formulas (1) to (3);

[0012] Gaussian distribution of scattering angle:

[0013] (1)

[0014] is the standard deviation of the Gaussian distribution;

[0015] Radiation length calculation formula:

[0016] (2)

[0017] Where A is the relative atomic mass, Z is the atomic number, and ρ is the medium density;

[0018] Standard deviation of Gaussian distribution:

[0019] (3)

[0020] Where βc is the incident velocity of the muon, p is the incident momentum of the muon, and L is the incident track length.

[0021] Furthermore, the specific process of obtaining the density distribution data in the third step is as follows:

[0022] Before muon imaging, the detection area will be divided into a finite number of elemental units using the PoCA algorithm. Each unit is called a voxel. The variance of the scattering angle distribution is calculated using the scattering angle data obtained above. The specific formula is:

[0023] (4)

[0024] Where λ is the linear scattering density;

[0025] Substitute λ ≈ Zρ into formula (4) to get ;

[0026] Use the principle of unbiased estimation of parameters to perform unbiased estimation of the scattering density. The specific formula is:

[0027] (5)

[0028] Is the main scattering angle component on the vertical projection plane;

[0029] Assign the obtained scattering density value to the voxel where the scattering point is located. Through continuous accumulation of muon incident events, threshold classification is performed on the scattering density value to identify the contour range of the abnormal density body, and finally the density distribution data of the required detection area is obtained.

[0030] Furthermore, the specific process of the fourth step using the cross-gradient algorithm based on discrete smooth interpolation constraints for joint inversion is as follows:

[0031] The cross-gradient function for physical property constraints on the resistivity distribution data obtained in the second step and the density distribution data obtained in the third step is defined as:

[0032] (6)

[0033] represent the density and resistivity involved in cross-gradient calculation; where the magnitude of the t value can reflect the structural similarity between the two at the same location, and the smaller the t value, the more similar the structures;

[0034] Use the discrete smooth interpolation algorithm to interpolate and constrain the partial derivatives in the gradient function, and substitute it to deduce:

[0035]

[0036]

[0037] Then, form an inversion objective function by adding a constraint term:

[0038]

[0039] In the formula, α is the introduced discrete factor;

[0040] Introduce the partial derivative matrix A of the cross-gradient with cross-gradient values and discrete smooth interpolation constraints and solve to obtain:

[0041]

[0042] Among them, G is the sensitivity matrix and M is the weighting matrix;

[0043] Finally, through the above process for joint inversion, the identification of density anomalies and resistivity anomalies in the required detection area is achieved.

[0044] Compared with the prior art, the present invention respectively uses transient electromagnetic detection and muon detection for the required detection area, and after processing the data obtained from the transient electromagnetic detection, the resistivity distribution data of the required detection area is obtained; at the same time, after processing the data obtained from the muon detection, the density distribution data of the required detection area is obtained; then, the cross-gradient algorithm based on discrete smooth interpolation constraints is used for joint inversion of the data of the two physical properties (i.e., resistivity and density). The advantage of the cross-gradient function is that it does not depend on the relationship between physical properties, but performs inversion through the similarity of the underground structure; however, there will be local oscillations and noise effects during the joint inversion process of resistivity and density, while the discrete smooth interpolation algorithm can suppress local oscillations and noise effects to a certain extent by smoothing the results of the constraints, making the final inversion result more accurate and achieving the precise identification of density anomalies and resistivity anomalies in the required detection area. Through a specific algorithm for joint inversion of transient electromagnetic detection data and muon detection data, combining the advantages of both and reducing noise effects and local oscillations during the joint inversion process, the geological detection accuracy is ultimately effectively improved. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Figure 1 is the overall layout schematic diagram of the present invention;

[0046] Figure 2 It is the anomaly identification map after the joint inversion of the present invention.

[0047] In the figure: 1 - transient electromagnetic probe, 2 - connecting wire, 3 - transient electromagnetic mainframe, 4 - borehole, 5 - muon detector, 6 - muon ray. Specific implementation manner

[0048] The present invention will be further described below.

[0049] As Figure 1 shown, the specific steps of the present invention are as follows:

[0050] Step 1: Install the muon detector 5 in the roadway to perform muon detection on the geological anomaly conditions in the required detection area around the roadway; construct a borehole 4 from the roadway to the required detection area, connect the transient electromagnetic probe 1 and the transient electromagnetic mainframe 3 through the connecting wire 2, and then install the transient electromagnetic probe 1 in the borehole 4 to perform transient electromagnetic detection on the geological conditions in the required detection area;

[0051] Step 2: Stop the transient electromagnetic probe 1 after it penetrates a certain distance in the borehole 4. At this time, start the transient electromagnetic mainframe 3 and rotate the transient electromagnetic probe 1 to perform an in-situ transient electromagnetic detection once, and return the detected data to the transient electromagnetic mainframe 3 until multiple in-situ transient electromagnetic detections of the entire depth of the borehole 4 are completed. Then, the transient electromagnetic mainframe 3 preprocesses the collected data, and performs transient electromagnetic resistivity inversion on the processed induced electromotive force data through the Occam algorithm to obtain the resistivity distribution data of the required detection area;

[0052] Step 3: The muon detector 5 collects data on the required detection area based on the muon angle scattering imaging principle, so as to obtain the distribution of muon scattering points inside the roadway; the data collection includes obtaining the incident track data of the muon incident imaging area, the outgoing track data exiting from the imaging area, and the coordinate information of the muons at different positions by using the muon detector 5; the muon scattering angle is constrained by the Gaussian distribution of formulas (1) to (3);

[0053] Gaussian distribution of scattering angle:

[0054] (1)

[0055] is the standard deviation of the Gaussian distribution;

[0056] Formula for calculating the radiation length:

[0057] (2)

[0058] Among them, A is the relative atomic mass, Z is the atomic number, and ρ is the medium density;

[0059] Standard deviation of Gaussian distribution:

[0060] (3)

[0061] Among them, βc is the muon incident velocity, p is the muon incident momentum, and L is the incident track length.

[0062] Next, the data collected by the muon detector 5 is used to construct a muon scattering point cloud through the PoCA algorithm. The muon trajectory is determined by the information of the incident ray, scattering point, and outgoing ray. The scattering density is calculated using the parameter unbiased estimation method, and the muon imaging result is obtained through image reconstruction, so as to realize the density inversion of the required detection imaging area and obtain the density distribution data of the required detection area. The specific process is as follows:

[0063] Before realizing muon imaging using the PoCA algorithm, the detection area will be divided into a finite number of element units first. Each unit is called a voxel. The variance of the scattering angle distribution is calculated through the scattering angle data obtained above. The specific formula is:

[0064] (4)

[0065] Among them, λ is the linear scattering density;

[0066] Substitute λ ≈ Zρ into formula (4) to get

[0067] ;

[0068] The unbiased estimation of the scattering density is carried out using the parameter unbiased estimation principle. The specific formula is:

[0069] (5)

[0070] Is the main scattering angle component on the vertical projection plane;

[0071] The obtained scattering density value is assigned to the voxel where the scattering point is located. Through the continuous accumulation of muon incident events, the scattering density value is classified by threshold, so as to identify the contour range of the abnormal density body, and finally obtain the density distribution data of the required detection area.

[0072] Step 4: Use the cross-gradient algorithm based on discrete smooth interpolation constraints to jointly invert the resistivity distribution data obtained in Step 2 and the density distribution data obtained in Step 3. By coupling the two physical properties after discrete smooth interpolation constraint cross-gradient, the correlation and structure between different inversion results are enhanced, and the discrete smooth interpolation algorithm can reduce the local oscillation and noise influence in the cross-gradient algorithm processing, realizing the identification of density anomalies and resistivity anomalies in the required detection area. The specific process is as follows:

[0073] The cross-gradient function with physical property constraints for the resistivity distribution data obtained in Step 2 and the density distribution data obtained in Step 3 is defined as:

[0074] (6)

[0075] represent the density and resistivity participating in the cross-gradient calculation; where the value of t can reflect the structural similarity degree of the two at the same position, and the smaller the t value, the more similar the structure;

[0076] Use the discrete smooth interpolation algorithm to perform interpolation constraints on the partial derivatives in the gradient function

[0077] After substituting, it is deduced as:

[0078]

[0079]

[0080]

[0081]

[0082] Then, form the inversion objective function by adding a constraint term:

[0083]

[0084] In the formula, α is the introduced discrete factor;

[0085] Introduce the cross-gradient value and the partial derivative matrix A of the cross-gradient with discrete smooth interpolation constraints and solve to obtain:

[0086]

[0087] where G is the sensitivity matrix and M is the weighting matrix;

[0088] Finally, through the above process for joint inversion, the identification of density anomalies and resistivity anomalies in the required detection area is realized as Figure 2 shown. Finally, according to the joint inversion results, use the digital twin system platform to perform three-dimensional imaging on the required detection area.

[0089] The above are only the preferred embodiments of the present invention. It should be noted that for those of ordinary skill in the art, without departing from the principle of the present invention, several improvements and refinements can be made, and these improvements and refinements should also be regarded as the protection scope of the present invention.

Claims

1. A method for identifying deep geological anomalies based on transient electromagnetic and muon detection, characterized in that: The specific steps are: Step 1: a muon detector is arranged in the tunnel to perform muon detection on the geological conditions of the required detection area around the tunnel; a hole is drilled from the tunnel to the required detection area, and a transient electromagnetic probe is connected to a transient electromagnetic host, and then the transient electromagnetic probe is arranged in the drill hole to perform transient electromagnetic detection on the geological conditions of the required detection area; Step 2: Stop the transient electromagnetic probe after it goes a certain distance into the borehole, start the transient electromagnetic host and rotate the transient electromagnetic probe to perform an in-situ transient electromagnetic detection, and return the data obtained by the detection to the transient electromagnetic host until multiple in-situ transient electromagnetic detections of the entire borehole depth are completed. Then the transient electromagnetic host pre-processes the collected data, performs transient electromagnetic resistivity inversion on the processed induced electromotive force data through the Occam algorithm, and obtains the resistivity distribution data of the required detection area; Step 3: The muon detector collects data on the required detection area based on the muon angle scattering imaging principle, so as to obtain the distribution of muon scattering points inside the tunnel; then, the data collected by the muon detector is used to construct a muon scattering point cloud through the PoCA algorithm, and the muon trajectory is determined by the incident ray, scattering point, and outgoing ray information. The scattering density is calculated using the parameter unbiased estimation method, and the muon imaging result is obtained through image reconstruction, thereby realizing the density inversion of the required detection imaging area and obtaining the density distribution data of the required detection area; Step 4: Use the cross-gradient algorithm based on discrete smooth interpolation constraints to jointly invert the resistivity distribution data obtained in step 2 and the density distribution data obtained in step 3. By coupling the two physical properties after the cross-gradient based on discrete smooth interpolation constraints, the correlation and structure between different inversion results are enhanced. In addition, the discrete smooth interpolation algorithm can reduce the influence of local oscillation and noise in the cross-gradient algorithm processing, and realize the identification of density anomalies and resistivity anomalies in the required detection area. Finally, the digital twin system platform is used to perform three-dimensional imaging of the required detection area based on the joint inversion results.

2. The method for identifying deep geological anomalies based on transient electromagnetic and muon detection according to claim 1 is characterized in that: The data acquisition in step 3 includes using a muon detector to obtain incident track data of muons entering the imaging area, exit track data emitted from the imaging area, and coordinate information of muons at different positions; wherein the muon scattering angle is constrained by the Gaussian distribution of formulas (1) to (3); Scattering angle Gaussian distribution: （1） Among them, σ θ is the standard deviation of Gaussian distribution; Radiation length calculation formula: （2） Among them, A is the relative atomic mass, Z is the atomic number, and ρ is the medium density; Gaussian distribution standard deviation: （3） Among them, βc is the muon incident velocity, p is the muon incident momentum, and L is the incident track length.

3. The method for identifying deep geological anomalies based on transient electromagnetic and muon detection according to claim 2 is characterized in that: The specific process of obtaining density distribution data in step 3 is as follows: Before using the PoCA algorithm to realize muon imaging, the detection area is divided into a finite number of element units. Each unit is called a voxel. The variance of the scattering angle distribution is calculated using the scattering angle data obtained above. The specific formula is: （4） Where, λ is the linear scattering density; Substituting λ≈Zρ into formula (4) yields ; The scattering density is estimated unbiasedly using the unbiased parameter estimation principle. The specific formula is: （5） Among them, θ x and θ y is the main scattering angle component on the vertical projection plane; The obtained scattering density value is assigned to the voxel where the scattering point is located. Through the continuous accumulation of muon incident events, the scattering density value is threshold-classified to identify the contour range of the abnormal density body, and finally the density distribution data of the required detection area is obtained.

4. The method for identifying deep geological anomalies based on transient electromagnetic and muon detection according to claim 1, characterized in that: The specific process of the joint inversion using the cross gradient algorithm based on discrete smooth interpolation constraints in step 4 is as follows: The cross gradient function for physical property constraints on the resistivity distribution data obtained in step 2 and the density distribution data obtained in step 3 is defined as: （6） In the formula, ∇ represents the gradient operation symbol; m1 and m2 represent the density and resistivity involved in the cross-gradient calculation; the size of the t value can reflect the structural similarity of the two at the same position. The smaller the t value, the more similar the structure. The discrete smooth interpolation algorithm is used to calculate the gradient function The partial derivatives in are interpolated and constrained, and then the following is derived: Then, the inversion objective function is formed by adding constraints: In the formula, α is the introduced discrete factor; Introducing the partial derivative matrix A of the cross gradient with cross gradient value and discrete smooth interpolation constraints and solving it yields: Among them, G is the sensitivity matrix and M is the weighting matrix; Finally, a joint inversion is performed through the above process to realize the identification of density anomalies and resistivity anomalies in the required detection area.

Images