Well-ground cross-medium electromagnetic detection and surrounding rock parameter inversion method and system

Abstract

The application discloses a well-ground cross-medium electromagnetic detection and surrounding rock parameter inversion method and system, and relates to the field of advanced detection of underground engineering.According to the electrical property of the target stratum and the transmission condition, a grounding electrode group or a non-grounding loop is arranged in the well, and a composite excitation current containing a frequency domain stable section and a time domain off section is output; a surface array obtains a frequency domain electric field response and a time domain magnetic field decay response under a unified time base.According to the roadway geometry, the metal structure, the terrain elevation and the parameters of the transmitting and receiving device, an engineering space coupling constraint model is established, the observation data is preprocessed and near-source interference is decoupled, and three-dimensional resistivity and polarization rate distribution are inversed jointly.According to the polarization decay half-time of the low-resistance abnormal area, a target geological body and its boundary are identified, and drilling advancing speed, rotation torque and tunneling direction adjustment parameters are generated.The application can reduce the interference of the underground engineering structure on the well-ground electromagnetic data, and improve the reliability of the surrounding rock physical property identification and tunneling decision.

Description

Technical Field

[0001] This invention relates to the field of geophysical exploration and construction control technology for underground engineering, specifically to a method and system for cross-medium electromagnetic exploration and surrounding rock parameter inversion in wells and on the ground. Background Technology

[0002] During the excavation or mining of deep mine roadways, tunnels, and other underground engineering projects, it is necessary to promptly grasp the lithological changes, water-bearing fractured zones, abnormal conductors, and ore body boundaries in front of and around the working face. Current construction methods typically combine preliminary geological data, core drilling, and equipment load feedback to determine the drilling route and process parameters. However, as the operating depth increases, discrepancies may arise between the initial data and the actual exposed strata. Relying solely on mechanical feedback such as propulsion resistance and torque usually only allows for a response after the drill bit approaches or enters an abnormal area.

[0003] Well-to-ground electromagnetic detection can shorten the propagation distance between the excitation source and deep targets by utilizing underground excitation sources. However, engineering structures such as underground tunnel cavities, tracks, beams, anchor bolts, and pipelines can alter the electromagnetic field distribution. Surface elevation variations, industrial power frequency noise, differences in transmit and receive geometry, and differences in the dimensions and sensitivity of frequency domain data versus time domain data can also affect inversion stability. Without a unified time base, engineering spatial constraints, and a joint inversion mechanism adapted to signal quality, it is easy to misinterpret near-source metal responses or topographic effects as geological anomalies, making it difficult to directly generate surrounding rock properties and boundary information usable by tunneling equipment. Summary of the Invention

[0004] This invention aims to provide a method and system for cross-medium electromagnetic detection and surrounding rock parameter inversion in wells and on the ground. By using underground source composite excitation, synchronous acquisition of surface array, engineering spatial constraint decoupling, time-frequency joint inversion, and geological body discrimination, the influence of tunnel engineering structure and survey area environment on the observation results is reduced, and the inverted surrounding rock properties and boundaries are converted into drilling operation parameters.

[0005] In a first aspect, the present invention provides a method for well-to-surface cross-medium electromagnetic detection and surrounding rock parameter inversion, comprising the following steps:

[0006] Underground source deployment and composite excitation: Based on the location of the underground roadway in the target mining area, the electrical parameters of the target stratum, and the transmission power supply parameters, determine the deployment method of the grounding electrode group or the ungrounded return line, and output a composite excitation current that includes a frequency domain stable segment and a time domain turn-off segment;

[0007] Surface array time-frequency acquisition: A unified time reference is established between the underground source and the surface array. Frequency domain electric field amplitude and phase data corresponding to the frequency domain stable segment are acquired, as well as time domain magnetic field data after the time domain cutoff segment ends. Weighted observation data are formed based on the signal-to-noise ratio of each frequency point and each time window.

[0008] Engineering constraint decoupling and joint inversion: Obtain the three-dimensional geometry of the tunnel, the spatial location of the metal structure, the surface elevation, and the geometric parameters of the transceiver equipment. Correct the observation data for gross errors, baselines, trends, topography, and equipment coefficients. Based on the above spatial parameters, establish an engineering spatial coupling constraint model. Under the constraints of the model, decouple the frequency domain electric field data and the time domain magnetic field data from near-source interference and perform joint inversion to obtain the three-dimensional resistivity distribution and the three-dimensional polarizability distribution.

[0009] Geological body identification and operation parameter generation: Low-resistivity anomaly areas are identified by three-dimensional resistivity distribution, the corresponding polarization attenuation sequence is fitted and the polarization attenuation half-life is determined, the target geological body attributes are determined based on the half-life, the three-dimensional boundary is extracted, and adjustment parameters for drilling speed, rotation torque and tunneling direction are generated.

[0010] In one implementation, the allowable equivalent total grounding resistance is determined based on the transmission supply voltage and the target transmission current, and the number of parallel connections is determined based on the measured independent grounding resistance of each grounding electrode; the grounding electrode group is arranged along the roadway axis and maintains a preset distance from the metal structure, and an ungrounded return line is used when the grounding conditions do not meet the requirements.

[0011] In one implementation, the composite excitation current adopts a variable duty cycle sequence, with the frequency domain stable segment being a square wave flat-top segment and the time domain turn-off segment being a bipolar trapezoidal wave falling edge segment; the turn-off time is determined by combining the target detection depth and the apparent resistivity of the measurement area, and the turn-off time is controlled within 20 microseconds.

[0012] In one implementation, the surface array is gridded according to the effective range determined by the surface projection position of the underground roadway and the background noise at the site. Each acquisition station is equipped with mutually orthogonal horizontal electric field sensors and three-dimensional magnetic field sensors, and the transmission and reception time are aligned by satellite timing signals or high-stability crystal oscillators.

[0013] In one implementation, the preprocessing includes abnormal waveform removal, superposition of repeated observations, baseline zeroing, polynomial trend subtraction, power frequency and harmonic suppression, terrain elevation normalization, sensor attitude correction, and correction of transceiver geometric factors and device coefficients.

[0014] In one implementation, the engineering spatial coupling constraint model includes grid cells corresponding to tunnel cavities, surrounding rock interfaces, metal tracks, canopy beams, and support components, and sets initial physical properties and boundary conditions based on the electrical parameters of various media.

[0015] In one implementation, the joint inversion uses an objective functional consisting of a data fitting term and a model smoothing constraint term to weight the observed data based on the data covariance matrix determined by the signal-to-noise ratio, and updates the regularization parameter during the iteration process.

[0016] In one implementation, the polarization decay sequence of the low resistivity anomaly region is fitted using a single exponential relaxation model, the polarization decay half-life is calculated based on the fitted time constant, and then compared with the half-life intervals of different geological bodies calibrated by the samples.

[0017] In one implementation, an ellipsoidal envelope is fitted to the target geological body to be avoided to form a safe envelope surface. Candidate rays are generated from the current drill bit position, and the tunneling direction is determined based on the safe avoidance margin and the deflection angle relative to the original tunneling axis.

[0018] In one implementation, signal-to-noise ratio (SNR) assessment is performed for each frequency point and each time window; data below a hard threshold are discarded, and data in the transitional SNR range are weighted using a continuous weighting function to reduce the impact of low-quality data on the inversion convergence direction.

[0019] Secondly, the present invention provides a well-to-surface cross-medium electromagnetic detection and surrounding rock parameter inversion system, comprising:

[0020] The underground composite excitation subsystem is used to drive the output of composite excitation current from the grounded electrode group or the ungrounded return line and to provide a time reference.

[0021] The surface array acquisition subsystem is used to acquire frequency domain electric field data and time domain magnetic field data, and to perform signal quality assessment on the observed data;

[0022] The Engineering Constraint Inversion and Operation Control Center is used to establish engineering spatial coupling constraint models, complete data preprocessing, near-source interference decoupling, resistivity and polarizability joint inversion, geological body attribute identification and three-dimensional boundary extraction, and output operation parameter adjustment instructions to downhole tunneling equipment.

[0023] Compared with the prior art, the present invention has at least the following beneficial effects:

[0024] This invention places underground composite excitation and surface array time-frequency acquisition under a unified time reference, enabling the frequency domain electric field response and time domain magnetic field attenuation response to enter the same data processing chain. By introducing an engineering spatial coupling constraint model formed by tunnel cavities, surrounding rock interfaces, metal structures, terrain elevation, and transceiver geometry parameters, the field distortion caused by engineering structures can be explicitly limited during forward and inversion processes. Furthermore, signal-to-noise ratio weighting reduces the impact of anomalous frequencies and late-stage low-quality time windows. The resistivity and polarizability distributions obtained from the joint inversion are further combined with the polarization attenuation half-decay time to distinguish different types of low-resistivity anomalies and extract three-dimensional boundaries, thereby forming the basis for adjusting drilling speed, rotational torque, and tunneling direction, achieving a data closed loop between detection results and underground construction operations.

[0025] Other technical features and functions of the present invention will be further explained in the following specific embodiments in conjunction with the accompanying drawings.

[0026] The following examples are used to illustrate the technical solutions of the present invention, and are not intended to limit the scope of protection of the present invention. Attached Figure Description

[0027] The accompanying drawings are used to illustrate the device layout, spatial relationships, data acquisition, and inversion effects in the embodiments of the present invention; the same or corresponding objects in the drawings are identified by the same identifier.

[0028] Figure 1 A schematic diagram of the parallel arrangement of the positive and negative electrode groups of the underground emission source in the underground roadway;

[0029] Figure 2 This is a schematic cross-sectional view of the underground composite excitation, stratigraphic propagation, and surface array reception.

[0030] Figure 3 This is a schematic diagram showing the spatial relationship between the surface observation area, the topographic surface, and the A and B ends of the underground emission source.

[0031] Figure 4 Schematic diagram of the geological structure and apparent resistivity parameters at different depths in the target mining area;

[0032] Figure 5 A schematic diagram showing the relationship between the surface array grid, the downhole secondary field response, and the synchronous time synchronization.

[0033] Figure 6 The figures show a comparison of the electric field distribution and resistivity inversion results of ground-excited and well-excited methods. Figures (a) to (d) represent the electric field distribution under different excitation methods and frequencies, while figures (e) to (g) represent the depth curves of the actual model and different inversion results. Detailed Implementation

[0034] The embodiments of the present invention will be described below with reference to the accompanying drawings. These embodiments illustrate the technical relationship between underground source deployment, time-frequency observation, engineering constraint inversion, and operational parameter generation; without altering this technical relationship, each step can be parameterized according to the site conditions of the mine roadway, tunnel, or other underground engineering project.

[0035] The term "well-to-surface cross-medium electromagnetic detection" as used in this application refers to a detection process in which a controllable emission source is set up in an underground tunnel or underground working space, allowing the electromagnetic field to propagate through the surrounding rock and its internal target geological body and be collected by a surface array. The "frequency domain stable section" refers to the stable energized section used to extract the electric field amplitude and phase, and the "time domain turn-off section" refers to the waveform section used to establish rapid turn-off conditions and observe the attenuation of the secondary field. The "engineering spatial coupling constraint model" is a forward model and boundary constraints jointly formed by the tunnel cavity, surrounding rock interface, metal structures, topography, and transceiver geometry.

[0036] In this embodiment, satellite timing or other airborne communication nodes are used to provide time reference for both the transmitter and receiver, and are not used as an independent propagation medium for electromagnetic detection. The underground transmitter, the surface array, and the ground data processing center together constitute the well-to-surface detection link; the primary field generated by underground excitation propagates in the surrounding rock and forms a secondary field at the target geological body with electrical differences, and the surface array collects the superposition response of the primary and secondary fields.

[0037] The method in this embodiment sequentially performs the following steps: underground source deployment and composite excitation, surface array time-frequency acquisition, observation data preprocessing, engineering spatial constraint decoupling and joint inversion, geological body identification, and operation parameter generation. The output of the previous stage constitutes the input of the next stage, and the correspondence is maintained by a unified transmit / receive time base, spatial coordinate reference, and data field identifier.

[0038] First, the location of the underground emission source is determined. Based on existing geological data of the target mining area, the three-dimensional coordinates of the underground tunnels, and the apparent resistivity of the strata measured on-site, a tunnel section that meets the emission current requirements and is easy to isolate from metal structures is selected. In one example, a tunnel section with a depth of not less than 600 meters is selected, and mudstone, sandy mudstone, argillaceous sandstone, or karst fracture zone with an apparent resistivity of less than 500 ohm-meters is selected as the grounding electrode deployment area. The above depth and resistivity are only the on-site configuration values ​​for this embodiment, and can be re-determined according to the target depth and power supply capacity in actual engineering.

[0039] Within the selected section, long-conductor grounding electrodes are preferentially used as the transmitting end, and multiple grounding electrodes are connected in parallel to form positive and negative electrode groups. First, the independent grounding resistance of each candidate electrode point is measured. Then, the allowable equivalent total grounding resistance is determined based on the transmitting voltage and the target transmitting current to avoid directly equating the formation resistivity with the electrode grounding resistance. The allowable equivalent total grounding resistance is determined according to the following formula:

[0040] ;

[0041] in, This indicates the supply current, measured in amperes. This indicates the power supply voltage to the transmitter, measured in volts. This represents the equivalent total grounding resistance required by the system, expressed in ohms.

[0042] After obtaining the permissible equivalent total grounding resistance, calculate the equivalent grounding resistance of the parallel electrode group based on the measured independent grounding resistance of each electrode:

[0043] ;

[0044] in, This represents the equivalent total grounding resistance after multiple grounding electrodes are connected in parallel, and the unit is ohms. The total number of grounding electrodes connected in parallel is a positive integer greater than or equal to 2. The index variable representing the electrode (with values...) ); Indicates the first The independent grounding resistance of each grounding electrode (assuming that the grounding resistance of each electrode is equal and all are...) When, the formula can be simplified to ), the unit is ohms.

[0045] Accordingly, the number of grounding electrodes is gradually increased until the equivalent grounding resistance after parallel connection meets the transmission requirements. For example, when the transmission voltage is 1500 volts and the target transmission current is not less than 30 amperes, the allowable total equivalent grounding resistance is not greater than 50 ohms; when the measured independent grounding resistance of a single electrode point is about 500 ohms, ten or more parallel electrode points can be set. Figure 1 The diagram illustrates a layout where the positive (+A) and negative (-B) terminals are each composed of multiple electrodes connected in parallel. These electrodes can be vertically installed on the tunnel floor via electrode holes. If the floor is hardened, the measured independent grounding resistance is too high, or the grounding conditions cannot be met on-site, the system switches to an ungrounded return line excitation mode.

[0046] After the physical deployment of the transmitter is completed, a composite waveform current is injected into the transmitter. This composite waveform current is configured as a variable duty cycle waveform sequence consisting of a frequency-domain stabilization segment and a time-domain shutdown segment. Within a single transmission cycle, this waveform sequence includes a frequency-domain stabilization segment at the front end and a time-domain shutdown segment at the rear end. The output frequency of the frequency-domain stabilization segment is configured to match the near-field frequency-domain depth sounding resolution requirements; the shutdown time threshold of the time-domain shutdown segment is set to meet the requirements of the far-field time-domain diffusion depth. In the specific parameter configuration, the shutdown time of the time-domain shutdown segment is set to less than 20 microseconds. This critical value is selected based on the skin effect physical mechanism. Before operation, the initial apparent resistivity of the shallow surrounding rock in the test area is measured, and combined with the expected advanced detection depth, the critical time constraint is derived based on the early time-domain diffusion formula. The set shutdown time... (Unit: seconds) must meet the following requirements (in The value is the permeability of free space, expressed in Henry per meter. The depth of the target detection is expressed in meters. The value represents the initial apparent resistivity of the survey area (in ohms-meters). This ensures that the secondary field attenuation process can generate a sufficiently strong eddy current response at the target depth before the high-frequency signal is cut off, avoiding premature dissipation of high-frequency characteristic information in shallow strata or attenuation of deep detection signal strength due to excessively long cutoff time.

[0047] After the composite excitation current enters the transmitting end, a primary electromagnetic field is established under the boundary conditions formed by the interface between the tunnel cavity and the surrounding rock, and a secondary field is induced in the deep electrical anomaly body. The geometry of the tunnel and the electrical differences of the surrounding rock can guide or concentrate the field energy distribution. In this embodiment, this effect is regarded as a waveguide enhancement factor in the propagation process from the underground source to the surface array, rather than as a separate inversion conclusion. Figure 2 The diagram shows the cross-sectional relationship between the launch terminal, the target geological body, and the surface array at a depth of approximately 600 meters underground.

[0048] To unify the description of frequency domain and time domain observations, a propagation model of an underground electric dipole source in a homogeneous or layered medium is used as the basis for forward modeling. In subsequent calculations, an engineering spatial coupling constraint model is used to correct for tunnel cavities, metallic bodies, and terrain boundaries. The horizontal electric field in the frequency domain can be expressed as:

[0049] ;

[0050] ;

[0051] in, It represents the horizontal electric field strength in the frequency domain, that is, the apparent resistivity response function of the horizontal electric field across the entire frequency range, with units of volts per meter; This indicates the amplitude of the supply current intensity, measured in amperes. The effective physical length of the emitting electric dipole is expressed in meters. It represents the apparent resistivity across the entire frequency range, and the unit is ohm-meter; The angle between the direction of the transmitted electric dipole and the line connecting the center point of the electric dipole to the receiving point on the ground surface is expressed in radians. This indicates the transmission distance between the underground transmitter and the surface receiver, expressed in meters. Represents the wave number, expressed in meters to the power of negative one; Represents the imaginary unit; It represents the angular frequency of the signal, with the unit being radians per second; This represents the magnetic permeability of the formation medium, measured in Henry per meter.

[0052] The magnetic field in the frequency domain perpendicular direction can be represented as:

[0053] ;

[0054] in, It represents the magnetic field strength in the vertical direction in the frequency domain, that is, the apparent resistivity response function of the vertical magnetic field across the entire frequency range, and the unit is ampere / meter; This indicates the amplitude of the supply current intensity, measured in amperes. The effective physical length of the emitting electric dipole is expressed in meters. The angle between the direction of the transmitted electric dipole and the line connecting the center point of the electric dipole to the receiving point on the ground surface is expressed in radians. This indicates the transmission distance between the underground transmitter and the surface receiver, expressed in meters. The wave number is represented and its calculation method is the same as the aforementioned formula.

[0055] The electric field in the horizontal direction in the time domain can be expressed as:

[0056] ; ;

[0057] in, It represents the horizontal electric field strength in the time domain, that is, the apparent resistivity response function of the horizontal electric field over the entire time period, with units of volts per meter; This indicates the amplitude of the supply current intensity, measured in amperes. The effective physical length of the emitting electric dipole is expressed in meters. The free permeability is represented by a fixed value. The unit is Newtons per square ampere or Henry per meter; This indicates the transmission distance between the underground transmitter and the surface receiver, expressed in meters. The angle between the direction of the transmitting electric dipole and the line connecting the center point of the electric dipole to the receiving point is expressed in radians. Represents the integral variable; Represents a dimensionless time-depth parameter; This represents the apparent resistivity over the entire time period, and the unit is ohm-meter. This indicates the observation time, expressed in seconds.

[0058] The magnetic flux density in the perpendicular direction in the time domain can be expressed as:

[0059] ;

[0060] in, It represents the magnetic flux density in the vertical direction in the time domain, that is, the apparent resistivity response function of the vertical magnetic field over the entire time period, and the unit is Tesla. This indicates the amplitude of the supply current intensity, measured in amperes. The effective physical length of the emitting electric dipole is expressed in meters. The angle between the direction of the transmitting electric dipole and the line connecting the center point of the electric dipole to the receiving point is expressed in radians. This represents the permeability of free space, expressed in Henry per meter. Indicates the transmit / receive distance, in meters; Represents the integral variable; This represents a dimensionless time-depth parameter.

[0061] After the underground emission source is deployed, the observation area is determined based on the projection position of the underground tunnel onto the surface and the results of on-site background noise tests. For example... Figure 3 As shown, multiple five-channel acquisition stations can be deployed in a distributed grid within a range of approximately 2 to 5 kilometers from the projection location. The station spacing and boundaries can also be adjusted based on the effective signal coverage area. Combined with... Figure 4 The stratigraphic structure and apparent resistivity profile shown are used to set the key observation directions; each acquisition station is equipped with mutually orthogonal horizontal electric field sensors and three-dimensional magnetic field sensors.

[0062] Before data acquisition, activate the synchronization controller between the underground transmitter and the surface acquisition station. For example... Figure 5 As shown, a common time axis is established using satellite timing signals or a high-stability crystal oscillator, and the on, stabilization, and off times of the transmitted waveform are correlated with the sampling windows of each acquisition station. During the frequency-domain stabilization period, the electric field amplitude and phase are recorded to form frequency-domain electric field data; after the time-domain off period, the amplitude and attenuation slope of the secondary magnetic field are recorded to form time-domain magnetic field data. Synchronization errors are controlled within the system calibration allowable range, and timestamps and transmission cycle numbers are recorded in the data file.

[0063] After obtaining the raw time-frequency data from the surface array, observation quality control and geometric normalization are first performed, followed by the establishment of an engineering spatial coupling constraint model and the implementation of joint inversion. This processing sequence avoids directly incorporating significant instrument saturation, terrain differences, or device coefficient differences into the three-dimensional inversion.

[0064] Before inputting the data into the inversion model, the collected frequency-domain electric field data and time-domain magnetic field data are first subjected to multi-level denoising preprocessing. The preprocessing process includes four sub-stages. The first sub-stage is the removal of gross errors from field observations. The time-domain waveforms and frequency-domain spectra of each measuring point are thoroughly checked. By setting waveform amplitude thresholds, distorted waveforms that exceed the instrument's dynamic range and cause saturation clipping are removed, as are abnormal measuring point data whose baseline drift amplitude exceeds the set standard. The baseline mean of the received signal is extracted during the early background period when the transmit current has no field source. If the peak-to-peak value of the transient drift during this period exceeds 15% of the maximum amplitude of the normal received signal, it is determined that the front-end amplifier circuit of the sensor in that channel is in a transient saturation or unlocking state caused by a strong near-source electromagnetic pulse. The entire measuring point data is then marked as a discarded channel and removed. Further, the statistical method of 3 The criteria (where σ represents the standard deviation of multiple samples from the same measuring point, with the same unit as the amplitude unit of the corresponding sampled data) and the median absolute deviation algorithm are used to calculate the distribution characteristics of multiple sampled data, and single superimposed abrupt outliers exceeding the statistical threshold are removed. For abandoned tunnel data that are strongly affected by residual metal in the goaf or large metal structures in the local tunnels, resulting in strong human interference, ordinary kriging interpolation is performed using data from spatially adjacent normal measuring points, with a constant spatial semivariogram using a spherical model. The range parameter of the spherical model is set to the straight-line spatial distance from the current abandoned measuring point to the farthest effective reference measuring point participating in the interpolation; the sill parameter is set to the statistical variance of the electromagnetic amplitudes of all adjacent normal measuring points participating in the interpolation. After gross error removal, data from the early period of signal reception without field source background are extracted, the average electromagnetic baseline value is calculated, and the entire data is subtracted from this average baseline value to complete baseline zeroing and eliminate the instrument's zero-point drift error. A polynomial fitting algorithm with the highest power fixed at 2 or 3 was used to extract the slowly varying trend term from the data and subtract it from the original data to remove low-frequency error signals caused by temperature drift. Simultaneously, data from the inclinometer built into the sensor was read, and the amplitudes of each orthogonal component of the electromagnetic field were corrected by three-dimensional spatial rotation projection based on the measured inclinometer deviation values.

[0065] The second sub-stage involves superposition and frequency band noise reduction. Repeated observations at the same measurement point are superimposed and averaged. For time-domain magnetic field data, moving average, median filtering, wavelet thresholding, or empirical mode decomposition are used sequentially to separate the low-frequency effective response. For frequency-domain electric field data, bandpass and power frequency comb notch filtering are applied. The signal-to-noise ratio (SNR) for each frequency point and time window is calculated based on the effective signal amplitude and the root mean square value of the background noise. Data with an SNR below 0 dB are set to zero in the covariance matrix; data between 0 dB and 3 dB are weighted according to the square of the ratio of the measured SNR to 3 dB; and the remaining data retain normal weights.

[0066] The third sub-stage is topographic and elevation correction. The digital elevation model of the survey area is read, the surface undulation boundary is established, and the spatial field is extended or equivalently corrected using a unified reference elevation surface, so that the observations of different elevation acquisition stations are converted to the same elevation datum; for irregular scattered points, two-dimensional spline interpolation can be used to resample into regular grid or equally spaced profile data.

[0067] The fourth sub-stage is device and geometric correction. Based on the actual transmit / receive distance, receiving electrode distance, transmit current, and transmit dipole moment, the geometric factor and device constant are calculated to normalize the observed field values ​​to uniform transmit current, magnetic moment, and pole distance conditions. Combining the field component characteristics of the near-field, transition, and far-field regions, the near-field induction components that do not satisfy the plane wave approximation are identified for subsequent constrained inversion.

[0068] Near-source interference is first mitigated through physical deployment. Three-dimensional distance measurements are taken of structures such as tracks, metal beams, and pipelines, prioritizing a spatial distance of 3 to 5 meters between the parallel electrode group and these structures, with a lateral offset of at least 2 meters from the track. When space is limited, the long axis of the electrode group should be parallel to the extension direction of the roadway and main metal structures to avoid strong orthogonal coupling between the electrode group and long, narrow metal bodies. For metal connections that can be temporarily disconnected, insulation can be used to weaken the closed induced eddy current loop, provided it does not compromise construction safety.

[0069] Engineering structures that cannot be eliminated through on-site layout are incorporated into mathematical constraints. By reading 3D mine engineering drawings or laser scanning point clouds, the spatial topology of roadway centerlines, cross-sectional dimensions, cavity boundaries, metal tracks, canopies, anchor bolts, and pipelines is extracted. The target space is divided into a non-uniform hexahedral mesh, and initial physical properties are assigned based on measured or calibrated electrical parameters of the rock and metal materials. The mesh corresponding to the metal structures can be used as a known region with high electrical conductivity for locking, and the surrounding rock interface is treated according to the corresponding boundary continuity conditions, thereby assembling and forming an engineering spatial coupling constraint model.

[0070] During joint inversion, preprocessed frequency-domain electric field data and time-domain magnetic field data are combined into a target observation vector. The engineering spatial coupling constraint model is then input into a three-dimensional finite element forward modeling program to obtain a prediction vector. The inversion objective function consists of weighted data residuals and model smoothing constraints. The data covariance matrix carries the aforementioned signal-to-noise ratio weights, and the prior model vector carries the initial physical properties of the tunnel, metallic body, and known strata. The objective functional can be expressed as:

[0071] ;

[0072] in, For the target functional; The vector of underground three-dimensional spatial physical property parameters to be inverted (including resistivity and polarizability). The target observation vector; The predicted vector for the three-dimensional finite element forward model; This is the data covariance matrix (used for weighted allocation based on signal-to-noise ratio). This is a regularization parameter used to balance data fit and model smoothness; The roughness operator matrix of the model; This represents the prior initial model vector (i.e., the prior physical properties extracted from the engineering space coupling constraint model). The model roughness operator matrix is ​​described below. It is constructed using a first-order forward difference operator in three-dimensional space, specifically along the spatial coordinate system. , , It is composed of difference matrices in three orthogonal directions; the regularization parameter The initial value is dynamically updated using a strategy of nonlinear cooling decay with the number of iterations. It equals the initial data residual term obtained from the first forward modeling (i.e. ) and initial model constraint terms (i.e. The ratio of ) to; the subsequent Second-rate( The regularization parameter update formula for the iteration (where the integer is a positive integer greater than or equal to 1) is as follows: ,in For the first The regularization parameter for the next iteration For the first The regularization parameter for the next iteration, the cooling coefficient This is a dimensionless constant with a fixed value range between 0.5 and 0.8. In the above regularization parameter update process, the symbol k represents the iteration number, which has a different meaning from the symbol k used to represent the wavenumber. The nonlinear conjugate gradient method is used to update the target functional. The derivative is calculated and the direction of steepest gradient descent is found through iterative solutions.

[0073] Each iteration updates the resistivity and polarizability in the 3D grid based on the residual between the observed vector and the predicted vector, until the target functional is below the convergence threshold or the maximum number of iterations is reached, and outputs a two-parameter 3D distribution. Figure 6 A comparison of field distribution and inversion curves for ground-excited and well-ground-excited conditions is presented: the energy distribution under well-ground-excited conditions is more concentrated near the target depth, and the obtained resistivity change interface corresponds more closely to the real model; this figure is used to illustrate the test results of this embodiment and does not constitute an absolute limitation on the effect under all geological conditions.

[0074] After obtaining the two-parameter three-dimensional distribution, low-resistivity anomaly regions were first extracted according to the bedrock resistivity threshold calibrated on-site. The polarizability variation sequence of each grid node within the anomaly region over time was then correlated with spatial coordinates. Resistivity was used to delineate the candidate anomaly range, and polarization attenuation characteristics were used to further distinguish different geological bodies with similar low-resistivity features.

[0075] The polarizability decay sequence in the low-resistivity anomaly region is fitted using a single-exponential relaxation model with nonlinear least squares fitting. The baseline function is:

[0076] ;

[0077] in, Indicates time The fitted value of the changing apparent polarizability is a dimensionless quantity, where t represents the decay observation time calculated from the end of the time domain cutoff segment, in seconds; is the magnitude parameter of the initial polarizability maximum value, and is a dimensionless quantity; The time constant characterizing the dielectric polarization discharge rate, expressed in seconds; The constant bias term representing the background free polarization level is dimensionless. This is used to obtain the time constant in the solution. Then, based on the mathematical relationship of half-life The physical time required for the polarizability to decay from its initial maximum to half its value was calculated, which is the half-life of the excitation polarization effect. This indicates the half-life of the excitation polarization effect, expressed in seconds.

[0078] A geological lithology sample library can be established using calibration data from the local or similar strata. In the exemplary calibration of this embodiment, the half-life interval for water-bearing fault fracture zones is 0.5 to 2 seconds, the half-life interval for carbonaceous non-mineralized rock strata is 2 to 8 seconds, and the half-life interval for sulfide-containing metallic targets is greater than 10 seconds. The above intervals are used to demonstrate the discrimination process. When applied to other survey areas, recalibration should be performed based on rock sample testing, known boreholes, or exposed geological bodies to avoid directly using the example thresholds as universal constants.

[0079] The actual half-life is compared with the sample library intervals, and attribute labels are output: when it falls into the water-bearing fractured zone interval, it is marked as a water-bearing anomaly requiring deceleration and early warning; when it falls into the carbonaceous non-mineralized rock layer interval, it is marked as a non-mineralized low-resistivity rock layer; and when it falls into the sulfide target interval, it is marked as a metallic target. Adjacent meshes with the same attribute labels are aggregated into connected components, and three-dimensional boundary point sets are extracted along the attribute or half-life isosurface.

[0080] The geological assessment results and boundary point set are sent to the control unit of the underground tunneling equipment. For water-bearing fractured zones or non-mineral anomalies that need to be avoided, the control unit generates instructions to reduce the advance speed and rotation torque, and enters the safety envelope and candidate trajectory calculation. For sulfide-bearing targets, the control unit can generate a tracking mining direction based on its boundary and target abundance direction, while still being constrained by equipment capacity and construction safety.

[0081] After obtaining the three-dimensional coordinate boundary set of the target geological body to be avoided, the system uses a three-dimensional ellipsoid parameterized equation and least squares method to fit the initial envelope surface covering the anomaly. This is then multiplied by a safety expansion coefficient greater than 1. (Preferably 1.2 to 1.5), generating a safe envelope convex surface. The equation of the ellipsoid is:

[0082] ;

[0083] in, For each three-dimensional coordinate point on the boundary set, the unit of each coordinate component is meters; The coordinates of the centroid of the target anomaly obtained by fitting using the least squares method are given, with each coordinate component in meters. These are the semi-axis lengths obtained from the fitting along the three mutually orthogonal principal axes, in meters; The dimensionless safety expansion coefficient is set to a value greater than 1 (preferably 1.2 to 1.5); the semi-axis length Multiply by the aforementioned safety expansion coefficient respectively Then, the final semi-axis length is used to generate the safety envelope convex surface. Subsequently, using the current drill bit center three-dimensional coordinates as the origin, a diverging simulated ray set is generated using a spatial polar coordinate system with a step size of 1 degree. During this determination process, the equipment safety avoidance margin threshold is... (Unit: meters) Based on the pre-defined physical dimensions of the underground tunneling equipment, the specific definition formula is as follows: .in, This indicates the maximum physical radius of rotation of the current underground tunneling equipment's cutting head or outer contour, in meters. This represents the preset buffer margin for rock mass cutting and breaking, in meters, and its value is a fixed constant between 0.5 and 1.5.

[0084] For each candidate ray, the intersection state and minimum normal distance with the safety envelope are calculated. Directions that do not meet the equipment safety avoidance margin are eliminated. Then, the direction with the smallest angle with the original tunneling axis is selected from the remaining directions, and the horizontal deflection angle and vertical pitch angle are output. For anomalies that need to be avoided, the drilling trajectory can be adjusted or reinforcement processes such as pre-grouting can be triggered accordingly. For mineable targets, tracking commands can be generated along the spatial extension direction of the target body while meeting the equipment safety boundary, thereby completing the closed loop between detection, discrimination and construction control.

[0085] The system in this embodiment includes an underground composite excitation subsystem, a surface array acquisition subsystem, and an engineering constraint inversion and operation control center. Each component exchanges data using a unified timestamp, measurement point number, coordinate reference, and launch cycle number, and executes the methods and steps described above.

[0086] The underground composite excitation subsystem is installed within the underground mining tunnel and is encased in an explosion-proof shell. It includes a programmable rectifier-inverter unit, a power switch array, a transmitter switching unit, and a high-stability synchronization module. The programmable rectifier-inverter unit forms a frequency-domain stabilization segment and a time-domain shutdown segment according to control commands; in this embodiment, its rated output current capability is not less than 100 amperes. The transmitter switching unit selects the excitation mode between the grounded electrode group and the ungrounded return line. The high-stability synchronization module generates turn-on, stabilization sampling, and turn-off trigger signals based on an external clock or a temperature-controlled crystal oscillator.

[0087] The surface array acquisition subsystem consists of multiple acquisition stations. Each acquisition station includes an electric field preamplifier circuit, a three-dimensional magnetic field sensor, an analog-to-digital converter, a synchronization controller, and a digital signal processing unit. The digital signal processing unit sets a narrowband notch filter for the 50 Hz fundamental wave and its harmonics. In this embodiment, the quality factor of the comb-shaped notch filter is greater than 50, and the -3dB bandwidth of each stopband is less than 1 Hz, in order to reduce industrial power grid interference and preserve adjacent effective frequency bands as much as possible.

[0088] Because industrial power grids in mines or construction areas exhibit time-varying fluctuations, the system incorporates an adaptive pole placement method to dynamically track small-range, unexpected frequency deviations within the range of 49.8 Hz to 50.2 Hz. Specifically, the built-in digital phase-locked loop extracts the current actual frequency of the power grid in real time. The system dynamically updates the pole and zero positions of the Infinite Impulse Response (IIR) adaptive notch filter in the z-domain based on this frequency, and its system transfer function is configured as follows:

[0089] ;

[0090] in, The system transfer function represents a digital filter; Represents the complex variable operator of the z-transform; The normalized digital angular frequency is given by f0 in the formula, which represents the current actual frequency of the power grid extracted in real time by the digital phase-locked loop, in Hertz. The sampling rate for analog-to-digital conversion is expressed in Hertz (Hz). The pole radius parameter is used to control the notch stopband width (quality factor); the symbol 'r' here has a different meaning from the symbol 'r' used to represent the transmit / receive distance, and its value is configured as follows: Through dynamic updates The parameters ensure that while deeply suppressing power frequency spike noise, the energy integrity of the effective secondary field attenuation signal generated by well-to-ground excitation in the extremely low frequency band is preserved to the maximum extent.

[0091] The magnetic field sensor can be equipped with a composite shield, which includes a structural frame and a composite shielding layer formed by stacking high-permeability materials and high-conductivity meshes. The shielding layer utilizes magnetic flux bypass and eddy current loss to attenuate higher-frequency spatial electromagnetic noise. The material, thickness, and opening form of the shielding structure should be calibrated in conjunction with the target's low-frequency detection bandwidth to avoid excessive attenuation of the effective magnetic field attenuation signal.

[0092] The engineering constraint inversion and operation control center can consist of a server or industrial control computer, and is connected to the underground composite excitation subsystem and the surface array acquisition subsystem via fiber optic or wireless links. The center includes a data quality control module, a spatial model construction module, a finite element forward and joint inversion module, a geological body identification module, and an operation parameter generation module. Its outputs include three-dimensional resistivity and polarizability distributions, target geological body attributes, three-dimensional boundary point sets, and drilling speed, rotational torque, and tunneling direction adjustment parameters.

[0093] The aforementioned data processing and control steps can be implemented by a processor executing a program in memory, or by a digital signal processor, a field-programmable gate array, an application-specific integrated circuit, or a combination thereof. The hardware partitioning of each module can be adjusted according to the sampling scale, real-time requirements, and explosion-proof requirements, but the data correspondence between underground excitation, surface acquisition, engineering constraint inversion, and operational parameter generation should be maintained.

[0094] The above embodiments are used to illustrate the present invention. The transmitted waveform, number of sensors, grid size, inversion algorithm, and control interface can be adjusted within the scope defined in the claims.

Claims

1. A method for cross-medium electromagnetic detection and surrounding rock parameter inversion in wells, characterized in that, include: The layout of the underground transmitter is determined based on the location of the underground roadway in the target mining area, the electrical parameters of the target stratum, and the transmission power supply parameters. Grounded electrode groups or ungrounded return lines are laid in the selected stratum, and the underground transmitter outputs a composite excitation current that includes a frequency domain stable segment and a time domain turn-off segment. A unified time reference is established between the underground emission source and the surface array. The surface array collects the frequency domain electric field amplitude and phase data corresponding to the frequency domain stable segment, as well as the time domain magnetic field data after the time domain shutdown segment ends. Weighted observation data are formed based on the signal-to-noise ratio of each frequency point and each time window. The three-dimensional coordinates and cross-sectional dimensions of the underground roadway, the spatial location of the metal structures in the roadway, the topographic elevation of the survey area, and the geometric parameters of the transceiver device are obtained. Gross errors are removed, baselines and trends are corrected, topographic elevations are corrected, and device coefficients are normalized on the weighted observation data. An engineering spatial coupling constraint model is established based on the obtained spatial information. Using the aforementioned engineering spatial coupling constraint model as the initial model and boundary constraints for forward modeling, near-source interference decoupling and joint inversion are performed on the preprocessed frequency domain electric field data and time domain magnetic field data to obtain the three-dimensional resistivity distribution and three-dimensional polarizability distribution of the target mining area. Based on the three-dimensional resistivity distribution, a low-resistivity anomaly region is determined. The polarization decay sequence of the low-resistivity anomaly region is fitted and the polarization decay half-life is determined. Based on the polarization decay half-life, the properties of the target geological body are determined and its three-dimensional boundary is extracted. Then, adjustment parameters for drilling speed, rotation torque and tunneling direction are generated according to the properties and three-dimensional boundary.

2. The well-to-surface cross-medium electromagnetic detection and surrounding rock parameter inversion method according to claim 1, characterized in that, The allowable equivalent total grounding resistance is determined based on the transmission supply voltage and the target transmission current, and the number of grounding electrodes is determined and connected based on the measured independent grounding resistance of each grounding electrode. The grounding electrode group is arranged along the axial direction of the underground roadway and maintains a preset distance from the metal structures in the roadway. When the grounding conditions do not meet the requirements of the equivalent total grounding resistance, the underground transmission source is switched to an ungrounded return line.

3. The well-to-surface cross-medium electromagnetic detection and surrounding rock parameter inversion method according to claim 1, characterized in that, The composite excitation current is a variable duty cycle sequence, the frequency domain stabilization segment adopts a square wave flat-top segment, and the time domain shutdown segment adopts a bipolar trapezoidal wave falling edge segment; the shutdown time of the time domain shutdown segment is determined based on the target detection depth and the apparent resistivity of the measurement area, and the shutdown time is made less than 20 microseconds.

4. The well-to-surface cross-medium electromagnetic detection and surrounding rock parameter inversion method according to claim 1, characterized in that, The surface array is gridded according to the effective observation range determined by the surface projection position of the underground roadway and the on-site background noise. Each acquisition station includes mutually orthogonal horizontal electric field sensors and three-dimensional magnetic field sensors. The turn-off time of the underground emission source is aligned with the sampling time window of each acquisition station through satellite timing signals or high-stability crystal oscillators.

5. The well-to-surface cross-medium electromagnetic detection and surrounding rock parameter inversion method according to claim 1, characterized in that, When preprocessing frequency domain electric field data and time domain magnetic field data, saturation clipping, abnormal drift and statistical outlier data are removed, repeated observation data are superimposed and averaged, and baseline zeroing, polynomial trend subtraction, power frequency and its harmonic suppression, terrain elevation normalization and transmit / receive geometric factor correction are performed.

6. The well-to-surface cross-medium electromagnetic detection and surrounding rock parameter inversion method according to claim 1, characterized in that, When establishing the engineering spatial coupling constraint model, the tunnel cavity, surrounding rock interface, metal track, metal canopy beam and support components are discretized into a three-dimensional non-uniform grid, and the initial physical properties and boundary conditions are set according to the electrical parameters of the corresponding medium to limit the current density distribution in the forward modeling calculation.

7. The well-to-surface cross-medium electromagnetic detection and surrounding rock parameter inversion method according to claim 1, characterized in that, The joint inversion employs an objective functional that includes data fitting terms and model smoothing constraint terms. The frequency points and time windows are weighted by the data covariance matrix determined based on the signal-to-noise ratio, and the regularization parameters are updated during the iteration process until the objective functional satisfies the convergence condition or reaches the maximum number of iterations.

8. The well-to-surface cross-medium electromagnetic detection and surrounding rock parameter inversion method according to claim 1, characterized in that, A single exponential relaxation model was used to fit the polarization decay sequence of the low-resistivity anomaly region. The polarization decay half-life was calculated based on the time constant obtained from the fitting, and the polarization decay half-life was compared with the half-life intervals of the pre-calibrated water-bearing fracture zone, carbonaceous non-mineralized rock strata, and sulfide-bearing metal target body.

9. The well-to-surface cross-medium electromagnetic detection and surrounding rock parameter inversion method according to claim 1, characterized in that, An ellipsoidal envelope is fitted to the three-dimensional boundary of the target geological body that is determined to need to be avoided, and a safety envelope surface is formed using a safety expansion coefficient greater than 1. Candidate rays are generated from the current drill bit position, and the candidate ray with the smallest deflection angle relative to the original tunneling axis and the distance between it and the safety envelope surface that meets the equipment safety avoidance margin is selected to determine the tunneling direction adjustment parameters.

10. A well-to-surface cross-medium electromagnetic detection and surrounding rock parameter inversion system, characterized in that, include: The underground composite excitation subsystem is used to drive the grounded electrode group or the ungrounded return line to output a composite excitation current containing a frequency domain stable segment and a time domain turn-off segment according to the transmission power supply parameters, and to provide a time reference corresponding to the surface array. The surface array acquisition subsystem is used to acquire frequency domain electric field data and time domain magnetic field data under the time reference, and to evaluate the signal-to-noise ratio at each frequency point and time window. The engineering constraint and joint inversion module is used to establish an engineering spatial coupling constraint model based on the spatial parameters of tunnels, metal structures, terrain elevation and transceiver equipment, and to perform preprocessing, near-source interference decoupling and joint inversion of resistivity and polarizability on the collected data. The geological identification and operation control module is used to determine the properties and three-dimensional boundaries of the target geological body based on resistivity anomalies and polarization decay half-life, and output adjustment parameters for drilling speed, rotation torque and tunneling direction to the downhole tunneling equipment.

Images