Water-inflow channel fusion recognition method based on nuclear magnetic resonance and ground-air electromagnetic cooperation

Abstract

This invention discloses a method for identifying water inflow channels based on nuclear magnetic resonance (NMR) and ground-to-air electromagnetic (GTEM) fusion, applicable to the fields of geophysical exploration and hydrogeology. It addresses the problems of multiple solutions and insufficient accuracy when using single geophysical methods to identify water inflow channels. Specifically, this invention targets open-pit coal mine exploration areas, simultaneously implementing GTEM scanning, NMR measurement points, and flow field profiling to acquire multi-source geophysical data. A three-dimensional integrated geological model incorporating resistivity, water content, and hydraulic conductivity parameters is constructed. Based on an inversion framework of coupled electromagnetic induction, NMR relaxation, and seepage mechanisms, a synergistic inversion is performed to generate a multi-parameter fused attribute volume. Furthermore, hydrogeological feature definitions and machine learning spatial pattern recognition are used to automatically delineate suspected water inflow channels. Finally, by optimizing the layout of verification boreholes and dynamically correcting the model and identification patterns using measured data, high-precision, high-efficiency integrated identification and three-dimensional visualization of water inflow channels are achieved.

Description

Technical Field

[0001] This invention belongs to the field of geophysical exploration and hydrogeology, specifically a method for identifying water inrush channels based on nuclear magnetic resonance and ground-to-space electromagnetic fusion. Background Technology

[0002] In the process of safe production and efficient mining in open-pit mines, accurate exploration of hydrogeological conditions and effective prevention and control of water hazards such as water inrush and mudslides are crucial core issues. Water inflow channels, as the preferred paths for groundwater to converge, migrate, and concentrate in the mining pit in complex geological bodies, are directly related to the reliability of drainage design optimization, slope stability assessment, and early warning of major water hazard risks.

[0003] Traditionally, the detection of underground water channels has relied primarily on drilling combined with hydrological tests. While this method is direct and reliable, it suffers from limitations such as "one-hole observation," high costs, low efficiency, and difficulty in constructing continuous and complete three-dimensional spatial distribution images of the channels between boreholes. Geophysical methods, due to their non-destructive, efficient, and wide spatial coverage advantages, have become an important supplement to and even a partial replacement for drilling. Currently, various geophysical methods are applied to hydrogeological surveys, each with its own characteristics: Frequency Domain Electromagnetic Spectrometry (FDEM) can rapidly perform large-area scanning by detecting differences in geoelectric structure and is sensitive to low-resistivity anomalies (such as water-rich zones); Nuclear Magnetic Resonance (NMR) groundwater detection technology is one of the few methods that can directly and quantitatively detect underground water content and pore structure information; while flow field methods (or natural electric field methods / charging methods, etc.) have unique responses to natural or artificial electric field anomalies generated by groundwater seepage movement and can directly indicate the activity of water-conducting channels.

[0004] While each of the aforementioned individual methods can provide valuable information, in the complex geological environment of actual open-pit coal mines, any single geophysical response exhibits multiple interpretations. For example, low resistivity anomalies may be caused by water-rich areas or by areas rich in conductive minerals such as clay; nuclear magnetic resonance signals are greatly affected by electromagnetic noise and are difficult to distinguish between static and dynamic reserves; the spatial positioning accuracy of flow field anomalies is relatively limited. Therefore, in recent years, a clear development trend in the field of geophysical exploration has been towards the synergy and integration of multiple methods and parameters. By integrating various physical field data with different sensitivities to the occurrence and movement of underground fluids, a unified geological-geophysical model can be constructed for joint interpretation and inversion, aiming to mutually constrain and complement each other, thereby reducing multiple interpretations and improving the detection accuracy and interpretation reliability of complex geological targets such as water inflow channels. However, how to achieve deep coupling and intelligent identification of multiple heterogeneous geophysical data at the levels of acquisition, processing, inversion, and even final geological interpretation remains a technical challenge that urgently needs in-depth research and resolution. Against this backdrop, developing a systematic method for three-dimensional fusion identification of water inrush channels that can synergistically utilize multi-source information such as nuclear magnetic resonance, geo-electromagnetic resonance, and flow field methods is of significant theoretical and engineering application value.

[0005] The following problems exist in the existing technology: Existing technologies often rely on a single method (such as using only ground-to-air electromagnetic or nuclear magnetic resonance methods) for detection. Ground-to-air electromagnetic methods are sensitive to low-resistivity bodies, but cannot distinguish whether the low resistance is caused by groundwater or conductive rock layers. Nuclear magnetic resonance methods can directly detect water content, but have limited resolution for static water content and dynamic water flow, and are easily affected by electromagnetic noise. Flow field methods can indicate seepage paths, but have low positioning accuracy and depth resolution. Using any one method alone makes it difficult to uniquely and accurately characterize complex water inflow channels that have both specific water-bearing and water-conducting properties, leading to a high risk of misjudgment or missed detection. Existing multi-method integrated exploration often stops at the stage of simple comparison and manual overlay interpretation of the resulting maps after data acquisition, or only performs loose joint inversion. Data from different geophysical methods are independent of each other during the inversion process, lacking unified modeling and collaborative constraints based on core physical processes such as electromagnetic induction, nuclear magnetic resonance relaxation, and groundwater seepage. This "soft integration" approach cannot fundamentally overcome the multiple solutions and is difficult to generate a three-dimensional geological model that is self-consistent in terms of physical mechanisms and can truly reflect the spatial correlation of multiple parameters such as resistivity, water content, and hydraulic conductivity. Existing technologies for delineating water inflow channels mainly rely on interpreters' visual judgment and empirical correlation of geophysical anomaly features, lacking quantitative and standardized identification criteria. This process is highly subjective, inefficient, and difficult to handle massive three-dimensional data volumes. It is also unable to systematically uncover complex patterns hidden in the spatial combination of multiple physical property parameters, and it is even more difficult to quantitatively evaluate the confidence level of the identification results. In traditional processes, geophysical interpretation results (such as anomalies) are used to guide the deployment of verification boreholes. However, the actual hydrogeological data obtained from the boreholes are usually only used to verify the conclusions and are rarely systematically fed back to the previous geophysical models for calibration and optimization. Once the model is generated, it remains fixed and cannot use new information to reduce the uncertainty of the model itself, nor can it adaptively optimize the identification algorithm parameters, which limits the method's ability to continuously improve accuracy in iterative exploration. Summary of the Invention

[0006] The present invention aims to at least solve one of the technical problems existing in the prior art; to this end, the present invention proposes a method for identifying water inrush channels based on nuclear magnetic resonance and ground-to-air electromagnetic synergy, in order to solve the above-mentioned technical problem.

[0007] The first aspect of the present invention provides a method for fusion identification of water inrush channels based on nuclear magnetic resonance and ground-to-air electromagnetic coordination, comprising the following steps: S1: Within the target exploration area, simultaneously implement frequency domain ground-to-space electromagnetic method planar scanning, nuclear magnetic resonance groundwater detection method point measurement, and flow field method profile measurement to obtain multi-source geophysical field data; S2: Construct a three-dimensional integrated initial geological model that includes resistivity, water content and hydraulic conductivity parameters; using the multi-source geophysical field data as joint constraints, iteratively optimize the three-dimensional integrated initial geological model through an inversion framework that couples electromagnetic induction, nuclear magnetic resonance relaxation and groundwater seepage physical mechanisms to generate a three-dimensional multi-parameter fused attribute body. S3: Based on hydrogeological principles, define the multi-dimensional feature vector combination pattern corresponding to the water inflow channel in the three-dimensional multi-parameter fusion attribute volume; use machine learning algorithms to search and classify spatial patterns in the three-dimensional multi-parameter fusion attribute volume, automatically identify and delineate abnormal spatial areas that conform to the multi-dimensional feature vector combination pattern, mark them as suspected water inflow channels, and output a three-dimensional map of the identification results with channel spatial morphology, category and confidence information. S4: Deploy verification boreholes at key locations of suspected water inflow channels indicated by the three-dimensional map of the identification results; use the measured hydrogeological data revealed by the boreholes to correct and optimize the three-dimensional multi-parameter fusion attribute volume and the feature vector combination mode.

[0008] Preferably, in step S1, within the target exploration area, frequency-domain ground-to-air electromagnetic method planar scanning, nuclear magnetic resonance groundwater detection method point measurement, and flow field method profile measurement are simultaneously implemented, including the following steps: The target exploration area of ​​open-pit coal mines is divided based on a unified two-dimensional basic grid. Based on the preset geological structure information and hydrological model, a survey line network for frequency domain ground-to-air electromagnetic method surface scanning is automatically generated on the two-dimensional basic grid. In particular, in known anomalous areas and anomalous areas inferred from prior geological information, encrypted verification survey lines are generated. The location of groundwater detection points is dynamically determined by optimizing a function, which is: ,in, Optimized set of measurement point locations, It is a collection of known geological structural location information. This is a set of potential anomaly regions predicted by frequency-domain ground-to-space electromagnetic methods based on geological structural information and hydrological models. This is the set of all possible measurement point locations within the target exploration area. For spatial density function, This is a spatial overlap function. For spatial coverage function, , , These are the weighting coefficients; based on the optimization function. From the optimization results, the measuring points with the highest overlap with geological structures and predicted anomalies were selected as key verification measuring points. Lay out the survey lines for flow field method profile measurement, ensuring that the survey lines pass through key verification points and the areas where the denser verification survey lines are located.

[0009] Preferably, the step S1, acquiring multi-source geophysical field data, includes the following steps: For the frequency domain ground-to-air electromagnetic method, a ground-based transmitting device transmits specifically coded electromagnetic waves into the ground; an aircraft carrying an electromagnetic receiving system flies along the designed survey line to collect the electromagnetic response signals of the underground medium; the transmitted electromagnetic waves adaptively select the fundamental frequency pseudo-random waveform coding group according to the environmental noise level; and the flight altitude of the aircraft is dynamically adjusted according to real-time terrain data. For nuclear magnetic resonance (NMR) groundwater detection, a ground-based transmitter and electromagnetic receiver system are used to acquire NMR signals and the resulting attenuation curves. The signal-to-noise ratio of the NMR signals is calculated in real time during the acquisition process. The decay curve was initially fitted using a single exponential method to obtain the correlation coefficient. Set the signal-to-noise ratio threshold. Correlation coefficient threshold ,achieve and If either of the two conditions is met, the operation of increasing the number of superpositions and switching the pulse moment sequence will be automatically triggered; if the threshold condition is still not met after the operation is triggered, the measurement point will be marked as an interference point. For the flow field method, before each profile measurement, the background potential gradient distribution is collected in an area assuming no concentrated leakage. After formal measurements, the abnormal potential gradient caused by the seepage channels was initially separated. ,in This represents the measured total potential gradient; The raw data collected by the three methods are all assigned a unified spatiotemporal reference coordinate. The observed values, the quality indicators generated based on the data quality judgment results, and the corresponding acquisition parameters are combined to package and generate a standard data package to form multi-source geophysical field data.

[0010] Preferably, in step S2, a three-dimensional integrated initial geological model is constructed, including resistivity, water content, and hydraulic conductivity parameters. Using the multi-source geophysical field data as joint constraints, the three-dimensional integrated initial geological model is iteratively optimized through an inversion framework that couples electromagnetic induction, nuclear magnetic resonance relaxation, and groundwater seepage physical mechanisms. This includes the following steps: Based on the multi-source geophysical field data and known hydrogeological information, a three-dimensional integrated initial geological model is constructed, and the target exploration area is discretized into... The model comprises three-dimensional mesh cells, each containing resistivity. Volumetric water content θ, longitudinal relaxation time The four physical properties are: hydroconductivity K; Within the inversion framework, a multiphysics collaborative inversion objective function is established, which is: Where m is the parameter vector of the three-dimensional integrated initial geological model. , , These are the data fitting differences for the frequency-domain ground-to-space electromagnetic method, the nuclear magnetic resonance groundwater detection method, and the flow field method, respectively. For model parameter smoothing constraints, For cross-gradient structure coupling constraint terms, and For regularization parameters; Among them, the cross-gradient structure coupling constraint term The calculation formula is: in, , , Let i and n represent the spatial gradient vectors of resistivity, water content, and hydraulic conductivity in the i-th 3D mesh cell, respectively, i=1,2,... , × represents the cross product of vectors; Using a three-dimensional integrated initial geological model as the starting point for iteration and multi-source geophysical field data as joint constraints, the objective function is minimized through optimization algorithms within the inversion framework. The model parameters are iteratively updated.

[0011] Preferably, in step S2, generating a three-dimensional multi-parameter fused attribute body includes the following steps: During the iterative process of the inversion framework, the physical property relationship between resistivity, water content, and relaxation time is established through a rock physics bridging equation, which is as follows: ;in, The predicted conductivity value, Porosity estimated from water content θ and lithology For saturation, Let A be the electrical conductivity of pore water, and A be an empirical coefficient. and For cementation index and saturation index; The data fitting term of the flow field method Forward modeling is performed based on the coupled equations of the seepage field and the current field. The coupled equations are as follows: ;in, The distribution of groundwater head. For potential distribution, For the injected current intensity, Let r be the location of the current source point, and r be the spatial position vector. For the Dirac function, For gradient operators, and The spatial distributions of hydraulic conductivity and electrical conductivity, defined by the model parameter vector m, are respectively used. By solving the coupling equations, the predicted potential gradient considering the groundwater flow effect is obtained, and then the fitting difference with the measured data is calculated. The iterative process of its inversion framework employs a constrained optimization algorithm, updating the model parameter vector m in each iteration until the objective function is reached. The change in resistivity satisfies one of the following two conditions: less than a preset threshold or reaching the maximum number of iterations. The output is a three-dimensional multi-parameter fused attribute body containing the spatial distribution of resistivity, water content, relaxation time, and hydroconductivity.

[0012] Preferably, in step S3, based on hydrogeological principles, the multi-dimensional feature vector combination pattern corresponding to the water inflow channel in the three-dimensional multi-parameter fused attribute volume is defined, including the following steps: The basic physical property parameters of each three-dimensional mesh unit i are extracted from the three-dimensional multi-parameter fused attribute volume, and the derived features of each three-dimensional mesh unit i are calculated based on the basic physical property parameters to construct a seven-dimensional basic feature vector. Among them, resistivity Volumetric water content Longitudinal relaxation time Hydraulic conductivity , The coefficient of conductivity Spatial gradient magnitude, The water conductivity index, The risk index for water inrush. Among them, the water conductivity index The calculation formula is: ,in, The hydraulic gradient at point i in the three-dimensional mesh. Water inrush risk index The calculation formula is: ,in, It is a very small constant; Based on hydrogeological principles, characteristic fingerprint patterns of several typical geological targets are predefined, including fingerprints of active water-conducting channels, fingerprints of water-rich aquitards, fingerprints of weakly water-bearing fractures, and fingerprints of dry and intact rock strata.

[0013] Preferably, in step S3, a machine learning algorithm is used to search and classify spatial patterns in a three-dimensional multi-parameter fusion attribute volume, automatically identifying and delineating abnormal spatial regions that conform to the multi-dimensional feature vector combination pattern, marking them as suspected water inrush channels, and outputting a three-dimensional image of the identification results with channel spatial morphology, category, and confidence information, including the following steps: The fundamental feature vectors of all three-dimensional mesh elements The vectors in the seven-dimensional feature space are standardized and projected into a seven-dimensional feature space. A density clustering algorithm is used to perform cluster analysis on the vectors in the seven-dimensional feature space. Each three-dimensional grid cell is assigned a cluster label, and cells without cluster labels are marked as noise points. For each non-noise cluster Calculate its cluster average eigenvector. and calculate Based on the similarity distance between each feature fingerprint pattern, each cluster is initially classified into one of the following categories according to the nearest neighbor principle: suspected water inflow channel, water-rich layer, weakly water-bearing fracture, or dry rock layer. For clusters initially classified as suspected water inflow channels, after three-dimensional spatial morphological filtering, those with a volume greater than a preset threshold are selected based on three-dimensional neighborhood connectivity. Furthermore, the spatial distribution can form a connected domain with a continuous path from the potential replenishment boundary to the discharge area, and then the spatial attitude parameters of the connected domain can be extracted, including the principal direction, dip angle and thickness. For each finally determined suspected water inflow channel connectivity region Calculate the overall confidence score The calculation formula is as follows: ;in, For connected components Average eigenvector fingerprints with active water-conducting channels similarity, For connected components volume, The total volume of all suspected water inflow channels connected domains. For connectivity score, , , Preset weighting coefficients; The final result is a 3D image of the recognition outcome, which includes information on the spatial morphology, category, and confidence level of the channels.

[0014] Preferably, step S4 includes the following steps: Based on the identified results, the connected domains corresponding to each suspected water inflow channel in the 3D image are... Confidence score And the posterior uncertainty estimation of each parameter in the three-dimensional multi-parameter fused attribute volume M; the optimal location of the verification borehole is determined through a multi-objective optimization function. : Where J(·) is the indicator function, when the optimal borehole position Located in the connected domain The value is 1 if the time condition is met, and 0 otherwise. , , For preset weighting coefficients, For parameter p at position Posterior uncertainty estimation at the location, The three-dimensional spatial boundary of the target exploration area. To characterize from the boundary The inferred seepage path to the mine pit was determined by the borehole location. Path connectivity function for verification degree; In the optimal position Drilling was carried out at the site to obtain measured hydrogeological data. This includes lithological columnar sections, stratified water inflow Q, hydrostatic pressure P, and hydraulic conductivity obtained through in-situ hydrogeological tests. Using measured hydrogeological data The three-dimensional multi-parameter fused attribute volume M is corrected using a Bayesian inversion framework to generate an updated three-dimensional multi-parameter fused attribute volume. Its Bayesian inversion framework constructs borehole data using a forward mapping function. The likelihood function of the 3D multi-parameter fused attribute volume M Using the three-dimensional multi-parameter fused attribute volume M as the prior model, the posterior model distribution is obtained according to Bayes' theorem. , For all measured hydrogeological data The dataset; by sampling the posterior model distribution, an updated model set is obtained, and its mean is calculated as... ; where, the forward mapping function The three-dimensional multi-parameter fused attribute volume M is located at... The parameters at a given location are mapped to predicted values ​​of the z-th type of observation data. For stratified inflow rate Q and hydrostatic pressure P, and The results were obtained through local groundwater numerical simulation and solving the seepage field equations, respectively. Based on the borehole verification results, the multi-dimensional feature vector combination pattern is adaptively optimized, and the corresponding areas are selected according to the actual geological conditions revealed by the borehole. The corresponding grid feature vector The samples are classified and categorized into the confirmed channel sample library and other corresponding predefined geological category sample libraries. Using the updated sample library, the feature fingerprint pattern definition in the multidimensional feature vector combination pattern is updated by updating the corresponding model parameters.

[0015] Compared with the prior art, the beneficial effects of the present invention are: This invention simultaneously implements ground-to-air electromagnetic, nuclear magnetic resonance, and flow field method measurements, and constructs an inversion framework that couples multiple physical mechanisms. Parameters such as resistivity, water content, relaxation time, and hydraulic conductivity are synergistically optimized in a unified three-dimensional model. This "hard fusion" ensures that the inversion results must simultaneously satisfy the responses of multiple physical fields. By utilizing the different sensitivities of different methods to water body endowment and transport information to mutually constrain each other, it uniquely and accurately characterizes the entity of a water-rich and highly conductive inflow channel, greatly improving the reliability of the identification. This invention provides a rich data foundation for channel identification by generating a three-dimensional multi-parameter fusion attribute body containing multiple types of physical properties. Furthermore, by predefining multi-dimensional feature vector combination patterns based on hydrogeological principles and using machine learning algorithms (such as density clustering) for spatial pattern search and classification, it achieves automatic, objective, and quantitative delineation of abnormal areas that conform to the characteristics of water inflow channels. This not only significantly improves identification efficiency but also provides a quantitative basis for the reliability of the results by calculating confidence scores. This invention incorporates the borehole verification process into the methodology and designs a verification borehole location method based on multi-objective optimization. Simultaneously, it utilizes measured data revealed by boreholes to correct the three-dimensional multi-parameter fusion attribute volume using a Bayesian inversion framework and adaptively optimizes the definition of feature vector combination patterns. This closed-loop process enables the entire system to learn and evolve, continuously iterating and updating the model and recognition algorithm as exploration progresses, thereby continuously improving the accuracy of subsequent predictions and achieving self-improvement of the method and dynamic adaptability for engineering applications. The final output of this invention is a three-dimensional image of the identification results with channel spatial morphology, category and confidence information. This result not only intuitively shows the three-dimensional spatial location, geometric shape and distribution path of the water inflow channel, but also includes quantitative attribute and confidence information, which can be directly used to guide the design of prevention and control projects (such as the precise layout of drainage holes and curtain grouting target areas), transforming complex geophysical exploration results into high-value information products that are easy to understand and facilitate engineering decision-making. Attached Figure Description

[0016] Figure 1 This is a schematic diagram of the method flow of the present invention. Detailed Implementation

[0017] The technical solution of the present invention will be clearly and completely described below with reference to the embodiments. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments.

[0018] Please see Figure 1 This invention is a method for identifying water inrush channels based on nuclear magnetic resonance and ground-to-air electromagnetic fusion, comprising the following steps: S1: Within the target exploration area, simultaneously implement frequency domain ground-to-space electromagnetic method planar scanning, nuclear magnetic resonance groundwater detection method point measurement, and flow field method profile measurement to obtain multi-source geophysical field data; S2: Construct a three-dimensional integrated initial geological model that includes resistivity, water content and hydraulic conductivity parameters; using the multi-source geophysical field data as joint constraints, iteratively optimize the three-dimensional integrated initial geological model through an inversion framework that couples electromagnetic induction, nuclear magnetic resonance relaxation and groundwater seepage physical mechanisms to generate a three-dimensional multi-parameter fused attribute body. S3: Based on hydrogeological principles, define the multi-dimensional feature vector combination pattern corresponding to the water inflow channel in the three-dimensional multi-parameter fusion attribute volume; use machine learning algorithms to search and classify spatial patterns in the three-dimensional multi-parameter fusion attribute volume, automatically identify and delineate abnormal spatial areas that conform to the multi-dimensional feature vector combination pattern, mark them as suspected water inflow channels, and output a three-dimensional map of the identification results with channel spatial morphology, category and confidence information. S4: Deploy verification boreholes at key locations of suspected water inflow channels indicated by the three-dimensional map of the identification results; use the measured hydrogeological data revealed by the boreholes to correct and optimize the three-dimensional multi-parameter fusion attribute volume and the feature vector combination mode.

[0019] Specifically, firstly, within the target exploration area of ​​the open-pit coal mine, simultaneous surface scanning using frequency-domain ground-to-space electromagnetic methods, point measurements using nuclear magnetic resonance methods, and profile measurements using flow field methods are conducted. These three methods acquire multi-source geophysical field data reflecting underground resistivity structure, water-bearing porosity characteristics, and groundwater seepage activity, laying a unified data foundation for subsequent fusion analysis.

[0020] Next, an initial three-dimensional geological model is constructed, which discretizes the exploration area into grid cells and defines four key parameters—resistivity, water content, relaxation time, and hydraulic conductivity—within each cell. Starting with this model, an inversion framework coupling electromagnetic induction, nuclear magnetic resonance relaxation, and the physical mechanism of groundwater seepage is established. Within this framework, the aforementioned multi-source data are used as joint constraints. By minimizing an objective function that includes the fitting differences between the data from each method and coupling constraints of the model structure, the initial model is iteratively optimized, ultimately generating a three-dimensional multi-parameter fused attribute volume that can consistently reflect the spatial distribution of multiple parameters.

[0021] Then, based on hydrogeological principles, multidimensional feature vectors for each grid cell are extracted and calculated from the three-dimensional attribute volume, thereby predefining the combination pattern of physical properties of the water inflow channels. Machine learning algorithms (such as density clustering) are used to automatically search and classify the entire three-dimensional space, identify and delineate abnormal connected regions that conform to the preset channel pattern, and calculate their spatial morphology, category attribution, and identification confidence. Finally, a visualized three-dimensional identification result map is output.

[0022] Finally, based on the high-confidence areas in the identification results map and the uncertainty distribution of the model, the optimal layout location of the verification boreholes is determined through an optimization algorithm. After drilling, the acquired measured hydrogeological data (such as lithology, water inflow, and measured hydraulic conductivity) are used to correct and optimize the aforementioned three-dimensional multi-parameter fusion attribute volume through a Bayesian update framework. The feature recognition mode is then adaptively updated using new data, thus forming a closed-loop technology of "acquisition-inversion-identification-verification-optimization" that continuously improves the accuracy of the entire method.

[0023] In one embodiment of the present invention, in step S1, within the target exploration area, frequency-domain ground-to-air electromagnetic surface scanning, nuclear magnetic resonance groundwater detection point measurement, and flow field profiling are simultaneously performed, including the following steps: The target exploration area of ​​open-pit coal mines is divided based on a unified two-dimensional basic grid. Based on the preset geological structure information and hydrological model, a survey line network for frequency domain ground-to-air electromagnetic method surface scanning is automatically generated on the two-dimensional basic grid. In particular, in known anomalous areas and anomalous areas inferred from prior geological information, encrypted verification survey lines are generated. The location of groundwater detection points is dynamically determined by optimizing a function, which is: ,in, Optimized set of measurement point locations, It is a collection of known geological structural location information. This is a set of potential anomaly regions predicted by frequency-domain ground-to-space electromagnetic methods based on geological structural information and hydrological models. This is the set of all possible measurement point locations within the target exploration area. For spatial density function, This is a spatial overlap function. For spatial coverage function, , , These are the weighting coefficients; based on the optimization function. From the optimization results, the measuring points with the highest overlap with geological structures and predicted anomalies were selected as key verification measuring points. Lay out the survey lines for flow field method profile measurement, ensuring that the survey lines pass through key verification points and the areas where the denser verification survey lines are located.

[0024] Specifically, based on a pre-defined, unified two-dimensional basic grid, the target exploration area of ​​an open-pit coal mine is systematically spatially discretized. This grid divides the entire planar exploration area into a series of regularly arranged cells. The determination of the cell size (i.e., grid spacing) of the two-dimensional basic grid requires comprehensive consideration of various factors, including the minimum spatial scale of the water inflow channels or other geological targets to be detected, the theoretical detection resolution of the frequency-domain ground-to-space electromagnetic method and nuclear magnetic resonance system used, and the computational requirements of the subsequent three-dimensional inversion model. In practical applications, for exploration at the mine scale (e.g., several square kilometers), the grid cell side length is typically set between 5 and 25 meters; for example, a typical value could be a 10-meter × 10-meter square grid. The coordinate system of this grid adopts the national geodetic coordinate system (e.g., CGCS2000) and elevation datum.

[0025] Based on the established two-dimensional basic grid, the flight measurement path (survey line network) for the frequency domain ground-to-air electromagnetic method (FDEM) is automatically planned. Data acquisition for the FDEM is achieved through an unmanned aerial vehicle (UAV) platform, specifically a multi-rotor UAV or a fixed-wing UAV equipped with an electromagnetic receiving system. In actual operations, the choice of UAV type is determined based on the area of ​​the open-pit coal mine exploration area, the degree of terrain undulation, and the required flight endurance. For larger areas (e.g., exceeding 2 square kilometers) with relatively flat terrain, fixed-wing UAVs can be used due to their strong endurance and high coverage efficiency. For areas with complex terrain, numerous obstacles, or requiring focused densification, multi-rotor UAVs are selected due to their greater maneuverability and ability to hover or fly at varying altitudes. The onboard electromagnetic receiving system typically includes a three-component magnetic induction sensor, a data acquisition unit, a Global Navigation Satellite System (GNSS) receiver, and an inertial measurement unit (IMU) to ensure obstacle avoidance and survey line tracking, while simultaneously recording flight attitude and position information in real time. The ground-based transmitting device is positioned at a suitable location outside or inside the open-pit coal mine exploration area to transmit specially coded electromagnetic signals, maintaining time synchronization with the aircraft's receiving system.

[0026] The generation of its survey line network is based on the intelligent design of pre-set geological structural information (such as known fault lines and stratigraphic boundaries) and hydrological models (such as inferred aquifer distribution and recharge boundaries). The geological structural information comes from existing geological exploration results in the exploration area, including but not limited to regional geological maps, borehole core logging, seismic exploration interpretation profiles, and existing geophysical interpretation maps. This information provides the spatial distribution of structural elements such as faults, folds, stratigraphic interfaces, and lithological boundaries. The hydrological model is constructed based on regional hydrogeological surveys, borehole pumping test data, long-term groundwater level observation data, and topographic analysis. It describes the spatial distribution of aquifers and impermeable layers, groundwater recharge and discharge boundaries, initial flow fields, and hydraulic parameter zoning. When designing the survey line network, the hydrological model is used to preliminarily predict the location of groundwater-rich areas or potential runoff channels (i.e., inferred anomaly areas), thereby guiding the layout of densified survey lines. First, based on the boundaries and main structural trends of the exploration area, a set of basic parallel survey lines covering the entire area is generated. The spacing between the survey lines is typically 2 to 4 times the side length of the basic grid cell (e.g., 20 to 40 meters). To improve the detection accuracy of key areas, densified verification survey lines are automatically added in known and inferred anomaly areas (i.e., reducing the spacing between survey lines or adding intersecting survey lines perpendicular to the main survey line direction). Known anomaly areas include locations where historical water inrushes have occurred, anomaly areas indicated by existing geophysical data, or fracture zones marked on geological maps. Inferred anomaly areas are potential water-conducting channels or water-rich areas inferred from prior geological information and hydrological models through numerical simulations or expert experience. The spacing of the densified survey lines can be reduced to half or less of the spacing between the basic survey lines (e.g., 10 meters), forming a denser local observation network to obtain data with higher spatial resolution.

[0027] Nuclear magnetic resonance (NMR) measurements are point-based measurements, which are relatively expensive and have limited detection range at a single point. Therefore, the placement of these points requires high optimization. A multi-objective optimization function is used to dynamically determine the optimal set of measurement points. The optimization function is defined as: .in, This is the desired, optimized set of nuclear magnetic resonance measurement point locations. It is a collection of known geological structure location information; the information comes from geological mapping, borehole data, seismic exploration interpretation results, etc., and is usually represented as points (such as borehole locations), lines (such as fault lines) or polygonal regions (such as known aquifer ranges); in the algorithm, this information is transformed and mapped onto a two-dimensional base grid. This is a set of potential anomaly regions predicted by frequency domain ground-to-space electromagnetic methods based on geological structural information and hydrological models, through preliminary forward modeling or empirical rules; for example, at the inferred fault location, because its resistivity may be different from the surrounding rock, it can be predicted as an FDEM anomaly region. This refers to the set of all technically feasible and permissible locations for nuclear magnetic resonance (NMR) measurement devices within the target exploration area. Let the set of measuring points P be relative to the set of known geological structures. The spatial density function, which is used to evaluate the spatial density function of NMR points relative to known geological structures. The coverage density; a specific calculation method is to calculate the reciprocal of the average nearest distance from the candidate measurement point set to the known geological structure point / line / surface, or to calculate the proportion of measurement points falling within a certain buffer zone of the known structure; the higher the function value, the better the "proximity" of the measurement point to the key geological structure. Let P be the set of measurement points and the set of predicted anomaly regions. A spatial overlap function is used to evaluate the spatial overlap between NMR measurement points and predicted FDEM anomaly regions. The degree of overlap; specifically, it can be calculated as the proportion of the number of measurement points falling within the grid of the predicted anomaly area in the candidate measurement point set to the total number of measurement points; the higher the function value, the greater the possibility that NMR measurements can directly verify FDEM anomalies. Let P be the spatial coverage function of the set of measurement points P within the target exploration area. This function is used to ensure that the measurement points are within the global exploration area. Uniformity and representativeness of the internal distribution; for example, the Thiessen polygon method can be used to calculate the uniformity index of the coverage of all feasible measurement points by the candidate measurement point set, or to calculate the entropy value of the candidate measurement point set in space; the higher the function value, the more uniform the spatial distribution of the measurement points, and the stronger the exploration of unknown areas. , , These are the weighting coefficients for the three functions, used to adjust the relative importance of different optimization objectives in the overall objective function. Their values ​​depend on the specific exploration objectives and prior knowledge; if the main purpose of exploration is to verify and thoroughly investigate known structures, then... A larger value can be set (e.g., 0.5). and Smaller; if the focus is on detecting unknown FDEM anomalies, then A larger value should be set (e.g., 0.4); if it is necessary to take into account the regional census, then... A certain weight (e.g., 0.1) needs to be assigned; typically, the weight coefficients need to meet the normalization condition: + + =1. Before project implementation, an initial set of weights can be determined through trial calculations or experience, or dynamically adjusted based on feedback from preliminary data collection.

[0028] The optimization function is solved using numerical optimization algorithms (such as genetic algorithms, simulated annealing algorithms, or greedy algorithms) to obtain a set of optimized NMR measurement point locations. Subsequently, key verification measurement points were further selected from the optimization results; the selection criteria were: those points that simultaneously correspond to known geological structures. and predicted FDEM anomaly areas The measurement point with the highest overlap. Specifically, the overlap between each optimized measurement point and... and Spatial correlation (such as inverse distance or membership degree) is used to select the top N (e.g., top 20%) measurement points with the highest comprehensive correlation as key verification measurement points.

[0029] Finally, the flow field method (usually referring to the natural electric field method or charging method for measuring seepage field) profile measurement lines are laid out; the layout of the flow field method measurement lines must ensure that they pass through the areas where the key verification measurement points of nuclear magnetic resonance are located and the areas covered by the densified verification measurement lines of the frequency domain ground-to-space electromagnetic method.

[0030] In one embodiment of the present invention, step S1, acquiring multi-source geophysical field data, includes the following steps: For the frequency domain ground-to-air electromagnetic method, a ground-based transmitting device transmits specifically coded electromagnetic waves into the ground; an aircraft carrying an electromagnetic receiving system flies along the designed survey line to collect the electromagnetic response signals of the underground medium; the transmitted electromagnetic waves adaptively select the fundamental frequency pseudo-random waveform coding group according to the environmental noise level; and the flight altitude of the aircraft is dynamically adjusted according to real-time terrain data. For nuclear magnetic resonance (NMR) groundwater detection, a ground-based transmitter and electromagnetic receiver system are used to acquire NMR signals and the resulting attenuation curves. The signal-to-noise ratio of the NMR signals is calculated in real time during the acquisition process. The decay curve was initially fitted using a single exponential method to obtain the correlation coefficient. Set the signal-to-noise ratio threshold. Correlation coefficient threshold ,achieve and If either of the two conditions is met, the operation of increasing the number of superpositions and switching the pulse moment sequence will be automatically triggered; if the threshold condition is still not met after the operation is triggered, the measurement point will be marked as an interference point. For the flow field method, before each profile measurement, the background potential gradient distribution is collected in an area assuming no concentrated leakage. After formal measurements, the abnormal potential gradient caused by the seepage channels was initially separated. ,in This represents the measured total potential gradient; The raw data collected by the three methods are all assigned a unified spatiotemporal reference coordinate. The observed values, the quality indicators generated based on the data quality judgment results, and the corresponding acquisition parameters are combined to package and generate a standard data package to form multi-source geophysical field data.

[0031] Specifically, the goal of the frequency-domain ground-to-air electromagnetic method is to acquire induced signals from the resistivity distribution of the subsurface medium. Its ground-based transmitting device is typically a long conductor or large loop transmitter that emits specifically coded electromagnetic wave signals into the subsurface. The transmitted waveform is composed of multiple sinusoidal components of specific frequencies, encoded using a fundamental frequency pseudo-random waveform. Multiple sets of fundamental frequency pseudo-random waveforms with different frequency combinations and power spectrum distributions are pre-stored, forming a coding group library. Before or during measurement, the environmental electromagnetic noise level of the current survey line is assessed in real time (primarily through short-term monitoring of background noise during periods of no transmission). Based on the assessment results, a waveform that best matches the current noise characteristics is adaptively selected from the coding group library for transmission. For example, in areas with strong industrial interference, a coding group with a dominant frequency that avoids interference bands and has more concentrated energy is selected; in quiet areas, a coding group with a wider bandwidth to acquire more information is selected. The electromagnetic receiving system, carried by an aircraft, flies along a pre-designed survey line network, continuously acquiring secondary electromagnetic field signals, i.e., electromagnetic response signals, induced by the subsurface medium. The aircraft's flight altitude is dynamically adjusted based on real-time terrain data (such as ground elevation obtained through lidar or real-time differential GPS), with the goal of maintaining a relatively constant altitude between the receiving system and the ground. For example, if the relative flight altitude is set to 50 meters, the aircraft's absolute altitude will rise when flying over hills to maintain a height of 50 meters above the ground, and will descend when flying over valleys.

[0032] The goal of nuclear magnetic resonance (NMR) groundwater detection is to directly detect the content and state of free groundwater. An excitation pulse at a Larmor frequency (determined by the local geomagnetic field strength) is applied to the ground via a ground-based transmitter loop (connected to a transmitting device). After the pulse is turned off, the same loop or another receiver loop acts as an antenna to receive the NMR signal generated by the relaxation of hydrogen protons in the groundwater. This signal decays over time, and its decay curve is recorded. First, the signal-to-noise ratio of the currently acquired signal is calculated in real time. The method involves selecting a time window in the early part of the attenuation curve (the portion where the signal is stronger) to calculate the average signal amplitude, and selecting a time window at the tail end of the curve (where the signal has attenuated to noise levels) to calculate the standard deviation of the noise. The ratio of the two is the average of the average signal amplitude and the average noise level. Then, an initial fit evaluation is performed, and a single exponential decay model is applied to the currently collected decay curve data points. The least squares fit, where It is the signal amplitude at time t. It is the initial amplitude. This is the longitudinal relaxation time; record the correlation coefficient of this fit. This is used to measure whether the decay curve conforms to an ideal single exponential decay pattern. Finally, two preset quality thresholds are used for threshold judgment; one is the signal-to-noise ratio threshold. The typical value range is between 3 and 10, for example, it is usually set to =5, to ensure the signal is clearly distinguishable; the other is the fitting correlation coefficient threshold. The typical value range is between 0.8 and 0.95, for example, setting... =0.9, to ensure reliable attenuation pattern.

[0033] Immediately after a single measurement sequence is completed, a judgment is made; if the condition is met... or If either of these two conditions is met, the currently acquired data quality is deemed substandard, automatically triggering operations to increase the number of superpositions and switch the pulse moment sequence. It automatically increases the average number of signal superpositions; for example, initially set to 4 superpositions, if this is insufficient, it increases to 8, 16, and so on, until the preset maximum superposition limit is reached (e.g., 64 superpositions). Superposition can effectively suppress random noise and improve... The pulse moment (the product of the excitation pulse duration and the current intensity) determines the detection depth and excitation volume; multiple pulse moment sequences are pre-stored (e.g., sequences varying from small to large for detecting shallow details, or uniformly distributed sequences for balanced detection); if the fit is poor... If the signal is too low, it indicates complex hydrogeological conditions within the current excitation volume, and the signal does not conform to single exponential decay. Try switching to another preset pulse moment sequence for remeasurement to observe whether a more regular decay curve can be obtained. If, after performing one or both of the above operations, the re-acquired data still cannot simultaneously meet the requirements... and If so, then the measuring point is marked as an interference point.

[0034] The flow field method (specifically, the natural electric field method or artificial source electric field method for seepage monitoring) aims to detect natural electric field or excitation potential anomalies generated by groundwater seepage. Before each pre-set profile measurement begins, the background potential gradient distribution is first measured; a section on the profile identified as having "no concentrated seepage zone" is selected, based on prior geological data or field investigation; within this section, the potential difference between two points is measured along the profile using non-polarized electrodes at a fixed point spacing (e.g., equal to the basic grid size), and divided by the point spacing to obtain the potential gradient. After the background measurement is completed, a formal potential gradient measurement is performed along the entire profile to obtain the total potential gradient distribution, including the background field and possible seepage anomaly fields. Preliminary identification of abnormal potential gradients that may be caused by seepage channels: The background scene here Considered as a stable electric field distribution under regional or normal conditions; calculated These are observations that are directly used in subsequent data inversion, directly reflecting information about local seepage activity.

[0035] All measurement data acquired by the three methods (FDEM flight trajectory points, NMR measurement point center coordinates, and flow field method electrode positions) were processed using Global Navigation Satellite System (GNSS) Real-Time Dynamic Differential (RTK) technology to obtain their three-dimensional coordinates, which were then unified to the coordinate system and elevation datum used by a pre-defined two-dimensional base grid. The time of all acquisition devices was synchronized using a high-precision time signal provided by a GNSS receiver, ensuring that data acquired by different methods had a unified time stamp. For the raw data acquired by each method at each measurement point or along each survey line, structured standard data packets were generated, forming multi-source geophysical field data. Each data packet contained at least observations, quality indicators, and acquisition parameters. Specifically, the physical quantity observations were processed after preliminary format conversion and unit unification, such as the complex electromagnetic field components of FDEM, the decay curve data sequence of NMR, and the potential gradient values ​​of the flow field method. Quality indicators were quality assessment information generated based on real-time or post-acquisition analysis during the acquisition process; for example, the flight attitude angle information of FDEM data and the signal-to-noise ratio of NMR data. Correlation coefficient This includes whether the data was marked as an "interference point," the electrode grounding resistance value for flow field method data, etc. The specific parameters collected include recording all relevant parameters for this observation; for example, the emission frequency, emission current, and flight altitude of the FDEM; the Larmor frequency, pulse moment sequence, and number of superpositions of the NMR; and the electrode arrangement and measurement time of the flow field method.

[0036] In one embodiment of the present invention, in step S2, a three-dimensional integrated initial geological model including resistivity, water content, and hydraulic conductivity parameters is constructed. Using the multi-source geophysical field data as joint constraints, the three-dimensional integrated initial geological model is iteratively optimized through an inversion framework that couples electromagnetic induction, nuclear magnetic resonance relaxation, and groundwater seepage physical mechanisms, including the following steps: Based on the multi-source geophysical field data and known hydrogeological information, a three-dimensional integrated initial geological model is constructed, and the target exploration area is discretized into... The model comprises three-dimensional mesh cells, each containing resistivity. Volumetric water content θ, longitudinal relaxation time The four physical properties are: hydroconductivity K; Within the inversion framework, a multiphysics collaborative inversion objective function is established, which is: Where m is the parameter vector of the three-dimensional integrated initial geological model. , , These are the data fitting differences for the frequency-domain ground-to-space electromagnetic method, the nuclear magnetic resonance groundwater detection method, and the flow field method, respectively. For model parameter smoothing constraints, For cross-gradient structure coupling constraint terms, and For regularization parameters; Among them, the cross-gradient structure coupling constraint term The calculation formula is: in, , , Let i and n represent the spatial gradient vectors of resistivity, water content, and hydraulic conductivity in the i-th 3D mesh cell, respectively, i=1,2,... , × represents the cross product of vectors; Using a three-dimensional integrated initial geological model as the starting point for iteration and multi-source geophysical field data as joint constraints, the objective function is minimized through optimization algorithms within the inversion framework. The model parameters are iteratively updated.

[0037] Specifically, based on the established two-dimensional basic grid, it is extended vertically and layered to form a three-dimensional spatial grid covering the entire volume of the open-pit coal mine exploration area. The target exploration area is discretized in three-dimensional space as follows: The three-dimensional grid cells can be regular or irregular, and their dimensions (such as side length) depend on the resolution of the target and the computing power. For the detection of water inflow channels in open-pit coal mines, the vertical grid layer thickness is typically between 2 and 10 meters, while the horizontal grid size is consistent with or a multiple of the aforementioned two-dimensional base grid (e.g., 10 meters × 10 meters × 5 meters). The total number of grid cells... The number could reach tens of thousands to hundreds of thousands, depending on the volume and resolution of the exploration area. After meshing, for each 3D mesh element i (where i = ... Initial physical property parameter values ​​are assigned to complete the construction of the three-dimensional integrated initial geological model, which simultaneously contains four physical property parameters to be inverted. Specifically, the parameters are resistivity... The unit is ohm-meter (Ω·m), which describes the electrical conductivity of a medium; volumetric water content. Dimensionless (or expressed as a percentage), describing the volumetric water content per unit volume of medium; longitudinal relaxation time. The unit is milliseconds (ms), which describes the rate of decay of nuclear magnetic resonance signals and is related to pore size and the state of water occurrence; hydroconductivity The units are meters per second (m / s) or meters per day (m / d), describing the ability of the medium to conduct groundwater. These parameters collectively constitute the components of the model parameter vector m in the three-dimensional grid cell. The assignment of parameters in the three-dimensional integrated initial geological model is based on prior information and multi-source geophysical field data, including regional hydrogeological survey reports, lithological columnar sections and well logging data from existing boreholes, and general rock physics empirical relationships; for example, the initial resistivity of a complete limestone area. Set at 500 Ω·m, initial water content Set to 0.05 (5%), initial relaxation time Set to 200ms, initial hydraulic conductivity K is set to m / s.

[0038] An inversion framework for collaborative inversion is established, the core of which is a multiphysics collaborative inversion objective function. The function is defined as follows: Where m represents all physical properties of all three-dimensional mesh elements ( , , The model parameter vector of (K). The data fitting difference in the frequency domain air-to-ground electromagnetic method measures the discrepancy between the electromagnetic response predicted by the current model parameter vector m (obtained through forward modeling of Maxwell's equations) and the actual observed electromagnetic data; it is typically expressed as a weighted least squares form, for example... ,in It is a weight matrix estimated based on data error. The data fit difference for nuclear magnetic resonance (NMR) is used to measure the water content in the current model parameter vector m. and relaxation time The difference between the predicted NMR decay curves (calculated using the NMR forward modeling formula) and the actual observed decay curve data. The data fit difference for the flow field method is measured by the conductivity in the current model parameter vector m. The distribution of groundwater conductivity K, and the predicted potential gradient and the actual observed anomalous potential gradient under the conditions of coupled groundwater seepage and current field. The differences between them. This is a model parameter smoothing constraint term, its purpose being to tend to obtain a model that varies gently in space, avoiding violent, non-physical oscillations, and thus stabilizing the inversion process; a common definition is the sum of squares of the differences in parameter values ​​between adjacent grid cells; for example, for resistivity... ,have ,in Represents the adjacent cells of 3D mesh cell i; overall smoothness constraint It is a weighted sum of the smoothing constraints of each of the four parameters. This is a cross-gradient structure coupling constraint term; it does not require different parameters (such as resistivity and water content) to be numerically correlated, but forces them to maintain a consistent spatial structure (i.e., variation trend); its specific calculation formula is as follows: ; where, for three-dimensional mesh element i, , , These are the spatial gradient vectors of resistivity, water content, and hydraulic conductivity at the center of the unit. These gradient vectors are approximated by the difference of model parameters in adjacent units, for example... ,in is the grid spacing; × represents the vector cross product. and The regularization parameter controls... and In the overall objective function The relative importance of them; and The value needs to be determined through trial and error or an adaptive method based on data error; a common strategy is to use the "L-curve" method, by plotting different ( , The relationship between the data fit difference and the magnitude of the model constraint terms under various values ​​is plotted, and the value at the inflection point is selected as the equilibrium solution; in typical cases, and The range of values ​​is within arrive Between; for example, the initial inversion setting =0.01, =0.05, and then fine-tuned according to the model update during the iteration process.

[0039] In one embodiment of the present invention, step S2, generating a three-dimensional multi-parameter fused attribute body, includes the following steps: During the iterative process of the inversion framework, the physical property relationship between resistivity, water content, and relaxation time is established through a rock physics bridging equation, which is as follows: ;in, The predicted conductivity value, Porosity estimated from water content θ and lithology For saturation, Let A be the electrical conductivity of pore water, and A be an empirical coefficient. and For cementation index and saturation index; The data fitting term of the flow field method Forward modeling is performed based on the coupled equations of the seepage field and the current field. The coupled equations are as follows: ;in, The distribution of groundwater head. For potential distribution, For the injected current intensity, Let r be the location of the current source point, and r be the spatial position vector. For the Dirac function, For gradient operators, and The spatial distributions of hydraulic conductivity and electrical conductivity, defined by the model parameter vector m, are respectively used. By solving the coupling equations, the predicted potential gradient considering the groundwater flow effect is obtained, and then the fitting difference with the measured data is calculated. The iterative process of its inversion framework employs a constrained optimization algorithm, updating the model parameter vector m in each iteration until the objective function is reached. The change in resistivity satisfies one of the following two conditions: less than a preset threshold or reaching the maximum number of iterations. The output is a three-dimensional multi-parameter fused attribute body containing the spatial distribution of resistivity, water content, relaxation time, and hydroconductivity.

[0040] Specifically, in the inversion iteration process, in order to utilize nuclear magnetic resonance data (reflecting water content θ and relaxation time) To constrain resistivity inversion and ensure physical consistency among different parameters, this method introduces a rock physics bridging equation connecting conductivity (the reciprocal of resistivity), water content, pore structure, and relaxation time within the inversion framework. This equation is expressed as follows: ,in, It is the conductivity value predicted by the model parameter vector m of the current iteration step, in Siemens per meter (S / m). The porosity of the underground medium, ranging from 0 to 1; in the inversion framework, Based on the volumetric water content in the current model When estimating lithological information using known or inferred lithological data, a common estimation relationship is... Alternatively, an empirical range can be assigned based on the lithology (e.g., sandstone porosity 0.15-0.35, limestone 0.01-0.20). Saturation refers to the volume percentage of water in the pore space, which is dimensionless (range 0-1). In most groundwater detection scenarios, the target is the saturation zone or capillary saturation zone, and it is usually assumed that... =1. Pore ​​water conductivity, measured in S / m, is controlled by groundwater salinity, temperature, and ionic composition. It can be obtained through laboratory testing of groundwater samples or by providing an empirical value based on regional hydrogeological data. Typical groundwater... The range is between 0.01 S / m (freshwater) and 1 S / m (saltwater). The cementation index is a dimensionless empirical constant that reflects the tortuosity of the pore structure. Its value is usually between 1.3 (unconsolidated sand) and 2.5 (highly cemented rock). is the saturation exponent, also a dimensionless empirical constant, reflecting the complexity of the current path in the unsaturated state, typically ranging from 1.5 to 2.5. A is an empirical coefficient with dimensions in S·s / m, used to express the relaxation time. The reciprocal of A is related to surface conductivity; the magnitude of A is related to the specific surface area of ​​the rock, clay mineral content, and cation exchange capacity, and is usually determined by fitting core experimental data or regional experience; the range of A is... The typical value is 0.5. This is the longitudinal relaxation time in the current model, measured in seconds (s). The predicted conductivity is calculated using the bridging equation. This predicted value can be compared with or used as a constraint by the conductivity model obtained by direct inversion and updating of frequency domain ground-to-air electromagnetic data, thereby introducing information on water content and pore structure into resistivity inversion and realizing the coupling of physical mechanisms.

[0041] Empirical parameters (cementation index) in rock physics bridging equations Saturation index (Empirical coefficient A) and pore water conductivity Calibration needs to be based on the specific geological and hydrological conditions of the exploration area. Calibration methods include: before inversion, conducting laboratory rock physics tests using core samples from existing boreholes within the survey area; and fitting data such as resistivity, porosity, saturation, and nuclear magnetic resonance relaxation time of the rock samples to determine... , The initial values ​​of A and B are given; or their ranges are given based on empirical relationships in regional rock physics. During the inversion process, these parameters can be used as known constants or as unknown variables in the inversion optimization. In addition, in step S4, measured hydrogeological data (such as lithology, groundwater salinity, in-situ hydraulic conductivity, etc.) obtained from verification boreholes can be used to correct these parameters through Bayesian updates, thereby improving their regional applicability.

[0042] Data fitting terms of flow field method The measurement compares the model-predicted seepage electric field response with the measured anomalous potential gradient. The differences; its forward modeling is based on a set of partial differential equations describing the coupling of the groundwater seepage field and the steady current field: The first equation is the steady-state groundwater seepage equation, where, denoted as the groundwater head distribution (unit: meters), which is the field variable to be determined; K(m) is the spatial distribution function of the hydraulic conductivity defined by the current model parameter vector m; The first equation is the gradient operator; this equation is solved under given boundary conditions (such as constant head boundary or impermeable boundary) to obtain the head distribution within the study area. The second equation is the steady-state current field equation under point current source conditions, where... The potential distribution (unit: volts) is another field variable to be determined; It is the spatial distribution function of conductivity defined by the current model parameter vector m; The intensity of the injected current (unit: amperes); is the spatial position vector of the current source point; r is the position vector in three-dimensional space; The Dirac function represents the current at the source point. Injection at the site. These two equations are derived through the medium properties K and It is related to m and solved spatially; its physical mechanism lies in the fact that the flow of groundwater in pores or fractures drags charged ions, generating a "filtering electric field" or "flow potential" related to the seepage velocity field; this additional electric field is superimposed on the artificial source The generated fields together constitute the observable total potential field; by solving these two equations simultaneously using numerical methods (such as the finite element method), the current hydraulic conductivity K(m) and electrical conductivity can be calculated. Under the model, the predicted value of the potential gradient caused by a specific seepage field distribution. Furthermore, Compared with the measured abnormal potential gradient Compare and calculate the data fit difference. .

[0043] The inversion process takes the constructed three-dimensional integrated initial geological model as the starting point for iteration. Using the multi-source geophysical field data as joint constraints, iterative optimization is performed within the inversion framework. In each iteration, based on the current model parameter vector... The theoretical predicted values ​​of the frequency domain ground-to-air electromagnetic response, nuclear magnetic resonance attenuation curve, and coupled seepage-electric field response were calculated, and the objective function was also calculated. The value of the objective function and its gradient are used. Optimization algorithms (such as the nonlinear conjugate gradient method, quasi-Newton method (L-BFGS), etc.) are employed to determine a direction for model updates based on the current objective function value and its gradient information. and step length This generates a new model parameter vector. The iteration continues until a preset termination condition is met. There are typically two termination conditions; iteration stops when either one is met: the objective function... The relative change is less than a preset threshold. ,For example This indicates that the optimization process has fully converged; the number of iterations has reached the preset maximum value. ,For example =100. The final model parameter vector obtained when the iteration terminates. Each three-dimensional mesh cell contains four defined physical property parameters: resistivity. (or conductivity) Volumetric water content θ, longitudinal relaxation time And the hydroconductivity K. This complete, internally tightly linked parametric three-dimensional data volume is the three-dimensional multi-parameter fused attribute volume.

[0044] In one embodiment of the present invention, step S3, based on hydrogeological principles, defines the multi-dimensional feature vector combination pattern corresponding to the water inflow channel in the three-dimensional multi-parameter fusion attribute volume, including the following steps: The basic physical property parameters of each three-dimensional mesh unit i are extracted from the three-dimensional multi-parameter fused attribute volume, and the derived features of each three-dimensional mesh unit i are calculated based on the basic physical property parameters to construct a seven-dimensional basic feature vector. Among them, resistivity Volumetric water content Longitudinal relaxation time Hydraulic conductivity , The coefficient of conductivity Spatial gradient magnitude, The water conductivity index, The risk index for water inrush. Among them, the water conductivity index The calculation formula is: ,in, The hydraulic gradient at point i in the three-dimensional mesh. Water inrush risk index The calculation formula is: ,in, It is a very small constant; Based on hydrogeological principles, characteristic fingerprint patterns of several typical geological targets are predefined, including fingerprints of active water-conducting channels, fingerprints of water-rich aquitards, fingerprints of weakly water-bearing fractures, and fingerprints of dry and intact rock strata.

[0045] Specifically, the four core inverted physical property parameters (resistivity) of the three-dimensional mesh element i are read from the three-dimensional multi-parameter fused property volume. Volumetric water content Longitudinal relaxation time Hydraulic conductivity Based on the fundamental physical property parameters, three derived features are further calculated: the spatial gradient modulus of the hydraulic conductivity, the hydraulic conductivity index, and the inrush risk index. Finally, for each three-dimensional mesh element i in this property volume, a seven-dimensional fundamental feature vector is constructed. .

[0046] Spatial gradient modulus of hydraulic conductivity This describes the drastic spatial variation of water conductivity. The calculation method involves first estimating the gradient components of the water conductivity K at the center of a three-dimensional mesh cell i in three spatial directions (x, y, z), for example, using the central difference method. ,in and It is the K value of adjacent grid cells along the positive and negative x-axis. The grid spacing is set to x-axis; the same applies to y and z-axis. Then, the magnitude (i.e., the size) of the gradient vector is calculated. High gradient modulus values ​​typically indicate abrupt boundaries in lithology or hydraulic properties, such as the edges of faults or fracture zones.

[0047] Hydraulic conductivity index This formula combines the water conductivity of the underground medium with the hydraulic gradient driving the water flow, and approximates the potential seepage intensity or flow capacity. The calculation formula is as follows: ;in, It is the magnitude of the hydraulic gradient of groundwater at point i in the three-dimensional grid. It is the spatial rate of change of the water head h, which can be obtained by solving the steady-state groundwater seepage equation established in step S2. Given boundary conditions, the head distribution h(r) across the entire model region is obtained, and then the gradient of each three-dimensional mesh element is calculated. (Hydroconductivity index) The higher the value, the more likely the location is to become an active groundwater runoff channel.

[0048] Water inrush risk index It is a comprehensive indicator used to preliminarily assess the potential risk of water inrush at a specific spatial location. Its calculation formula is as follows: Among them, the molecular part It combines water-bearing capacity (water content θ) and flow conductivity (hydraulic index). ); the denominator is resistivity Add a very small constant ;join in (For example, the value is) Ω·m) is to prevent when the resistivity A division-by-zero error occurs when the resistivity approaches zero. From a hydrogeophysical perspective, low resistivity is usually associated with high pore water content or high mineralization, which may imply stronger water abundance and a higher material basis for inrush. Therefore, the design logic of this index is that in areas rich in water and with high conductivity (large molecules) and low resistivity (potentially indicating good water-bearing channels), the inrush risk index... It will be even higher. This water inrush risk index calculation formula is designed for typical mining area hydrogeological conditions; for rare high-resistivity, pure freshwater aquifers, the water inrush risk is mainly reflected in the high water content. With high water conductivity index (i.e., the numerator of the formula), and can be derived through other characteristics (such as high hydraulic conductivity). It is then comprehensively represented in subsequent pattern recognition.

[0049] Ultimately, the features of each 3D mesh cell i are integrated into a seven-dimensional vector. .

[0050] Based on prior knowledge of hydrogeology and geophysics, several types of "characteristic fingerprint patterns" for typical underground geological-hydrological structures are predefined. These patterns are vectors or regions with specific numerical combinations and ranges in a seven-dimensional characteristic space, representing idealized geological conditions. The predefined typical patterns typically include fingerprints of active water-conducting channels, water-rich aquitards, weakly water-bearing fractures, and dry, intact rock strata. Active water-conducting channel fingerprints correspond to features such as water-bearing faults, dissolution conduits, and high-permeability fracture zones; their expected characteristics include low to medium resistivity (due to water abundance), high water content, medium to long relaxation time (indicating larger pores or fractures), high hydraulic conductivity, high hydraulic conductivity gradient (channel boundaries), high hydraulic conductivity index, and high risk of water inrush. The fingerprint of aquitards rich in water corresponds to strata such as mudstone and intact shale, which may have high water content but extremely poor permeability. Its expected characteristics include low resistivity (rich in water), high water content, short relaxation time (small pores), extremely low hydraulic conductivity, low hydraulic conductivity gradient (homogeneous), low hydraulic conductivity index (due to extremely low K value), and a medium-to-low risk index for water inrush (although rich in water, it cannot flow). The fingerprint of weakly water-bearing fractured strata corresponds to rock masses with well-developed microfractures, containing small amounts of water but with generally poor connectivity. Its expected characteristics include resistivity, water content, relaxation time, and hydraulic conductivity all within a moderate range; the hydraulic conductivity gradient may vary somewhat; and the hydraulic conductivity index and risk index for water inrush are moderate to low. The fingerprint of dry, intact rock strata corresponds to intact, dry bedrock. Its expected characteristics include high resistivity, extremely low water content, insignificant or extremely low relaxation time, extremely low hydraulic conductivity, low gradient, and extremely low hydraulic conductivity index and risk index for water inrush. The specific numerical range of these fingerprint patterns (i.e., the upper and lower bounds of each dimension in the seven-dimensional space) needs to be calibrated based on the specific rock physical properties, hydrogeological conditions, and historical data (if available) of the exploration area. For example, by using the logging and test data of a certain water-conducting fracture exposed by a known borehole, the statistical feature vector of the corresponding area in the three-dimensional multi-parameter fusion attribute volume can be inferred, which can be used as a reference benchmark for the fingerprint of active water-conducting channels.

[0051] In one embodiment of the present invention, in step S3, a machine learning algorithm is used to perform spatial pattern search and classification in a three-dimensional multi-parameter fusion attribute volume, automatically identify and delineate abnormal spatial regions that conform to the multi-dimensional feature vector combination pattern, mark them as suspected water inrush channels, and output a three-dimensional image of the identification results with channel spatial morphology, category, and confidence information, including the following steps: The fundamental feature vectors of all three-dimensional mesh elements The vectors in the seven-dimensional feature space are standardized and projected into a seven-dimensional feature space. A density clustering algorithm is used to perform cluster analysis on the vectors in the seven-dimensional feature space. Each three-dimensional grid cell is assigned a cluster label, and cells without cluster labels are marked as noise points. For each non-noise cluster Calculate its cluster average eigenvector. and calculate Based on the similarity distance between each feature fingerprint pattern, each cluster is initially classified into one of the following categories according to the nearest neighbor principle: suspected water inflow channel, water-rich layer, weakly water-bearing fracture, or dry rock layer. For clusters initially classified as suspected water inflow channels, after three-dimensional spatial morphological filtering, those with a volume greater than a preset threshold are selected based on three-dimensional neighborhood connectivity. Furthermore, the spatial distribution can form a connected domain with a continuous path from the potential replenishment boundary to the discharge area, and then the spatial attitude parameters of the connected domain can be extracted, including the principal direction, dip angle and thickness. For each finally determined suspected water inflow channel connectivity region Calculate the overall confidence score The calculation formula is as follows: ;in, For connected components Average eigenvector fingerprints with active water-conducting channels similarity, For connected components volume, The total volume of all suspected water inflow channels connected domains. For connectivity score, , , Preset weighting coefficients; The final result is a 3D image of the recognition outcome, which includes information on the spatial morphology, category, and confidence level of the channels.

[0052] Specifically, a seven-dimensional basic feature vector is constructed for each three-dimensional mesh cell i. Standardization is performed; due to the significant differences in the physical meaning and numerical range of each characteristic parameter (e.g., resistivity)... The values ​​range from tens to thousands of Ω·m, while the water content θ is a decimal between 0 and 1. Direct cluster analysis would lead to a large range of dominant feature distance calculations. Therefore, the Z-score normalization method is used to process each feature dimension separately. The calculation formula is as follows: ;in It is the original value of the i-th 3D mesh element in the j-th feature dimension. and These are the mean and standard deviation of all 3D mesh cells in this feature dimension, respectively. All standardized vectors constitute a seven-dimensional feature space point set.

[0053] Density-based spatial clustering (DBSCAN) is used to perform cluster analysis on a seven-dimensional feature space point set. DBSCAN does not require pre-specifying the number of clusters and can effectively identify noise points (isolated points), making it suitable for handling irregular, non-spherical cluster distributions in geological bodies. The algorithm requires two key parameters: neighborhood radius and minimum number of points. The neighborhood radius defines the range of neighborhoods within which a point is considered a core point, typically ranging from 0.1 to 3.0 (in the standardized feature space), and needs to be determined experimentally or empirically; for example, an initial value of 1.5 can be used. The minimum number of points defines the minimum number of points required in the neighborhood of a core point, typically ranging from 5 to 20, for example, set to 10. After the algorithm executes, each three-dimensional mesh cell is assigned a cluster label. (j=1,2,…, ),in It represents the number of clusters discovered; cells not assigned to any cluster are marked as noise points, which correspond to areas of poor data quality, attribute transitions, or extreme anomalies.

[0054] For each non-noise cluster (Include (3D mesh cells), calculate the average eigenvector within each cluster. : ;in This is the standardized feature vector. Then, the average vector is calculated. With each characteristic fingerprint pattern (active water channel fingerprint) Water-rich waterproof layer fingerprint Weakly water-bearing fracture fingerprints fingerprints of dry, intact rock strata The similarity distance between fingerprints is used here; Euclidean distance is used as the similarity metric, for example, the distance to the fingerprint of an active water channel is: The smaller the distance, the higher the similarity. Based on the nearest neighbor principle, clusters are... The preliminary classification is based on the geological category corresponding to the fingerprint pattern with the smallest distance, namely, one of the following: suspected water inflow channel, water-rich layer, weakly water-bearing fissure, or dry rock layer.

[0055] Clusters initially classified as potential water inflow channels may appear as multiple scattered clumps or complex shapes in three-dimensional space, requiring further screening to identify connected spatial volumes that conform to the geological significance of water inflow channels. First, three-dimensional morphological filtering is performed. All mesh elements initially identified as potential water inflow channels undergo a three-dimensional morphological opening operation. The opening operation first performs an erosion operation to remove small, isolated elements (potentially noise), then a dilation operation to restore the approximate shape of the remaining main body. This helps smooth boundaries and eliminate small spatial discrete points. The structural element is typically a 3×3×3 cube. Then, the morphologically filtered potential water inflow channel elements are labeled with three-dimensional connected components. A 26-neighborhood standard (i.e., adjacent elements in the top, bottom, left, right, front, back, and all diagonal directions of a three-dimensional mesh element are considered connected) is used to scan and mark all spatially connected three-dimensional mesh elements as the same connected component. (o=1,2,…, , (Total number of connected components). For each connected component... Calculate its volume (Equals the number of grid cells contained in the connected component multiplied by the volume of a single cell); Set a minimum volume threshold. ,For example =100 (Equivalent to a size of approximately 10 meters × 10 meters × 1 meter), used to filter out anomalies that are too small and have little geological significance. Simultaneously, it determines whether the spatial distribution of the connected domain can form a continuous path from the model region boundary (potential recharge area) to the target discharge area (such as a mine pit); this can be achieved by checking whether the connected domain intersects with both the preset recharge boundary grid and the discharge area boundary grid. Only when both conditions are met... Only connected regions with the potential to form continuous paths are retained. For each suspected inrush channel connected region ultimately retained... Using its spatial coordinate point cloud, its spatial attitude parameters are extracted through principal component analysis (PCA); the direction of the first principal component indicates the main direction (strike / dip) of the channel, and the corresponding eigenvalue ratio reflects its linear extension; based on the fitted spatial plane, its dip angle is calculated; the average span of the connected domain perpendicular to the main direction can be used as an estimate of its thickness.

[0056] For each finally determined suspected water inflow channel connectivity region Calculate a comprehensive confidence score This is used to quantify the reliability of the identification result, and its calculation formula is as follows: Among them, the similarity item It is a connected component Average eigenvector (Calculation method: cluster average vector) and active water channel fingerprint The similarity can be calculated using normalized inverse distance or cosine similarity, ensuring the value is between 0 and 1; a higher value indicates a greater resemblance to the typical channel pattern. Relative volume term. The volume of this connected region accounts for a portion of the total volume of all suspected water inflow channels. The proportion; a larger relative volume usually means the target is more salient. Connectivity score items This is a value between 0 and 1, used to assess the completeness of the connected domain in forming an effective seepage path. Its calculation can be based on factors such as the distance between the connected domain and the recharge and discharge boundaries, and the spatial continuity of its internal hydraulic conductivity. For example, if the connected domain connects both the recharge and discharge boundaries, it receives 1 point; if it only connects one side, it receives 0.5 points; otherwise, the score is lower. (Weighting coefficient) , , To adjust the relative importance of the three indicators, the following must be met: + + =1; a typical weight allocation could be... =0.5 (emphasizing the matching of physical properties) =0.2 (considering scale) =0.3 (emphasizing hydrogeological function). The final calculated value is... The value is between 0 and 1, with higher values ​​indicating a higher confidence level in identifying it as a genuine water inflow channel.

[0057] Finally, all processing results are integrated to generate a 3D map of the identification outcome. This map is a 3D visualization data volume, in which each 3D grid cell contains at least a category label and confidence information. The category label identifies the final geological category to which the 3D grid cell belongs (e.g., water inflow channel, aquifer, weakly aquifer fracture, dry rock layer, noise / unclassified). For 3D grid cells labeled as water inflow channels, the confidence information also stores the connected components to which they belong. Overall confidence score In addition, the connectivity of each water inflow channel is recorded as an additional layer or annotation. Spatial morphological parameters (principal orientation, dip angle, thickness) and their confidence scores.

[0058] In one embodiment of the present invention, step S4 includes the following steps: Based on the identified results, the connected domains corresponding to each suspected water inflow channel in the 3D image are... Confidence score And the posterior uncertainty estimation of each parameter in the three-dimensional multi-parameter fused attribute volume M; the optimal location of the verification borehole is determined through a multi-objective optimization function. : Where J(·) is the indicator function, when the optimal borehole position Located in the connected domain The value is 1 if the time condition is met, and 0 otherwise. , , For preset weighting coefficients, For parameter p at position Posterior uncertainty estimation at the location, The three-dimensional spatial boundary of the target exploration area. To characterize from the boundary The inferred seepage path to the mine pit was determined by the borehole location. Path connectivity function for verification degree; In the optimal position Drilling was carried out at the site to obtain measured hydrogeological data. This includes lithological columnar sections, stratified water inflow Q, hydrostatic pressure P, and hydraulic conductivity obtained through in-situ hydrogeological tests. Using measured hydrogeological data The three-dimensional multi-parameter fused attribute volume M is corrected using a Bayesian inversion framework to generate an updated three-dimensional multi-parameter fused attribute volume. Its Bayesian inversion framework uses the forward mapping function. Building borehole data The likelihood function of the 3D multi-parameter fused attribute volume M Using the three-dimensional multi-parameter fused attribute volume M as the prior model, the posterior model distribution is obtained according to Bayes' theorem. , For all measured hydrogeological data The dataset; by sampling the posterior model distribution, an updated model set is obtained, and its mean is calculated as... ; where, the forward mapping function The three-dimensional multi-parameter fused attribute volume M is located at... The parameters at a given location are mapped to predicted values ​​of the z-th type of observation data. For stratified inflow rate Q and hydrostatic pressure P, and The results were obtained through local groundwater numerical simulation and solving the seepage field equations, respectively. Based on the borehole verification results, the multi-dimensional feature vector combination pattern is adaptively optimized, and the corresponding areas are selected according to the actual geological conditions revealed by the borehole. The corresponding grid feature vector The samples are classified and categorized into the confirmed channel sample library and other corresponding predefined geological category sample libraries. Using the updated sample library, the feature fingerprint pattern definition in the multidimensional feature vector combination pattern is updated by updating the corresponding model parameters.

[0059] Specifically, based on the 3D image of the identification results and the 3D multi-parameter fusion attribute volume, the optimal placement of one or more verification boreholes is determined through multi-objective optimization. The optimization goal is to enable boreholes to verify suspected water inflow channels to the greatest extent possible, reduce the uncertainty of key model parameters, and reveal potential seepage paths from the area boundary to the mine pit as much as possible. Determining the optimal borehole location is crucial. The multi-objective optimization function is defined as: Among them, suspected channel verification items This encourages drilling to be placed within identified suspected water inflow channels, with priority given to channels with high confidence levels; This represents the connectivity region of the o-th suspected water inrush channel. It is its overall confidence score; It is an indicator function, when the candidate borehole location Located in the connected domain When it is within the three-dimensional space range, its value is 1; otherwise it is 0. It corresponds to the connected component. The preset weighting coefficients can be set to the same as... Proportional or set separately according to channel size, for example Model uncertainty reduction term This encourages borehole placement in regions of high posterior uncertainty in key physical property parameters within the 3D multi-parameter fused attribute volume M, aiming to minimize the uncertainty of the overall model through borehole measurement data. It is the parameter p (resistivity) One of the following three factors (water content θ, hydraulic conductivity K) is located at... The posterior uncertainty estimate can be used to estimate the variance of parameters obtained through covariance matrix analysis or Monte Carlo sampling. The larger the value, the higher the uncertainty of the parameter at that position; This is used to balance the importance of validating known objectives with exploring unknown uncertainties. Seepage path validation item. This encourages borehole placement at key locations where the inferred seepage path from the area boundary to the mine discharge zone can be verified; It represents the three-dimensional spatial boundary of the target exploration area, specifically the boundary portion related to groundwater recharge; It is a path connectivity function used to characterize the path from the boundary. The inferred seepage path to the mine pit was determined by the borehole location. The degree of verification; one specific calculation method is based on the hydraulic conductivity K and the connected domain. Several high-probability seepage paths from the recharge boundary to the discharge zone are generated through groundwater numerical simulation or graph theory methods; then, candidate borehole locations are calculated. The function value is the reciprocal of the shortest distance to these simulated paths, or the membership degree; the closer the distance, the higher the function value. If the function falls directly on a simulated path, it can take the maximum value (e.g., 1). This is the preset weighting coefficient for this item. Weighting coefficient , , The typical value range is between 0 and 1, and it usually needs to be normalized or its relative proportion set based on experience; for example, if high-confidence channels are to be verified first, it can be set to... As the primary weight (e.g., the total accounts for 0.6). and Each accounts for 0.2. The optimal verification borehole location can be obtained by solving a multi-objective optimization function (using grid search or optimization algorithms). .

[0060] At the determined optimal drilling location Drilling was carried out at the site to obtain a complete set of measured hydrogeological data. These data typically include lithological columnar sections (describing the lithology of strata at different depths), stratified water yield Q (obtained through segmented pumping or pressure tests, in cubic meters per day (m³ / d)), hydrostatic pressure P (obtained through borehole water level observations, in meters of water column (m)), and hydraulic conductivity obtained from in-situ hydrogeological tests. (For example, obtained through analysis of single-hole pumping tests, in meters per day (m / d)).

[0061] Using the hydrogeological data measured in these boreholes The original 3D multi-parameter fused attribute volume M is corrected using a Bayesian inversion framework to generate an updated attribute volume. The Bayesian framework first establishes the forward mapping function. The three-dimensional multi-parameter fused attribute volume M is located at the optimal drilling position. The parameter at point j is mapped to the predicted value of the j-th type of observation data. For the hydraulic conductivity... Forward mapping is the most direct method, which involves extracting M at the optimal borehole location. At the corresponding depth segment value as predicted value .

[0062] For stratified inflow rate Q and hydrostatic pressure P, the forward mapping and This needs to be achieved by constructing and running a separate three-dimensional groundwater numerical model (called a local model) specifically designed to simulate the borehole hydrological test. Specifically, this is done using the optimal borehole location... Centered on a local region (typically with a horizontal range) within the 3D multi-parameter fused attribute volume M, a local region is extracted. A local model is constructed within a 50-100 meter radius around the borehole (vertically consistent with M). This local model features a finer mesh near the borehole (e.g., horizontal spacing 1-2 meters, vertical spacing 0.5-1 meter). Model parameters (hydraulic conductivity, porosity, etc.) are obtained from the parameter values ​​of M in this region through three-dimensional linear interpolation. Boundary conditions are determined based on the regional hydrogeological conditions: if the local model boundary is far from the borehole, the head distribution obtained from the inversion of M is used as a constant head boundary; otherwise, a zero-flow boundary or a given head gradient boundary is used. In the local model, source and sink terms are set according to the actual conditions of the field experiment. For the stratified inflow Q, pumping wells are set at the corresponding depth, using the actual pumping flow rate as a constant flow condition. The drawdown in the well is simulated by solving the steady-state (or unsteady-state) seepage equation, thereby calculating the predicted inflow. For the hydrostatic pressure P, the steady-state seepage equation under the natural flow field is directly solved to obtain the head value at the borehole location and convert it into hydrostatic pressure. Based on forward mapping, a data... The likelihood function of the 3D multi-parameter fused attribute volume M It is usually assumed that the observation error follows a Gaussian distribution, then the likelihood function is expressed as the product of the prediction errors of each data point.

[0063] Using the original three-dimensional multi-parameter fused attribute volume M as the prior model, the prior distribution is obtained. According to Bayes' theorem, the posterior distribution of the three-dimensional multi-parameter fused attribute volume M is: ;in This represents the dataset of all verified boreholes. Since the 3D multi-parameter fusion attribute volume M is high-dimensional and nonlinear, directly solving the posterior distribution analytically is extremely difficult. Therefore, Markov Chain Monte Carlo (MCMC) methods or ensemble smoothers (such as Ensemble Smoother) are used to sample the posterior distribution. These algorithms generate a series of model samples that conform to the posterior distribution. The mean (or median) of these posterior model samples is calculated as the updated 3D multi-parameter fused attribute body that incorporates borehole data constraints. Compared with the prior model M, The model better matches the measured data at the borehole location and within its influence range, and the uncertainty of the model is reduced.

[0064] Based on the actual geological conditions revealed by borehole verification, the predefined multi-dimensional feature vector combination pattern (i.e., feature fingerprint pattern) is adaptively optimized to make the pattern definition closer to the actual geological conditions. For each stratigraphic or section with clear hydrogeological significance revealed by the verification borehole (such as fracture zones confirmed as water-conducting channels, water-rich but impermeable mudstone layers, dry and intact rock layers, etc.), its updated three-dimensional multi-parameter fusion attribute volume is extracted. The three-dimensional grid cells corresponding to the spatial locations are identified. The seven-dimensional basic feature vectors of these three-dimensional grid cells are calculated. These feature vectors are then added to the corresponding confirmed geological category sample library according to their actual geological category. For example, samples confirmed as water-conducting channels are added to the "Confirmed Channel Sample Library," and samples confirmed as aquatic-rich impermeable layers are added to the "Aquatic-Rich Impermeable Layer Sample Library."

[0065] Using the updated and expanded sample databases for each geological category, the statistical characteristics of each category are recalculated. Typically, the mean vector of all feature vectors in the sample database for each category is calculated and used as the updated feature fingerprint pattern for that category. For example, the updated "active water-conducting channel fingerprint". This is the average of all feature vectors in the "confirmed channel sample library". Furthermore, the covariance matrix can be calculated to describe the dispersion of the fingerprint pattern. In this way, the definition of the feature fingerprint pattern evolves from a "predefined" one based on prior knowledge to a "data-driven definition" based on actual verification data, thus more accurately reflecting the specific geological and geophysical relationships of the exploration area.

[0066] The above embodiments are only used to illustrate the technical methods of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical methods of the present invention without departing from the spirit and scope of the technical methods of the present invention.

Claims

1. A method for recognizing water inrush channel fusion based on nuclear magnetic resonance and ground-air electromagnetic cooperation, characterized in that, Includes the following steps: S1: Within the target exploration area, simultaneously implement frequency domain ground-to-space electromagnetic method planar scanning, nuclear magnetic resonance groundwater detection method point measurement, and flow field method profile measurement to obtain multi-source geophysical field data; S2: Construct a three-dimensional integrated initial geological model that includes resistivity, water content and hydraulic conductivity parameters; using the multi-source geophysical field data as joint constraints, iteratively optimize the three-dimensional integrated initial geological model through an inversion framework that couples electromagnetic induction, nuclear magnetic resonance relaxation and groundwater seepage physical mechanisms to generate a three-dimensional multi-parameter fused attribute body. S3: Based on hydrogeological principles, define the multi-dimensional feature vector combination pattern corresponding to the water inflow channel in the three-dimensional multi-parameter fusion attribute body; use machine learning algorithms to search and classify spatial patterns in the three-dimensional multi-parameter fusion attribute body, automatically identify and delineate abnormal spatial areas that conform to the multi-dimensional feature vector combination pattern, mark them as suspected water inflow channels, and output a three-dimensional map of the identification results with channel spatial morphology, category and confidence information. S4: Deploy verification boreholes at key locations of suspected water inflow channels indicated by the three-dimensional map of the identification results; use the measured hydrogeological data revealed by the boreholes to correct and optimize the three-dimensional multi-parameter fusion attribute volume and the feature vector combination mode.

2. The water inrush channel fusion identification method based on nuclear magnetic resonance and ground-to-air electromagnetic synergy as described in claim 1, characterized in that, In step S1, within the target exploration area, frequency-domain ground-to-air electromagnetic surface scanning, nuclear magnetic resonance groundwater detection point measurement, and flow field profiling are simultaneously implemented, including the following steps: The target exploration area of ​​open-pit coal mines is divided based on a unified two-dimensional basic grid. Based on the preset geological structure information and hydrological model, a survey line network for frequency domain ground-to-air electromagnetic method surface scanning is automatically generated on the two-dimensional basic grid. In particular, in known anomalous areas and anomalous areas inferred from prior geological information, encrypted verification survey lines are generated. The location of groundwater detection points is dynamically determined by optimizing a function, which is: ,in, Optimized set of measurement point locations, It is a set of known geological structural location information. This is a set of potential anomaly regions predicted by frequency-domain ground-to-space electromagnetic methods based on geological structural information and hydrological models. This is the set of all possible measurement point locations within the target exploration area. Let be the spatial density function. This is a spatial overlap function. For spatial coverage function, , , These are the weighting coefficients; based on the optimization function. From the optimization results, the measuring points with the highest overlap with geological structures and predicted anomalies were selected as key verification measuring points. Lay out the survey lines for flow field method profile measurement, ensuring that the survey lines pass through key verification points and the areas where the denser verification survey lines are located.

3. The water inrush channel fusion identification method based on nuclear magnetic resonance and ground-to-air electromagnetic synergy as described in claim 2, characterized in that, In step S1, acquiring multi-source geophysical field data includes the following steps: For the frequency domain ground-to-air electromagnetic method, a ground-based transmitting device transmits specifically coded electromagnetic waves into the ground; an aircraft carrying an electromagnetic receiving system flies along the designed survey line to collect the electromagnetic response signals of the underground medium; the transmitted electromagnetic waves adaptively select the fundamental frequency pseudo-random waveform coding group according to the environmental noise level; and the flight altitude of the aircraft is dynamically adjusted according to real-time terrain data. For nuclear magnetic resonance (NMR) groundwater detection, a ground-based transmitter and electromagnetic receiver system are used to acquire NMR signals and the resulting attenuation curves. The signal-to-noise ratio of the NMR signals is calculated in real time during the acquisition process. The decay curve was initially fitted using a single exponential method to obtain the correlation coefficient. Set the signal-to-noise ratio threshold. Correlation coefficient threshold ,achieve and If either of the two conditions is met, the operation of increasing the number of superpositions and switching the pulse moment sequence will be automatically triggered; if the threshold condition is still not met after the operation is triggered, the measurement point will be marked as an interference point. For the flow field method, before each profile measurement, the background potential gradient distribution is collected in an area assuming no concentrated leakage. After formal measurements, the abnormal potential gradient caused by the seepage channels was initially separated. ,in This represents the measured total potential gradient; The raw data collected by the three methods are all assigned a unified spatiotemporal reference coordinate. The observed values, the quality indicators generated based on the data quality judgment results, and the corresponding acquisition parameters are combined to package and generate a standard data package to form multi-source geophysical field data.

4. The water inrush channel fusion identification method based on nuclear magnetic resonance and ground-to-air electromagnetic coordination according to claim 1, characterized in that, In step S2, a three-dimensional integrated initial geological model is constructed, including resistivity, water content, and hydraulic conductivity parameters. Using the multi-source geophysical field data as joint constraints, the three-dimensional integrated initial geological model is iteratively optimized through an inversion framework that couples electromagnetic induction, nuclear magnetic resonance relaxation, and groundwater seepage physical mechanisms. This optimization includes the following steps: Based on the multi-source geophysical field data and known hydrogeological information, a three-dimensional integrated initial geological model is constructed, and the target exploration area is discretized into... The model comprises three-dimensional mesh cells, each containing resistivity. Volumetric water content θ, longitudinal relaxation time The four physical properties are: hydroconductivity K; Within the inversion framework, a multiphysics collaborative inversion objective function is established, which is: Where m is the parameter vector of the three-dimensional integrated initial geological model. , , These are the data fitting differences for the frequency-domain ground-to-space electromagnetic method, the nuclear magnetic resonance groundwater detection method, and the flow field method, respectively. For model parameter smoothing constraints, For cross-gradient structure coupling constraint terms, and For regularization parameters; Among them, the cross-gradient structure coupling constraint term The calculation formula is: in, , , Let i and n represent the spatial gradient vectors of resistivity, water content, and hydraulic conductivity in the i-th 3D mesh cell, respectively, i=1,2,... , × represents the cross product of vectors; Using a three-dimensional integrated initial geological model as the starting point for iteration and multi-source geophysical field data as joint constraints, the objective function is minimized through optimization algorithms within the inversion framework. The model parameters are iteratively updated.

5. The water inrush channel fusion identification method based on nuclear magnetic resonance and ground-to-air electromagnetic coordination according to claim 4, characterized in that, In step S2, generating a three-dimensional multi-parameter fused attribute volume includes the following steps: During the iterative process of the inversion framework, the physical property relationship between resistivity, water content, and relaxation time is established through a rock physics bridging equation, which is as follows: ;in, The predicted conductivity value, Porosity estimated from water content θ and lithology For saturation, Let A be the electrical conductivity of pore water, and A be an empirical coefficient. and The cementation index and saturation index are used for bonding and saturation. The data fitting term of the flow field method Forward modeling is performed based on the coupled equations of the seepage field and the current field. The coupled equations are as follows: ;in, The distribution of groundwater head. For potential distribution, For the injected current intensity, Let r be the location of the current source point, and r be the spatial position vector. For the Dirac function, For gradient operators, and The spatial distributions of hydraulic conductivity and electrical conductivity, defined by the model parameter vector m, are respectively used. By solving the coupling equations, the predicted potential gradient considering the groundwater flow effect is obtained, and then the fitting difference with the measured data is calculated. The iterative process of its inversion framework employs a constrained optimization algorithm, updating the model parameter vector m in each iteration until the objective function is reached. The change in resistivity satisfies one of the following two conditions: less than a preset threshold or reaching the maximum number of iterations. The output is a three-dimensional multi-parameter fused attribute body containing the spatial distribution of resistivity, water content, relaxation time, and hydroconductivity.

6. The water inrush channel fusion identification method based on nuclear magnetic resonance and ground-to-air electromagnetic coordination according to claim 1, characterized in that, In step S3, based on hydrogeological principles, the multi-dimensional feature vector combination pattern corresponding to the water inflow channel in the three-dimensional multi-parameter fused attribute volume is defined, including the following steps: The basic physical property parameters of each three-dimensional mesh unit i are extracted from the three-dimensional multi-parameter fused attribute volume, and the derived features of each three-dimensional mesh unit i are calculated based on the basic physical property parameters to construct a seven-dimensional basic feature vector. Among them, resistivity Volumetric water content Longitudinal relaxation time Hydraulic conductivity , The coefficient of conductivity Spatial gradient magnitude, The water conductivity index, The risk index for water inrush; Among them, the water conductivity index The calculation formula is: ,in, The hydraulic gradient at point i in the three-dimensional mesh. Water inrush risk index The calculation formula is: ,in, It is a very small constant; Based on hydrogeological principles, characteristic fingerprint patterns of several typical geological targets are predefined, including fingerprints of active water-conducting channels, fingerprints of water-rich aquitards, fingerprints of weakly water-bearing fractures, and fingerprints of dry and intact rock strata.

7. The water inrush channel fusion identification method based on nuclear magnetic resonance and ground-to-air electromagnetic coordination according to claim 6, characterized in that, In step S3, machine learning algorithms are used to search and classify spatial patterns in a three-dimensional multi-parameter fusion attribute volume, automatically identifying and delineating abnormal spatial regions that conform to multi-dimensional feature vector combination patterns, marking them as suspected water inrush channels, and outputting a three-dimensional image of the identification results with channel spatial morphology, category, and confidence information, including the following steps: The fundamental feature vectors of all three-dimensional mesh elements The vectors in the seven-dimensional feature space are standardized and projected into a seven-dimensional feature space. A density clustering algorithm is used to perform cluster analysis on the vectors in the seven-dimensional feature space. Each three-dimensional grid cell is assigned a cluster label, and cells without cluster labels are marked as noise points. For each non-noise cluster Calculate its cluster average eigenvector. and calculate Based on the similarity distance between each feature fingerprint pattern, each cluster is initially classified into one of the following categories according to the nearest neighbor principle: suspected water inflow channel, water-rich layer, weakly water-bearing fracture, or dry rock layer. For clusters initially classified as suspected water inflow channels, after three-dimensional spatial morphological filtering, those with a volume greater than a preset threshold are selected based on three-dimensional neighborhood connectivity. Furthermore, the spatial distribution can form a connected domain with a continuous path from the potential replenishment boundary to the discharge area, and then the spatial attitude parameters of the connected domain can be extracted, including the principal direction, dip angle and thickness. For each finally determined suspected water inflow channel connectivity region Calculate the overall confidence score The calculation formula is as follows: ;in, For connected components Average eigenvector fingerprints with active water-conducting channels similarity, For connected components volume, The total volume of all suspected water inflow channels connected to each other. For connectivity score, , , Preset weighting coefficients; The final result is a 3D image of the recognition outcome, which includes information on the spatial morphology, category, and confidence level of the channels.

8. The method for fusion identification of water inrush channels based on nuclear magnetic resonance and ground-to-air electromagnetic coordination according to claim 1, characterized in that, S4 includes the following steps: Based on the identified results, the connected regions corresponding to each suspected water inrush channel in the 3D image are... Confidence score And the posterior uncertainty estimation of each parameter in the three-dimensional multi-parameter fused attribute volume M; the optimal location of the verification borehole is determined through a multi-objective optimization function. : Where J(·) is the indicator function, when the optimal borehole position Located in the connected domain The value is 1 if the time condition is met, and 0 otherwise. , , For preset weighting coefficients, For parameter p at position Posterior uncertainty estimation at the location, The three-dimensional spatial boundary of the target exploration area. To characterize from the boundary The inferred seepage path to the mine pit was determined by the borehole location. Path connectivity function for verification degree; In the optimal position Drilling was carried out at the site to obtain measured hydrogeological data. This includes lithological columnar sections, stratified water inflow Q, hydrostatic pressure P, and hydraulic conductivity obtained through in-situ hydrogeological tests. Using measured hydrogeological data The three-dimensional multi-parameter fused attribute volume M is corrected using a Bayesian inversion framework to generate an updated three-dimensional multi-parameter fused attribute volume. Its Bayesian inversion framework uses the forward mapping function. Building borehole data The likelihood function of the 3D multi-parameter fused attribute volume M Using the three-dimensional multi-parameter fused attribute volume M as the prior model, the posterior model distribution is obtained according to Bayes' theorem. , For all measured hydrogeological data The dataset; by sampling the posterior model distribution, an updated model set is obtained, and its mean is calculated as... Among them, the forward mapping function maps the three-dimensional multi-parameter fused attribute volume M to the position... The parameters at a given location are mapped to predicted values ​​of the z-th type of observation data. For stratified inflow rate Q and hydrostatic pressure P, and The results were obtained through local groundwater numerical simulation and solving the seepage field equations, respectively. Based on the borehole verification results, the multi-dimensional feature vector combination pattern is adaptively optimized, and the corresponding area is determined according to the actual geological conditions revealed by the borehole. The corresponding grid feature vector The samples are classified and categorized into the confirmed channel sample library and other corresponding predefined geological category sample libraries. Using the updated sample library, the feature fingerprint pattern definition in the multidimensional feature vector combination pattern is updated by updating the corresponding model parameters.

Images