Deep geological anomalous body identification method based on transient electromagnetism and muon detection

Abstract

The invention discloses a deep geological anomalous body identification method based on transient electromagnetic and muon detection, which comprises the following steps of: respectively carrying out transient electromagnetic detection and muon detection on a region needing to be detected, and processing data obtained by the transient electromagnetic detection to obtain resistivity distribution data of the region needing to be detected; meanwhile, processing data obtained by muon detection to obtain density distribution data of the area needing to be detected; joint inversion is carried out on the data of the two physical properties through a cross gradient algorithm based on discrete smooth interpolation constraint, and the advantage of a cross gradient function does not depend on a relational expression between the physical properties, but inversion is carried out through the similarity of the underground structure; the discrete smooth interpolation algorithm performs smooth constraint on the inversion result of the cross gradient function, so that the influence of local oscillation and noise can be suppressed to a certain extent, the precision of the final inversion result is higher, and accurate identification of density anomaly and resistivity anomaly in the area needing to be detected is realized.

Description

Technical Field

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

[0002] Joint inversion has always been a new idea for geophysical inversion. Its principle is to use different types of geophysical parameters to jointly constrain inversion in order to improve the reliability and accuracy of the inversion results. Although transient electromagnetic technology has been widely used in the fields of internal tunnel detection, mineral and groundwater resource exploration, the receiving point is located inside the tunnel. Due to the limitations of space and time observation, the positioning of the internal structure of the tunnel is not accurate, and observation errors are inevitable, which will lead to multiple solutions in the inversion process.

[0003] Cosmic ray muons (muons) are a type of natural ray. Based on their wide energy range and strong penetration, muon detection technology has been widely used in multiple fields. Among them, muon imaging inversion using plastic scintillator detectors can detect the density distribution inside the tunnel. The principle is to use the PoCA algorithm to construct a muon scattering point cloud, determine the muon trajectory through 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, thereby realizing the density inversion of the imaging area inside the tunnel and obtaining the density distribution map inside the tunnel. Traditional density distribution maps are usually obtained by gravity detection. Since the points of gravity detection are relatively sparse, 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 rock mass density can reduce the error of density data and present a better inversion effect. However, at present, muon detection alone cannot effectively improve the inversion accuracy. At the same time, the joint inversion of muon detection and other detection methods will be affected by noise, which ultimately leads to low joint inversion accuracy.

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

[0005] In view of 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 transient electromagnetic detection data and muon detection data, combines the advantages of both and can reduce the impact of noise in the joint inversion process, ultimately effectively improving the accuracy of geological detection.

[0006] In order to achieve the above-mentioned purpose, the technical solution adopted by the present invention is: a method for identifying deep geological anomalies based on transient electromagnetic and muon detection, the specific steps are: Step 1: a muon detector is arranged in the tunnel to perform muon detection on the geological anomalies in 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 anomalies in 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.

[0007] Furthermore, 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) 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.

[0008] Further, the specific process of obtaining density distribution data in step 3 is: 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) 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.

[0009] Furthermore, 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) Represents 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.

[0010] Compared with the prior art, the present invention adopts transient electromagnetic detection and muon detection for the required detection area, and obtains the resistivity distribution data of the required detection area after processing the data obtained by transient electromagnetic detection; at the same time, the density distribution data of the required detection area is obtained after processing the data obtained by muon detection; then the cross-gradient algorithm based on discrete smooth interpolation constraints is used to jointly invert the data of the two physical properties (i.e., resistivity and density), wherein the cross-gradient function has the advantage of not relying on the relationship between the physical properties, but inverting through the similarity of the underground structure; however, there will be local oscillation and noise influence in the joint inversion process of resistivity and density, and the discrete smooth interpolation algorithm can suppress the local oscillation and noise influence to a certain extent by smoothing the results, so that the final inversion result is more accurate, and the density anomaly and resistivity anomaly in the required detection area are accurately identified. Through a specific algorithm, the transient electromagnetic detection data and the muon detection data are jointly inverted, combining the advantages of both and reducing the noise influence and local oscillation in the joint inversion process, and finally effectively improving the accuracy of geological detection. BRIEF DESCRIPTION OF THE DRAWINGS

[0011] Figure 1 It is a schematic diagram of the overall layout of the present invention; Figure 2 This is an abnormal body identification diagram after joint inversion of the present invention.

[0012] In the figure: 1-transient electromagnetic probe, 2-connecting line, 3-transient electromagnetic host, 4-drilling hole, 5-muon detector, 6-muon ray. DETAILED DESCRIPTION

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

[0014] like Figure 1 As shown, the specific steps of the present invention are: Step 1: a muon detector 5 is arranged in the tunnel to perform muon detection on the geological anomalies in the required detection area around the tunnel; a borehole 4 is constructed from the tunnel to the required detection area, and a transient electromagnetic probe 1 is connected to a transient electromagnetic host 3 through a connecting line 2, and then the transient electromagnetic probe 1 is arranged in the borehole 4 to perform transient electromagnetic detection on the geological conditions in the required detection area; Step 2: The transient electromagnetic probe 1 is stopped after each penetration of a certain distance in the borehole 4. At this time, the transient electromagnetic host 3 is started and the transient electromagnetic probe 1 is rotated to perform an in-situ transient electromagnetic detection, and the data obtained by the detection is returned to the transient electromagnetic host 3 until multiple in-situ transient electromagnetic detections of the entire depth of the borehole 4 are completed. Then, the transient electromagnetic host 3 pre-processes 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; 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 tunnel; wherein the data collection includes using the muon detector 5 to obtain the incident track data of the muon incident imaging area, the exit track data emitted from the imaging area, and the coordinate information of the muon at different positions; wherein the muon scattering angle is constrained by the Gaussian distribution of formulas (1) to (3); Scattering angle Gaussian distribution: (1) 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.

[0015] Then, the PoCA algorithm is used to construct the muon scattering point cloud of the data collected by the muon detector 5. 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. 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. The specific process 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) 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.

[0016] 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. The two physical properties coupled after the cross-gradient based on discrete smooth interpolation constraints enhance the correlation and structure between different inversion results. In addition, the discrete smooth interpolation algorithm can reduce the local oscillation and noise effects in the cross-gradient algorithm processing, and realize the identification of density anomalies and resistivity anomalies in the required detection area. The specific process 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) Represents 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, the above process is used to perform joint inversion to identify density anomalies and resistivity anomalies in the required detection area. Figure 2 As shown, finally, the digital twin system platform is used to perform three-dimensional imaging of the required detection area based on the joint inversion results.

[0017] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principle of the present invention. These improvements and modifications should also be regarded as the scope of protection 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