Method for imaging background suppression of natural muons under complex surface environment with shield dynamic

Abstract

The application discloses a method for imaging background suppression of natural muon under shield dynamic in complex surface environment, which comprises the following steps: obtaining geological exploration data and surface building data, constructing a three-dimensional forward grid model, continuously collecting natural muon detection data in shield tunneling, and obtaining measured muon survival rate; extracting the direction interval which is not blocked by surface buildings to calculate the forward muon survival rate, and establishing a correlation model between the muon penetration energy and the equivalent mass thickness according to the forward and measured survival rates; decoupling and calibrating the equivalent density and geometric height of the surface building according to the correlation model; updating the calibrated parameters to the three-dimensional forward grid model, screening the detection data, and performing three-dimensional inversion of geological anomalies, and outputting the underground density distribution image. The application realizes the endogenous separation of surface building background interference and high-fidelity inversion of underground geological anomalies, and provides reliable support for shield tunneling geological disaster warning.

Description

Technical Field

[0001] This invention relates to the field of underground engineering geological exploration technology, and in particular to a method for suppressing the background of dynamic natural muon imaging under complex surface environments. Background Technology

[0002] With the acceleration of urbanization, urban rail transit systems are increasingly being built in densely populated cities, where underground space is occupied by various utility pipelines, building foundations, and heterogeneous fill. During shield tunneling, geological and structural risks are distributed throughout the entire project lifecycle, including preliminary geological exploration, on-site construction, and subsequent operation and maintenance. Hidden geological hazards, such as ground cavities or pipeline conflicts, may be amplified during excavation due to stress redistribution, ultimately manifesting as surface subsidence, pipeline damage, and other hazards. To prevent such situations, modern tunnel engineering often uses geophysical methods for advanced geological prediction (AGP). However, the complex environment of shield tunneling presents unique obstacles to the application of dynamic monitoring methods during drilling. Electromagnetic methods such as ground-penetrating radar (GPR) suffer severe attenuation in water-rich strata and are significantly affected by the shielding of the metal cutterhead and shield machine structure, making it difficult to guarantee the accuracy of the results; resistivity and transient electromagnetic methods are susceptible to stray currents from the city and electromagnetic interference caused by the shield machine itself; seismic or microseismic techniques are affected by vibrations during construction, resulting in low signal-to-noise ratios.

[0003] Natural muon imaging technology utilizes muons—secondary cosmic ray particles with extremely strong penetrating power—as a radiation source, offering a highly promising technical solution to the aforementioned problems. At sea level, the average energy of a muon is approximately 3-4 GeV, with a flux of about 10,000 m⁻²·min⁻¹. Based on the principle of muon-matter interaction, muon imaging technology can be divided into scattering imaging based on multiple Coulomb scattering and transmission imaging based on flux attenuation. The former is suitable for imaging studies of small to medium-sized objects, such as nuclear materials detection and concrete structures. The latter has already demonstrated excellent deep-sea detection performance in fields such as mineral exploration, archaeological excavation, and industrial production, providing a highly promising solution for these engineering needs. In recent years, the application of muon imaging in underground engineering has achieved initial empirical results. Thompson et al. used a portable, field-deployable natural muon detector to rapidly scan and image the entire overburden of a railway tunnel in the UK, discovering a potential cavity in the strata. Li Zhiwei et al. first verified the technology's ability to identify density anomalies within a 25 m overburden layer in a tunnel section of the Shenzhen Metro. Subsequently, they constructed an experimental system combining "münite detection + algorithm optimization" for the overburden layer of the tunnel between Apollo South Station and Yuanshan Station in the Shenzhen Metro, and obtained the density distribution of the overlying strata. SiangKhaw et al. used a mobile system to obtain the spatial distribution of muon flux above the Shanghai Outer Ring Tunnel, studying the impact of the overburden layer on muon flux attenuation.

[0004] However, existing technologies are mostly limited to static, fixed-point observations of existing tunnels or discrete surveys under simple surface conditions. When muon imaging technology truly evolves towards dynamic monitoring of tunnel boring machines (TBMs) during excavation and is applied to highly complex urban surface environments, it systematically exposes three major technological application bottlenecks and gaps in existing research: Spatiotemporal dynamic evolution contradiction (no existing solution): The low natural muon flux requires a long exposure time, which inherently contradicts the continuous spatial displacement of the tunnel boring machine. Existing static algorithms cannot handle the cumulative error in a dynamic displacement coordinate system.

[0005] False positive interference from heterogeneous surface buildings: Urban high-rise buildings introduce drastic non-uniform flux attenuation. Without an accurate background equivalence strategy, traditional inversion algorithms are prone to misinterpreting the mass load of surface buildings as density anomalies in underground soil and rock, resulting in artificial density artifacts in the corresponding zenith angle direction.

[0006] The challenge of multivariate physical coupling between height and density: Real buildings include heterogeneous floor slabs, beams, columns, and cavities, and their equivalent density fluctuates drastically with the incident angle of rays; moreover, there is a significant discrepancy between the number of floors on the construction drawings and the actual geometric height (often by several to tens of meters). Current technologies cannot effectively separate and independently calibrate the two variables that affect flux attenuation: the equivalent density and the geometric height of the building.

[0007] In summary, there is an urgent need for a systematic solution that can be integrated into the modern tunnel boring machine workflow. Summary of the Invention

[0008] This invention provides a method for suppressing the background of dynamic natural muon imaging under complex surface environments, which solves the problems of dynamic displacement and long exposure, false positive interference from surface buildings, and imaging artifacts caused by the difficulty in decoupling density-height coupling in the prior art. It achieves endogenous stripping of surface background interference and high-fidelity inversion of underground geological anomalies.

[0009] This invention provides a method for suppressing the background of dynamic natural muon imaging under complex surface environments, the method comprising: Obtain geological exploration data and surface building data of the target area, and construct a three-dimensional forward mesh model including underground geological structure and surface building geometric parameters and equivalent density; During the tunnel boring machine (TBM) excavation process, natural muon detection data is continuously collected. The natural muon detection data collected within a single segment assembly cycle is mapped to the spatial position of the TBM when it is stationary within the cycle, and the measured muon survival rate at the spatial position is obtained. Extract the directional intervals that are not obstructed by the surface buildings from the three-dimensional forward mesh model, calculate the forward muon survival rate within the directional intervals, and establish a correlation model between muon penetration energy and equivalent mass thickness based on the forward muon survival rate and the measured muon survival rate. Based on the association model, the equivalent density and geometric height of surface buildings in the three-dimensional forward mesh model are decoupled and calibrated. The calibrated equivalent density and geometric height are updated to the three-dimensional forward modeling grid. After filtering the natural muon detection data, three-dimensional geological anomaly inversion is performed to output a subsurface density distribution image.

[0010] In one possible implementation, the construction includes a three-dimensional forward mesh model comprising the underground geological structure, surface building geometry parameters, and equivalent density, including: A subsurface geological structure grid is established based on the geological exploration data from the geological borehole data. A surface building grid is established based on the engineering drawings in the surface building data, and the underground geological structure grid and the surface building grid are merged to form the three-dimensional forward model.

[0011] In one possible implementation, the continuous acquisition of natural muon detection data during shield tunneling, mapping the natural muon detection data acquired within a single segment assembly cycle to the spatial location of the shield machine when it is stationary within that cycle, and obtaining the measured muon survival rate at the spatial location, includes: The first muon detector is installed behind the cutterhead of the tunnel boring machine, and the second muon detector is installed on the supporting trolley behind the tunnel boring machine. Natural muon detection data are continuously collected during the shield tunneling process using the first muon detector and the second muon detector. The natural muon detection data collected within a single segment assembly cycle is equivalently mapped to the spatial position of the tunnel boring machine when it is in a static state during segment assembly within the cycle, and the mapped data is used as the static observation data at the spatial position. The measured muon flux is calculated based on the static observation data to obtain the measured muon flux under unobstructed aerial conditions. The measured muon flux is then normalized based on the measured muon flux to obtain the measured muon survival rate.

[0012] In one possible implementation, extracting directional intervals from the three-dimensional forward mesh model that are not obstructed by the surface buildings, and calculating the forward muon survival rate within the directional intervals, includes: Ray tracing is performed on the three-dimensional forward mesh model, traversing each directional interval, and determining whether there are mesh cells of the surface buildings on each ray path; The directional intervals in which no surface building grid cells are present on the ray path are marked as unobstructed directional intervals; Extract the formation density and thickness data in each direction within the unobstructed directional interval, and calculate the theoretical flux decay of muons after they pass through the formation using forward modeling to obtain the forward modeling muon survival rate in the unobstructed directional interval.

[0013] In one possible implementation, establishing a correlation model between muon penetration energy and equivalent mass thickness based on the forward muon survival rate and the measured muon survival rate includes: Based on the forward muon survival rate and the measured muon survival rate in each direction within the unobstructed direction interval, calculate the minimum penetration energy of the measured muon in each direction, and calculate the equivalent mass thickness in each direction based on the three-dimensional forward mesh model. Using the equivalent mass thickness as the independent variable and the measured minimum muon penetration energy as the dependent variable, the data in each direction within the unobstructed direction interval are fitted to establish a correlation model between muon penetration energy and equivalent mass thickness.

[0014] In one possible implementation, the decoupling calibration of the equivalent density and geometric height of surface buildings in the three-dimensional forward mesh model according to the associated model includes: The zenith angle direction range through which the ray path passes through the surface building is selected. Within the azimuth angle range blocked by the surface building, the equivalent density of the surface building is used as a variable, and the objective is to minimize the residual between the measured muon survival rate and the forward muon survival rate. The equivalent density is solved by iterative inversion. Based on the convergence of the equivalent density, the first decay inflection point of the measured muon survival rate as a function of the zenith angle and the second decay inflection point of the forward modeling muon survival rate as a function of the zenith angle are obtained. The geometric height of the surface buildings is adjusted until the first decay inflection point and the second decay inflection point coincide on the angle axis. Repeat the iterative inversion solution of the equivalent density and the iterative adjustment of the geometric height until both the equivalent density and the geometric height satisfy the preset convergence condition.

[0015] In one possible implementation, updating the calibrated equivalent density and geometric height to the 3D forward mesh model includes: Replace the initial equivalent density of the corresponding surface building in the three-dimensional forward mesh model with the calibrated equivalent density; Replace the initial geometric height of the corresponding surface building in the three-dimensional forward mesh model with the calibrated geometric height.

[0016] In one possible implementation, the step of filtering the natural muon detection data and then performing three-dimensional geological anomaly inversion to output a subsurface density distribution image includes: Forward simulation was performed based on the updated three-dimensional forward mesh model to obtain the forward muon survival rate in each directional interval; Calculate the deviation between the measured muon survival rate and the forward muon survival rate in each directional interval, compare the deviation with a preset deviation threshold, and discard the natural muon detection data corresponding to the directional intervals where the deviation is greater than the preset deviation threshold; Based on the three-dimensional forward mesh model, it is determined whether the ray path of each directional interval is blocked by the ground building. For the directional interval that is not blocked, the corresponding natural muon detection data is retained. For the directional interval that is blocked, the muon count in the directional interval is counted. When the muon count is less than a preset count threshold, the corresponding natural muon detection data is removed. The remaining natural muon detection data after screening are used to perform three-dimensional inversion of geological anomalies and output an image of underground density distribution.

[0017] In one possible implementation, the association model is a linear model, represented as follows: ; in, Indicates the first fitting coefficient; Indicates the second fitting coefficient; Indicates equivalent mass thickness; This represents the minimum penetrating energy of a muon.

[0018] One or more technical solutions provided in this invention have at least the following technical effects or advantages: This invention integrates preliminary geological drilling data and surface engineering drawings to construct a three-dimensional forward modeling mesh model, providing a unified benchmark model including underground geological structure and surface building parameters for subsequent muon imaging forward modeling calculations. This eliminates initial deviations caused by insufficient prior information or coarse models, ensuring the accuracy of background suppression from the source. By mapping natural muon detection data collected within a single tunnel segment assembly cycle to the spatial position of the tunnel boring machine (TBM) when it is stationary, this invention eliminates the cumulative spatial coordinate errors caused by continuous displacement of the TBM while ensuring the statistical quantity of muon data within a single cycle, making the static inversion algorithm applicable to dynamic monitoring scenarios during tunneling. Furthermore, by establishing muon penetration energy and... The correlation model between equivalent mass thickness eliminates the background error introduced by the non-uniform fluctuation of the energy loss rate of the formation medium, providing an accurate physical parameter benchmark for forward modeling. By decoupling and calibrating the equivalent density and geometric height of surface buildings, the coupling interference of surface building mass load and spatial height on muon attenuation is eliminated, preventing the inversion algorithm from misjudging surface buildings as underground soil and rock density anomalies. By updating the calibrated equivalent density and geometric height to the three-dimensional forward model and performing inversion after data filtering, low-confidence distorted ray data is eliminated, and false positive artifacts caused by residual errors in the engineering drawings are eliminated, outputting a high-fidelity underground density distribution image, providing reliable data support for geological disaster early warning in shield tunneling. Attached Figure Description

[0019] Figure 1 A flowchart illustrating the steps of the method for suppressing background in dynamic natural muon imaging under complex surface environments provided in this embodiment of the invention; Figure 2 A schematic diagram of shield tunneling stage detection provided in an embodiment of the present invention; Figure 3 This is a schematic diagram of the detector position provided in an embodiment of the present invention; Figure 4 A schematic diagram of the subway geology and the building model along the line provided in this embodiment of the invention; Figure 5a This is a schematic diagram illustrating the impact of low building density on survival rate, provided in an embodiment of the present invention. Figure 5b This is a schematic diagram illustrating the impact of medium building density on survival rate, provided in an embodiment of the present invention. Figure 5cThis is a schematic diagram illustrating the impact of high building density on survival rate, provided in an embodiment of the present invention. Figure 6 This is a schematic diagram of the physical model fitting results provided in an embodiment of the present invention; Figure 7 This is a schematic diagram of buildings along the route provided in an embodiment of the present invention; Figures 8a-8d This is a schematic diagram illustrating the synchronization analysis of the survival rate decay inflection point at the correct altitude provided in an embodiment of the present invention. Figure 8 e~ Figure 8h This is a schematic diagram illustrating the synchronization of the inflection point of survival rate decline after a 20m reduction in building height, as provided in an embodiment of the present invention. Detailed Implementation

[0020] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.

[0021] This invention provides a method for suppressing the background of dynamic natural muon imaging under complex surface environments, see [link to relevant documentation]. Figure 1 The method includes the following steps S101 to S105.

[0022] S101, acquire geological exploration data and surface building data of the target area, and construct a three-dimensional forward mesh model including underground geological structure and surface building geometric parameters and equivalent density; Specifically, in step S101, a three-dimensional forward mesh model is constructed, including the geometric parameters of underground geological structures and surface buildings and their equivalent densities, including the following steps S1011 to S1012.

[0023] S1011, Establish an underground geological structure grid based on geological exploration borehole data in the geological exploration data; S1012, establish a surface building grid based on the engineering drawings in the surface building data, and merge the underground geological structure grid and the surface building grid to form a three-dimensional forward modeling grid.

[0024] For example, based on preliminary geological drilling data and engineering drawings, a three-dimensional forward mesh model including geological strata and surface buildings is constructed. Surface buildings are classified by function and structure (such as office buildings, residential buildings, and construction projects), and initial equivalent density estimates are assigned to each type of building. The initial geometric height is derived by combining the number of floors recorded in the engineering drawings with building codes.

[0025] S102, during the tunnel boring machine excavation process, natural muon detection data is continuously collected, and the natural muon detection data collected in a single segment assembly cycle is mapped to the spatial position when the tunnel boring machine is stationary within the cycle to obtain the measured muon survival rate at the spatial position. Specifically, in step S102, during the tunnel boring machine excavation process, natural muon detection data is continuously collected, and the natural muon detection data collected within a single segment assembly cycle is mapped to the spatial position when the tunnel boring machine is stationary within the cycle to obtain the measured muon survival rate at the spatial position, including the following steps S1021 to S1024.

[0026] S1021, A first muon detector is installed behind the cutterhead of the tunnel boring machine, and a second muon detector is installed on the supporting trolley behind the tunnel boring machine. S1022, natural muon detection data are continuously collected during the shield tunneling process through the first muon detector and the second muon detector; S1023, the natural muon detection data collected during a single segment assembly cycle is equivalently mapped to the spatial position of the tunnel boring machine when the segment assembly is in a static state during the cycle, and the mapped data is used as the static observation data at the spatial position. S1024. Calculate the measured muon flux based on static observation data, obtain the measured muon flux under unobstructed airborne conditions, normalize the measured muon flux based on the airborne muon flux, and obtain the measured muon survival rate.

[0027] For example, muon detectors (including a main detector mounted behind the tunnel boring machine cutterhead and an auxiliary detector on a matching trolley) are deployed to continuously acquire raw muon tracking data during tunnel boring. For the alternating cycles of tunnel boring and segment assembly, a dynamic coordinate system is established. Data from a single or multiple excavation cycles (approximately 3-4 hours) is projected and equivalently represented to spatial nodes at the stationary position of the segment assembly, achieving quasi-static data fusion under dynamic displacement to obtain the measured muon survival rate. .

[0028] S103, extract the directional intervals unobstructed by surface buildings from the 3D forward mesh model, calculate the forward muon survival rate within the directional intervals, and establish a correlation model between muon penetration energy and equivalent mass thickness based on the forward muon survival rate and the measured muon survival rate; here, the correlation model is a linear model, expressed as: ; in, Indicates the first fitting coefficient; Indicates the second fitting coefficient; Indicates equivalent mass thickness; This represents the minimum penetrating energy of a muon.

[0029] Specifically, in step S103, the directional intervals that are not obstructed by surface buildings are extracted from the three-dimensional forward mesh model, and the survival rate of forward muons is calculated within the directional intervals, including the following steps S1031 to S1033.

[0030] S1031, perform ray tracing on the 3D forward mesh model, traverse the intervals in each direction, and determine whether there are mesh cells of surface buildings on each ray path; S1032 marks the directional intervals where there are no surface building grid cells on the ray path as unobstructed directional intervals; S1033: Extract the formation density and thickness data in each direction within the unobstructed direction interval, and calculate the theoretical flux decay of muons after they pass through the formation through forward modeling to obtain the forward modeling muon survival rate in the unobstructed direction interval.

[0031] Specifically, in step S103, a correlation model between muon penetration energy and equivalent mass thickness is established based on the forward muon survival rate and the measured muon survival rate, including: (1) Based on the forward muon survival rate and the measured muon survival rate in each direction within the unobstructed direction interval, calculate the minimum penetration energy of the measured muon in each direction, and calculate the equivalent mass thickness in each direction based on the three-dimensional forward mesh model. (2) Using the equivalent mass thickness as the independent variable and the measured minimum penetration energy of the muon as the dependent variable, the data in each direction within the unobstructed direction range are fitted to establish a correlation model between the muon penetration energy and the equivalent mass thickness.

[0032] For example, extract the azimuth data of the pure stratigraphic space unobstructed by buildings, and calculate the forward-modeled survival rate and the measured survival rate. Establish the minimum penetration energy of muons. Equivalent length (mass thickness) Linear approximate physical model ( By using multi-parameter fitting, the baseline coefficient of energy loss rate under the current working condition is obtained, thereby eliminating the background error caused by fluctuations in the energy loss rate of formation materials.

[0033] S104, Based on the correlation model, decouple and calibrate the equivalent density and geometric height of surface buildings in the three-dimensional forward mesh model; Specifically, in step S104, the equivalent density and geometric height of surface buildings in the three-dimensional forward mesh model are decoupled and calibrated according to the association model, including the following steps S1041 to S1043.

[0034] S1041, select the zenith angle direction interval through which the ray path passes through the surface building, and within the azimuth angle interval blocked by the surface building, take the equivalent density of the surface building as the variable, and take minimizing the residual between the measured muon survival rate and the forward model muon survival rate as the objective, and solve the equivalent density through iterative inversion. S1042, based on the convergence of equivalent density, obtain the first decay inflection point of the measured muon survival rate as a function of zenith angle, and the second decay inflection point of the forward modeling muon survival rate as a function of zenith angle, and adjust the geometric height of the surface buildings until the first decay inflection point and the second decay inflection point coincide on the angle axis. S1043, repeat the iterative inversion solution of the equivalent density and the iterative adjustment of the geometric height until both the equivalent density and the geometric height meet the preset convergence conditions.

[0035] For example, equivalent density isolated calibration: to avoid interference from height errors, extract the large zenith angle (e.g. The data includes observations of rad (radius) and the rays that necessarily penetrate the main body of the building. Within the azimuth interval where there is obstruction, a least-squares residual objective function is constructed based on the measured survival rate and the forward modeling theory survival rate. The residual is minimized through iterative optimization, and the global optimal equivalent density of various types of buildings is solved independently.

[0036] Geometric height synchronization calibration: Based on density convergence, the measured and forward-modeled survival rate curves with zenith angle are used. The model height is used to calibrate the synchronous characteristic of the inflection point of muon survival rate decay, which is generated by incremental increases. The abrupt decay phase at which measured data penetrates the actual roofline is identified. If the forward decay phase lags, the model building height is iteratively increased until the decay inflection point of both curves is reached. The axes are completely aligned. By combining continuous dynamic sweep data from multiple shield tunneling rings, the geometric height of the main building complexes along the route is reconstructed in three dimensions.

[0037] S105 updates the calibrated equivalent density and geometric height to the three-dimensional forward modeling grid, filters the natural muon detection data, performs three-dimensional inversion of geological anomalies, and outputs a subsurface density distribution image.

[0038] Specifically, in step S105, the calibrated equivalent density and geometric height are updated to the 3D forward mesh model, including: (1) Replace the initial equivalent density of the corresponding surface building in the three-dimensional forward mesh model with the calibrated equivalent density; (2) Replace the initial geometric height of the corresponding surface building in the three-dimensional forward mesh model with the calibrated geometric height.

[0039] Specifically, in step S105, after screening the natural muon detection data, three-dimensional inversion of geological anomalies is performed to output a subsurface density distribution image, including the following steps S1051 to S1054.

[0040] S1051, perform forward modeling simulation based on the updated three-dimensional forward mesh model to obtain the forward muon survival rate in each directional interval; S1052, calculate the deviation between the measured muon survival rate and the forward muon survival rate in each directional interval, compare the deviation with a preset deviation threshold, and remove the natural muon detection data corresponding to the directional interval with a deviation greater than the preset deviation threshold. S1053, Based on the three-dimensional forward mesh model, determine whether the ray path of each directional interval is blocked by the ground building. For the directional interval that is not blocked, retain the corresponding natural muon detection data. For the directional interval that is blocked, count the muon count in the directional interval. When the muon count is less than the preset count threshold, remove the corresponding natural muon detection data. S1054 utilizes the remaining natural muon detection data after screening to perform three-dimensional inversion of geological anomalies and outputs an image of underground density distribution.

[0041] For example, the corrected and accurate building parameters are updated in the 3D forward model. In the final 3D inversion of geological anomalies, adaptive ray filtering is implemented: high signal-to-noise ratio intervals not obscured by buildings are extracted first; for obscured areas, a confidence threshold for the minimum count value of a single interval (e.g., bin count) is set. It also forcibly truncates distorted rays where the deviation between the measured and forward models still exceeds the preset tolerance. This completely eliminates false positive artifacts caused by residual errors in the engineering drawings, outputting a high-precision underground three-dimensional density structure.

[0042] In a specific embodiment provided by the present invention, the present invention is applicable to Figure 2 The shield tunneling environment shown below is used to clearly and completely describe the technical solution of this invention in order to facilitate understanding and demonstrate the feasibility and practicality of this technology.

[0043] To further verify the effectiveness of the method of the present invention, an embodiment of long-term dynamic monitoring applied to complex urban surface environments is provided. In this embodiment, the thickness of the overburden layer of the tunnel is between 9 m and 11.5 m, and the strata mainly include geological media such as fill, cohesive soil, and weathered rock layers. The surface environment in which this embodiment is located includes various typical urban buildings and intricate underground pipe networks, exhibiting typical surface time-varying mass load interference. Based on known geological exploration and prior data on underground pipeline distribution, this embodiment quantitatively compares and calibrates the muon-based dynamic imaging and interference suppression method to verify its reliability in complex engineering scenarios.

[0044] During the system deployment phase, see Figure 3 Two highly integrated muon detection systems were deployed. Detector No. 1 was installed on the main working platform behind the tunnel boring machine (TBM) cutterhead, responsible for detecting the tunnel face and overlying strata. Detector No. 2 was mounted on a supporting trolley at the rear. In practical applications, due to the alternation of continuous displacement and static states within a single segment assembly cycle (approximately 3-4 hours), this embodiment introduces a quasi-static approximate data fusion strategy. Ignoring the 1.5 m spatial displacement zone during excavation, all muon data within a single ring excavation cycle are equivalent to static observation data of the segment assembly positions, thus ensuring sufficient muon data for subsequent calculations in the case of a single ring.

[0045] To address the complex urban surface environment, this study constructed a detailed three-dimensional forward model of stratigraphy and surface structures based on detailed survey data, such as... Figure 4 As shown, the top view and building numbering method are as follows: Figure 7 As shown in Table 1, the model delineates the main geological structures from bottom to top, including completely weathered granite, gravelly cohesive soil, and plain fill. The reference density of the strata is directly taken from the average value of samples taken from previous geological exploration boreholes. For surface buildings, they are classified into three types: office buildings, residential buildings, and buildings under construction. Their initial density is obtained through estimation, and their geometric height is roughly calculated based on the number of floors in the engineering drawings. Detailed data requires inverse iterative calibration and constraints using measured muon data; specific parameters are shown in Table 1.

[0046] Table 1

[0047] In practical applications of natural muon imaging, background flux under unobstructed conditions needs to be obtained through aerial measurements to calculate the measured muon survival rate Rprd. ; In the formula, Itarget is the measured muon flux, and Ifree-sky is the free-sky measured muon flux. The measured muon flux corresponds to the muon energy spectrum obtained by integrating from the minimum energy Emin of the muon in the corresponding direction. The free-sky flux is obtained by integrating from the sum of the energy Edet required to pass through the detector in the corresponding direction and the muon rest energy Eμ.

[0048] Because for muons with energies less than 100 GeV, their energy loss rate in the formation can be approximated as a constant of 2 MeVcm²g. 1. By extracting data from areas without building obstruction, a minimum penetration energy is established. With equivalent length linear physical model (like Figure 6 (As shown), thereby eliminating the background fluctuations caused by the non-uniformity of energy loss rate of formation materials.

[0049] Then, muon data from areas without obstruction by the measured object were selected for physical model fitting. Combined with multi-ring measured data, the following can be obtained: Figure 6 The fitting results shown demonstrate that this physical model can effectively characterize the energy loss of muons in such scenarios. This step is indispensable because accurate physical models are needed to calculate survival rate information during the forward modeling process.

[0050] After obtaining an accurate physical model, surface building parameters need to be corrected. First, the floor density needs to be corrected. The floor density values ​​are then calibrated. To eliminate the influence of height errors, data with a relatively large zenith angle (0.534 rad) is used for analysis. In the direction without buildings, i.e. Figures 5a-5c The φ angle range is approximately -0.6 rad to 0.4 rad, and the measured values ​​agree well with the simulation. However, within the azimuth angle range where buildings are present, the forward model survival rate exhibits extremely high sensitivity to the equivalent building density, compared to... Figure 5a and Figure 5b It can be seen that even if the density changes by 0.1 g / cm³ 3 This can also lead to a significant deviation between the forward model survival rate and the measured value. Subsequently, within a specific azimuth angle interval where building shading exists, a least-squares residual objective function based on the measured survival rate Rexp and the forward model theoretical survival rate Rpre is constructed: ; Where Ω represents the set of azimuth angles obscured by the building. ρ eq Let be the equivalent density of the building to be determined. Iterative adjustments are made by introducing gradient descent or quasi-Newton methods. ρ eq until the objective function F( ρ eq The convergence was brought within the allowable error range. In this embodiment, after calibration of a certain shield tunneling ring position, the average equivalent density of a certain residential building was finally accurately corrected to 0.2 g / cm³. 3 The accurate density information is shown in Table 2.

[0051] Table 2

[0052] After determining the density, the geometric height of the buildings needs to be calibrated, mainly due to the inconsistent standards for building floor heights. The core logic lies in analyzing the synchronicity of the inflection point between forward modeling and measured survival rate decline as the zenith angle θ increases. For example... Figures 8a-8d As shown, if the building height is relatively accurate, then as θ gradually increases from 0, the forward modeling and measured data will decrease synchronously. If the building height is too low, such as... Figures 8e to 8h As shown, the measured data attenuates due to penetration of the actual building rooftops as the θ angle increases, but the forward modeling data exhibits a lag in attenuation; conversely, the forward modeling survival rate decreases earlier. By combining the azimuth sweep characteristics of multi-ring data, this study ultimately determined the geometric height of the main building complexes along the route.

[0053] Calibration results show that this correction strategy can effectively eliminate macroscopic background system errors caused by surface buildings at the current system resolution scale, as shown in Table 2. It can also overcome flux fluctuation interference caused by the inherent local heterogeneity of building internal structures (such as floor slabs, beams, columns, and cavities) to a certain extent. Through effective correction of complex surface backgrounds, this method extracts more accurate muon survival rate changes, providing a high-quality data foundation for screening ray flux changes caused by geological anomalies, and thus providing strong support for accurate early warning of geological hazards during tunnel boring machine (TBM) excavation.

[0054] The various embodiments described in this specification are presented in a progressive manner. Similar or identical parts between embodiments can be referred to interchangeably. Each embodiment focuses on its differences from other embodiments. All or part of this invention can be used in numerous general-purpose or special-purpose computer system environments or configurations. Examples include: personal computers, server computers, handheld or portable devices, tablet devices, mobile communication terminals, multiprocessor systems, microprocessor-based systems, programmable electronic devices, network PCs, minicomputers, mainframe computers, and distributed computing environments including any of the above systems or devices, etc.

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

Claims

1. A method for suppressing the background of dynamic natural muon imaging under complex surface environments, characterized in that, include: Obtain geological exploration data and surface building data of the target area, and construct a three-dimensional forward mesh model including underground geological structure and surface building geometric parameters and equivalent density; During the tunnel boring machine (TBM) excavation process, natural muon detection data is continuously collected. The natural muon detection data collected within a single segment assembly cycle is mapped to the spatial position of the TBM when it is stationary within the cycle, and the measured muon survival rate at the spatial position is obtained. Extract the directional intervals that are not obstructed by the surface buildings from the three-dimensional forward mesh model, calculate the forward muon survival rate within the directional intervals, and establish a correlation model between muon penetration energy and equivalent mass thickness based on the forward muon survival rate and the measured muon survival rate. Based on the association model, the equivalent density and geometric height of surface buildings in the three-dimensional forward mesh model are decoupled and calibrated. The calibrated equivalent density and geometric height are updated to the three-dimensional forward modeling grid. After filtering the natural muon detection data, three-dimensional geological anomaly inversion is performed to output a subsurface density distribution image.

2. The method for suppressing the background of dynamic natural muon imaging under complex surface environments as described in claim 1, characterized in that, The construction includes a three-dimensional forward mesh model comprising underground geological structure, surface building geometry parameters, and equivalent density, including: A subsurface geological structure grid is established based on the geological exploration data from the geological borehole data. A surface building grid is established based on the engineering drawings in the surface building data, and the underground geological structure grid and the surface building grid are merged to form the three-dimensional forward model.

3. The method for suppressing the background of dynamic natural muon imaging under complex surface environments as described in claim 1, characterized in that, During the tunnel boring machine (TBM) excavation process, natural muon detection data is continuously collected. The natural muon detection data collected within a single segment assembly cycle is mapped to the spatial location of the TBM when it is stationary within that cycle, obtaining the measured muon survival rate at that spatial location. This includes: The first muon detector is installed behind the cutterhead of the tunnel boring machine, and the second muon detector is installed on the supporting trolley behind the tunnel boring machine. Natural muon detection data are continuously collected during the shield tunneling process using the first muon detector and the second muon detector. The natural muon detection data collected within a single segment assembly cycle is equivalently mapped to the spatial position of the tunnel boring machine when it is in a static state during segment assembly within the cycle, and the mapped data is used as the static observation data at the spatial position. The measured muon flux is calculated based on the static observation data to obtain the measured muon flux under unobstructed aerial conditions. The measured muon flux is then normalized based on the measured muon flux to obtain the measured muon survival rate.

4. The method for suppressing the background of dynamic natural muon imaging under complex surface environments as described in claim 1, characterized in that, The step of extracting directional intervals from the 3D forward mesh model that are not obstructed by the surface buildings, and calculating the forward muon survival rate within the directional intervals, includes: Ray tracing is performed on the three-dimensional forward mesh model, traversing each directional interval, and determining whether there are mesh cells of the surface buildings on each ray path; The directional intervals in which no surface building grid cells are present on the ray path are marked as unobstructed directional intervals; Extract the formation density and thickness data in each direction within the unobstructed directional interval, and calculate the theoretical flux decay of muons after they pass through the formation using forward modeling to obtain the forward modeling muon survival rate in the unobstructed directional interval.

5. The method for suppressing the background of dynamic natural muon imaging under complex surface environments as described in claim 1, characterized in that, The establishment of a correlation model between muon penetration energy and equivalent mass thickness based on the forward muon survival rate and the measured muon survival rate includes: Based on the forward muon survival rate and the measured muon survival rate in each direction within the unobstructed direction interval, calculate the minimum penetration energy of the measured muon in each direction, and calculate the equivalent mass thickness in each direction based on the three-dimensional forward mesh model. Using the equivalent mass thickness as the independent variable and the measured minimum muon penetration energy as the dependent variable, the data in each direction within the unobstructed direction interval are fitted to establish a correlation model between muon penetration energy and equivalent mass thickness.

6. The method for suppressing the background of dynamic natural muon imaging under complex surface environments as described in claim 1, characterized in that, The step of decoupling and calibrating the equivalent density and geometric height of surface buildings in the three-dimensional forward mesh model according to the association model includes: The zenith angle direction range through which the ray path passes through the surface building is selected. Within the azimuth angle range blocked by the surface building, the equivalent density of the surface building is used as a variable, and the objective is to minimize the residual between the measured muon survival rate and the forward muon survival rate. The equivalent density is solved by iterative inversion. Based on the convergence of the equivalent density, the first decay inflection point of the measured muon survival rate as a function of the zenith angle and the second decay inflection point of the forward modeling muon survival rate as a function of the zenith angle are obtained. The geometric height of the surface buildings is adjusted until the first decay inflection point and the second decay inflection point coincide on the angle axis. Repeat the iterative inversion solution of the equivalent density and the iterative adjustment of the geometric height until both the equivalent density and the geometric height satisfy the preset convergence condition.

7. The method for suppressing the background of dynamic natural muon imaging under complex surface environments as described in claim 1, characterized in that, The step of updating the calibrated equivalent density and geometric height to the 3D forward mesh model includes: Replace the initial equivalent density of the corresponding surface building in the three-dimensional forward mesh model with the calibrated equivalent density; Replace the initial geometric height of the corresponding surface building in the three-dimensional forward mesh model with the calibrated geometric height.

8. The method for suppressing the background of dynamic natural muon imaging under complex surface environments as described in claim 1, characterized in that, The process of filtering the natural muon detection data and performing three-dimensional inversion of geological anomalies to output a subsurface density distribution image includes: Forward simulation was performed based on the updated three-dimensional forward mesh model to obtain the forward muon survival rate in each directional interval; Calculate the deviation between the measured muon survival rate and the forward muon survival rate in each directional interval, compare the deviation with a preset deviation threshold, and discard the natural muon detection data corresponding to the directional intervals where the deviation is greater than the preset deviation threshold; Based on the three-dimensional forward mesh model, it is determined whether the ray path of each directional interval is blocked by the ground building. For the directional interval that is not blocked, the corresponding natural muon detection data is retained. For the directional interval that is blocked, the muon count in the directional interval is counted. When the muon count is less than a preset count threshold, the corresponding natural muon detection data is removed. The remaining natural muon detection data after screening are used to perform three-dimensional inversion of geological anomalies and output an image of underground density distribution.

9. The method for suppressing the background of dynamic natural muon imaging under complex surface environments as described in claim 1, characterized in that, The correlation model is a linear model, represented as follows: ； in, Indicates the first fitting coefficient; Indicates the second fitting coefficient; Indicates equivalent mass thickness; This represents the minimum penetrating energy of a muon.

Images