A karst cave exploration method based on three-dimensional geoelectricity and fissure water identification

Abstract

The present application relates to the technical field of intelligent sensor, in particular to a karst cave exploration method based on three-dimensional geoelectricity and fissure water identification, comprising the following steps: according to three-dimensional geoelectricity induction array, exciting pulse magnetic field, controlling induction coil to collect secondary field induction electromotive force of stratum feedback, generating transient electromagnetic voltage sequence, according to the transient electromagnetic voltage sequence, extracting frequency component and amplitude information of different time nodes. The present application strips the mixed signal into a series of independent inherent modal function components, can effectively separate background geological noise and key fluid activity signal, finally carries out orthogonal transformation and phase derivation to these components, generates instantaneous frequency amplitude vector, and calculates the determination index capable of quantifying the complexity of energy distribution on different frequency components, compares with the preset threshold value, completes the discrimination of fissure water static or dynamic flow transfer attribute, and provides decision support for engineering risk assessment.

Description

Technical Field

[0001] This invention relates to the field of intelligent sensor technology, and in particular to a method for karst cavity exploration based on three-dimensional geoelectricity and fissure water identification. Background Technology

[0002] Intelligent sensors are composite sensor technologies that integrate information acquisition, signal processing, logic judgment, and communication functions. Unlike traditional sensors that only provide analog signals of raw physical or chemical quantities, intelligent sensors have built-in microprocessors, memory, and communication interfaces, enabling them to automatically process the raw data acquired by their sensing units.

[0003] The shortcomings of existing technologies lie in the fact that while smart sensor technology integrates signal processing and logical judgment functions, its core operating mode is based on processing the direct physical quantities collected by its own sensing unit. When faced with complex tasks such as underground exploration that require imaging and diagnosing vast, heterogeneous spaces, the processors built into smart sensors primarily perform data correction tasks such as zero-drift subtraction and gain normalization, aiming to improve the signal quality of a single measurement point rather than constructing a global, multi-dimensional physical field model. Therefore, they cannot fundamentally obtain the data needed to construct underground three-dimensional structures. For example, a smart pressure sensor deployed underground can accurately report water pressure changes at its location, but it is difficult to use this data to map the structure of the entire aquifer, let alone determine the dynamic properties of water flow within distant karst conduits. Therefore, improvements are needed. Summary of the Invention

[0004] The purpose of this invention is to overcome the shortcomings of existing technologies and propose a karst cavity exploration method based on three-dimensional geoelectricity and fissure water identification.

[0005] To achieve the above objectives, the present invention adopts the following technical solution: a method for karst cavity exploration based on three-dimensional geoelectricity and fracture water identification, comprising the following steps: Based on the three-dimensional geoelectric induction array, a pulsed magnetic field is excited, and the induction coil is controlled to collect the secondary field induced electromotive force fed back by the stratum, generating a transient electromagnetic voltage sequence. Based on the transient electromagnetic voltage sequence, the frequency components and amplitude information at different time points are extracted to generate a time-frequency complex feature matrix. The instantaneous amplitude and phase data of the target frequency slice are extracted based on the time-frequency complex feature matrix. The frequency and conductivity parameters are substituted into the skin depth calculation formula to generate an instantaneous depth conversion dataset. Based on the instantaneous depth conversion dataset, the energy clusters in the time-frequency domain are mapped to the spatial geometric coordinate system to generate a karst geoelectric layer impedance map. Based on the karst geoelectric layered impedance spectrum, the differential regions are identified and the corresponding local fluctuation signals are extracted. The local fluctuation signals are decomposed to generate a set of intrinsic mode function components. The intrinsic mode function components are subjected to orthogonal transformation and phase differentiation calculation to generate an instantaneous frequency amplitude vector. Based on the instantaneous frequency amplitude vector, a judgment index is calculated and generated. The judgment index is then compared with a preset fluid dynamics noise threshold to generate a water body dynamic stability discrimination result.

[0006] Preferably, the steps for obtaining the time-frequency complex feature matrix are as follows: Based on the three-dimensional geoelectric induction array, the transmission sequence number and pulse time mark are set, the induction coil is controlled to synchronously collect the secondary field induced electromotive force according to the channel, the channel number, timestamp, sampling step and gain level are recorded, the induced electromotive force is zero drift subtraction and gain normalization is performed, samples exceeding the dynamic range threshold are removed, and transient electromagnetic voltage sequence is generated. Based on the transient electromagnetic voltage sequence, configure the center time sequence, window width sequence, sliding step and frequency scan list of the variable-width moving Gaussian window, slide the window on the time axis according to the center time sequence and calculate the spectral amplitude of each window according to the frequency scan list, record the time node, frequency index and spectral amplitude triplet, and generate a set of frequency components and amplitude information at different time nodes. Based on the frequency components and amplitude information set at different time points, the phase markers in the time-frequency joint calculation results are read, and the amplitude information and phase markers are paired one by one according to the time node index and the frequency index. Complex entries are filled with time as the row index and frequency as the column index to form a time-frequency complex feature matrix.

[0007] Preferably, the steps for obtaining the instantaneous depth transformation dataset are as follows: Based on the time-frequency complex feature matrix, a target frequency index and a time node sequence are set. The corresponding column is read according to the target frequency index and the complex entries are parsed one by one according to the time node sequence. The magnitude of the complex entries is calculated as the instantaneous amplitude and the argument of the complex entries is parsed as the phase. The channel marker and frequency value are bound according to the time node sequence and the order is sorted to generate the instantaneous amplitude and phase data of the target frequency slice. Based on the instantaneous amplitude and phase data of the target frequency slice, frequency parameters and conductivity parameters are introduced, magnetic permeability parameters are set and expressed as angular frequency parameters, and the frequency parameters and conductivity parameters are substituted into the skin depth calculation formula to calculate the depth value for each time node and pair it with the instantaneous amplitude and phase to generate an instantaneous depth conversion dataset.

[0008] Preferably, the steps for obtaining the karst geoelectric stratification impedance spectrum are as follows: Based on the instantaneous depth conversion dataset, the time node index and depth value of the energy cluster are read, the three-dimensional coordinates are converted according to the survey line coordinates and station numbers, and the spatial position of the energy cluster is assigned. The amplitude and phase values ​​are aggregated according to the spatial position to form a layered unit, and the impedance interval is calculated and filled with color marks to generate a karst geoelectric layered impedance map.

[0009] Preferably, the step of obtaining the intrinsic mode function component set is as follows: Based on the karst geoelectric layered impedance map, impedance gradient threshold, impedance ratio threshold, connected pixel threshold and depth index range are set. The interlayer impedance difference and neighborhood impedance gradient are calculated pixel by pixel in the karst geoelectric layered impedance map. Candidate pixels that meet the thresholds are marked and connected domain search is performed to lock the difference region. The amplitude sequence is extracted along the survey line number and depth index of the difference region and spliced ​​according to the station order to obtain the local fluctuation signal. Based on the local fluctuation signal, the extreme value detection step size, envelope interpolation accuracy, symmetry ratio threshold, and energy concentration threshold are set. Local maxima and local minima are located by time index and upper and lower envelopes are generated. The envelope average is calculated and approximate single components are successively stripped from the local fluctuation signal until the symmetry ratio threshold and energy concentration threshold are satisfied simultaneously, thereby generating an intrinsic mode function component set.

[0010] Preferably, the step of obtaining the instantaneous frequency amplitude vector is as follows: Based on the intrinsic mode function component group, set the orthogonal constraint threshold, phase expansion step size and sampling interval index, construct mutually orthogonal reference vector groups for the intrinsic mode function component group in sequence and perform projection to obtain the orthogonal component sequence, calculate the phase increment of adjacent sampling points according to the phase expansion result of the orthogonal component sequence and convert the frequency according to the sampling interval index, and assemble the corresponding amplitude label into a vector structure according to the time index to generate the instantaneous frequency amplitude vector.

[0011] Preferably, the step of obtaining the determination index is as follows: Based on the instantaneous frequency amplitude vector, the amplitude sequence and the corresponding frequency sequence are extracted according to the component number. The square of the amplitude of each component is calculated as the energy component. The sum of all energy components is obtained and the energy distribution probability is obtained by dividing each component's energy by the sum. The component number, amplitude, and frequency information are retained to generate the energy distribution probability. The judgment index is calculated based on the energy distribution probability and the frequency parameters in the instantaneous frequency amplitude vector.

[0012] Preferably, the steps for obtaining the hydrodynamic stability determination result are as follows: Based on the judgment index, each sample is numerically compared with the fluid dynamics noise threshold. Samples with judgment indices greater than or equal to the fluid dynamics noise threshold are marked as dynamic flow states, and samples with judgment indices lower than the fluid dynamics noise threshold are marked as static flow states, thus forming the water body dynamic stability discrimination result.

[0013] Compared with the prior art, the advantages and positive effects of the present invention are as follows: This invention generates a raw transient electromagnetic voltage sequence by exciting a pulsed magnetic field and collecting the secondary field-induced electromotive force fed back from the strata. Then, it performs time-frequency joint calculations using a widened moving Gaussian window, transforming the one-dimensional time-domain signal into a complex feature matrix containing both time and frequency dimensions. This process preserves the complete amplitude and phase information of the signal. Furthermore, by introducing frequency and conductivity parameters into the skin depth formula, the abstract time-frequency domain energy cluster is mapped to a spatial geometric coordinate system, generating an intuitive karst geoelectric stratification impedance map. This achieves a visual transformation from electromagnetic response data to underground spatial structure, improving the understanding of anomalies. The positioning accuracy of the spatial location and morphology of the mass is improved. For the differential regions in the map, local fluctuation signals are further extracted and adaptive mode decomposition is performed to separate the mixed signals into a series of independent intrinsic mode function components. This can effectively separate background geological noise from key fluid activity signals. Finally, orthogonal transformation and phase differentiation are performed on these components to generate instantaneous frequency amplitude vectors. A judgment index that can quantify the complexity of energy distribution on different frequency components is calculated. By comparing with a preset threshold, the static or dynamic flow properties of fracture water are judged, providing decision support for engineering risk assessment. Attached Figure Description

[0014] Figure 1 This is a schematic diagram of the steps of the present invention. Detailed Implementation

[0015] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0016] Please see Figure 1 This invention provides a technical solution: a method for karst cavity exploration based on three-dimensional geoelectricity and fracture water identification, comprising the following steps: Based on the three-dimensional geoelectric induction array, a pulsed magnetic field is excited to control the induction coil to collect the secondary field induced electromotive force fed back by the stratum, generating a transient electromagnetic voltage sequence. Based on the transient electromagnetic voltage sequence, the frequency components and amplitude information at different time points are extracted to generate a time-frequency complex feature matrix. The instantaneous amplitude and phase data of the target frequency slice are extracted based on the time-frequency complex feature matrix. The frequency and conductivity parameters are substituted into the skin depth calculation formula to generate an instantaneous depth conversion dataset. Based on the instantaneous depth conversion dataset, the energy clusters in the time-frequency domain are mapped to the spatial geometric coordinate system to generate a karst geoelectric layer impedance map. Based on the karst geoelectric layered impedance spectrum, the differential regions are identified and the corresponding local fluctuation signals are extracted. The local fluctuation signals are decomposed to generate a set of intrinsic mode function components. The intrinsic mode function components are orthogonally transformed and phase derivatives are calculated to generate an instantaneous frequency amplitude vector. Based on the instantaneous frequency amplitude vector, a judgment index is calculated and generated. The judgment index is then compared with a preset fluid dynamics noise threshold to generate a water body dynamic stability judgment result.

[0017] The steps for obtaining the time-frequency complex characteristic matrix are as follows: Based on the three-dimensional geoelectric induction array, the transmission sequence number and pulse time mark are set, the induction coil is controlled to synchronously collect the secondary field induced electromotive force according to the channel, the channel number, timestamp, sampling step and gain level are recorded, the induced electromotive force is zero drift subtraction and gain normalization is performed, samples exceeding the dynamic range threshold are removed, and transient electromagnetic voltage sequence is generated. Based on the transient electromagnetic voltage sequence, configure the center time sequence, window width sequence, sliding step and frequency scan list of the variable-width moving Gaussian window. Slide the window on the time axis according to the center time sequence and calculate the spectral amplitude for each window according to the frequency scan list. Record the time node, frequency index and spectral amplitude triplet to generate a set of frequency components and amplitude information for different time nodes. Based on the frequency components and amplitude information sets at different time points, the phase markers in the time-frequency joint calculation results are read, and the amplitude information and phase markers are paired one by one according to the time node index and the frequency index. Complex entries are filled with time as the row index and frequency as the column index to form a time-frequency complex feature matrix.

[0018] Specifically, based on the three-dimensional geoelectric induction array, the transmission sequence number and pulse time marker are set. First, the waveform parameters of the transmission current need to be defined. The power supply period of the bipolar rectangular wave is set to 40 milliseconds, with the positive pulse duration being 10 milliseconds, the negative pulse duration being 10 milliseconds, and the turn-off time being 10 milliseconds, thus forming a complete time control logic. Through a high-precision crystal oscillator clock synchronization controller, a trigger signal is sent to the transmission bridge, and the precise moments of the pulse rise and fall edges are recorded as pulse time markers. The data acquisition channel of the induction coil is started, and the sampling frequency of the analog-to-digital converter is configured to be 250kHz, i.e., the sampling interval is 4 microseconds. The acquisition window length is set to cover the entire turn-off period and part of the pulse period. The analog voltage signal output by the induction coil is synchronously read in channel order and converted into digital quantity. During the recording process, the channel number, such as the X component, Y component, and Z component marker, is written synchronously. At the same time, the current timestamp and the specific level value of the hardware gain amplifier, such as the gain factor, are recorded. The factor was set to 128. The raw data was then preprocessed, and a zero-drift subtraction operation was performed. The average value of the last 100 sampling points after the pulse turn-off was calculated as the baseline drift. The calculation formula is: ,in The number of background noise sampling points selected. For the first The voltage amplitude at each sampling point is calculated by subtracting the baseline drift from each sample value in the original voltage sequence to complete zero-drift correction. Subsequently, gain normalization is performed, and the corrected voltage value is divided by the recorded gain factor. The system restores the true signal amplitude, then performs outlier removal and sets a dynamic range threshold. This threshold is set based on the range of the analog-to-digital converter, for example, a range of [missing information]. If V, then set the threshold to 95% of the range, i.e., 9.5V, and check each sample point one by one. If the absolute value of the sample... ,in If the voltage value at the current sampling point is an abnormal data point caused by electromagnetic interference or circuit saturation, it is either discarded or replaced by linear interpolation of adjacent points. After completing all preprocessing steps, the data stream is organized to generate a transient electromagnetic voltage sequence.

[0019] Based on the transient electromagnetic voltage sequence, the center time sequence, window width sequence, sliding step size, and frequency scan list of a variable-width moving Gaussian window are configured. First, based on the diffusion characteristics of the transient electromagnetic field, the center time sequence is generated using a logarithmic equal-interval method, and the start time is set. The end time is 0.1 milliseconds. The time density coefficient is 100 milliseconds. Assuming each attenuation factor is 20 points, calculate the first... The center time of each window ,in For window number, The start time, For the density coefficient, a corresponding window width sequence is calculated, such that the window width increases linearly or logarithmically over time to match the attenuation law of the signal frequency. The nth... Width of each window ,in The window width coefficient, set to 0.2, is used to control the balance between time resolution and frequency resolution. The sliding step size is determined to be the difference between adjacent center times. A frequency scan list is constructed, with the frequency range set from 10Hz to 10kHz. Fifty characteristic frequency points are selected according to a logarithmic distribution. Then, a sliding window is applied along the time axis based on the center time series. A Gaussian window function is applied to each selected time segment for weighting. The expression for the Gaussian window function is: ,in For time variables, The current time is the center time of the window. Given the current window width, perform a frequency scan on the weighted signal sequence for each frequency in the frequency list. Perform short-time Fourier transform or Gabor transform to calculate the spectral amplitude. The calculation formula is: ,in The input is a transient electromagnetic voltage sequence. For integration variables, For the above Gaussian window function, For the first in the scan list One frequency value, The imaginary unit is used, and the calculated modulus value is the spectral amplitude at that time and frequency. All windows and all frequency points are traversed sequentially, and the calculation results are indexed and recorded in time order and frequency order. The time node, frequency index and spectral amplitude triplet are recorded to generate a set of frequency components and amplitude information at different time nodes.

[0020] Based on the frequency components and amplitude information sets at different time points, the phase markers in the time-frequency joint calculation results are read. First, each record in the set is traversed to obtain its corresponding time point index. and frequency index Extract the corresponding phase angle from the phase dataset generated synchronously in the preceding signal processing steps. This phase angle represents the phase lag information of the secondary field response relative to the transmitted waveform, and its value range is typically within... to Between, then create a dimension as An empty matrix container, in which The total number of time points. Given the length of the frequency scan list, begin filling the matrix entry by entry, and apply the amplitude information according to Euler's formula. With phase marker Combined into complex form, the calculation formula is: ,in The first complex characteristic matrix of time and frequency is the... Line 1 Plural entries of a column For the records in the set at time nodes and frequency The spectral amplitude below, For the corresponding phase angle, The imaginary unit, and These are the cosine and sine trigonometric functions, respectively. This calculation converts the frequency domain information in polar coordinates into a complex expression in rectangular coordinates, indexed by time nodes. As the row address of the matrix, the frequency index As the column address of the matrix, the calculated Fill in the corresponding positions in the matrix. If data is missing at a certain time or frequency point, use bilinear interpolation to estimate and fill in the missing data based on the values ​​of the four adjacent grid points. After mapping and filling all data points, a time-frequency complex feature matrix is ​​formed.

[0021] The steps to obtain the instantaneous depth transformation dataset are as follows: Based on the time-frequency complex feature matrix, the target frequency index and time node sequence are set. The corresponding column is read according to the target frequency index and the complex entries are parsed one by one according to the time node sequence. The magnitude of the complex entries is calculated as the instantaneous amplitude and the argument of the complex entries is parsed as the phase. The channel marker and frequency value are bound according to the time node sequence and the order is sorted to generate the instantaneous amplitude and phase data of the target frequency slice. Based on the instantaneous amplitude and phase data of the target frequency slice, frequency parameters and conductivity parameters are introduced, and magnetic permeability parameters are set and expressed as angular frequency parameters. The frequency parameters and conductivity parameters are substituted into the skin depth calculation formula to calculate the depth value for each time node and pair it with the instantaneous amplitude and phase to generate an instantaneous depth conversion dataset.

[0022] Specifically, based on the time-frequency complex feature matrix, the target frequency index and time node sequence are set. First, the required frequency range for analysis needs to be determined. Based on experience from previous geological exploration, 10Hz to 1kHz is selected as the effective frequency band. The corresponding column index set in the time-frequency complex feature matrix is ​​then located as the target frequency index. Simultaneously, a time node sequence covering the entire sampling period is generated, for example, containing 1024 time points from 0 to 100 milliseconds. The matrix data is read column by column according to the set target frequency index. For each read frequency column, the complex entries are parsed row by row according to the time node sequence. Each complex entry is denoted as... ,in For the real part, For the imaginary part, use the formula for calculating the modulus. Extract the instantaneous amplitude at that frequency at that moment. At the same time, the argument calculation formula is used. Phase obtained by analysis ,in The arctangent function is used. If the calculated phase jumps, phase expansion is performed to ensure phase continuity. Then, each parsed amplitude and phase data packet is bound to the corresponding channel tag, such as the X component, Y component, or Z component, and associated with the specific frequency value. All processed data are rearranged and organized according to the chronological order of the time nodes to ensure that the data time sequence is monotonically increasing. Invalid or redundant header and tail data are removed. Finally, a structured data table is built in memory to generate the instantaneous amplitude and phase data of the target frequency slice.

[0023] Based on the instantaneous amplitude and phase data of the target frequency slice, frequency parameters and conductivity parameters are introduced. First, the average conductivity parameter of the formation is set. This parameter is set based on the resistivity test results of rock samples from the test area, for example, a value of 0.02 S / m. Simultaneously, the vacuum permeability parameter is also set. for H / m, reads the frequency value corresponding to the currently processed frequency slice. Convert it into angular frequency parameters The calculation formula is: ,in Pi The frequency is then used as the input, and the frequency and conductivity parameters are substituted into the skin depth calculation formula, which is as follows: ,in The skin depth is the penetration depth of an electromagnetic wave at that frequency. The conductivity of the formation, Permeability, Given the angular frequency, this formula is used to calculate the corresponding depth value for each time node in the current frequency slice. If a time-varying conductivity model is used, the apparent conductivity of each time node needs to be read. Substitute the values ​​into the formula and perform point-by-point calculations. Match the calculated depth values ​​with the original instantaneous amplitude and phase data one by one. Construct a record row containing five dimensions of information: time, depth, amplitude, phase, and frequency. Repeat the above calculation process for all target frequency slices. Summarize all calculation results and sort them by depth to generate an instantaneous depth conversion dataset.

[0024] The steps for obtaining the geoelectric impedance spectrum of karst formation are as follows: Based on the instantaneous depth conversion dataset, the time node index and depth value of the energy cluster are read, the three-dimensional coordinates are converted according to the survey line coordinates and station numbers, and the spatial location of the energy cluster is assigned. The amplitude and phase values ​​are aggregated according to the spatial location to form layered units, and the impedance interval is calculated and filled with color marks to generate a karst geoelectric layered impedance map.

[0025] Specifically, based on the instantaneous depth conversion dataset, the time node index and depth value of the energy cluster are read. First, statistical analysis is performed on the instantaneous amplitude in the dataset, and the average amplitude is calculated. and standard deviation Set energy cluster recognition threshold Iterate through the dataset and find values ​​with amplitude greater than 100%. The data points are marked as candidate energy clusters. The time node index and corresponding depth value of these candidate points are read, and combined with the geographic coordinate information of the survey line, including the latitude and longitude of the survey line's starting point. and survey line azimuth and the station number of each measuring point. and spacing The three-dimensional spatial coordinates of each energy cluster are calculated using the coordinate transformation formula, which is as follows: , , ,in For planar coordinates, This refers to altitude and depth. The coordinates of the starting point of the survey line are in the plane. This refers to the surface elevation. For the calculated depth value, To determine the azimuth of the survey line, the energy clusters are mapped onto a three-dimensional grid based on their calculated spatial locations. The amplitude and phase values ​​within the same grid cell are aggregated, and the average wave impedance of the grid cell is calculated. A color mapping table is set according to the calculated impedance range, for example, mapping low impedance regions to red and high impedance regions to blue. Color filling and rendering are performed on each layered unit to generate a karst geoelectric layered impedance map.

[0026] The steps for obtaining the intrinsic mode function component set are as follows: Based on the karst geoelectric layered impedance map, impedance gradient threshold, impedance ratio threshold, connected pixel threshold and depth index range are set. The interlayer impedance difference and neighborhood impedance gradient are calculated pixel by pixel in the karst geoelectric layered impedance map. Candidate pixels that meet the threshold are marked and connected domain search is performed to lock the difference region. The amplitude sequence is extracted along the survey line number and depth index of the difference region and spliced ​​according to the station order to obtain the local fluctuation signal. Based on the local fluctuation signal, the extreme value detection step size, envelope interpolation accuracy, symmetry ratio threshold, and energy concentration threshold are set. Local maxima and local minima are located by time index and upper and lower envelopes are generated. The envelope average is calculated and approximate single components are successively stripped from the local fluctuation signal until the symmetry ratio threshold and energy concentration threshold are satisfied simultaneously, generating a set of intrinsic mode function components.

[0027] Specifically, based on the karst geoelectric impedance stratification map, impedance gradient threshold, impedance ratio threshold, connected pixel threshold, and depth index range are set. First, parameters are set according to the electrical differences between karst cavities and surrounding rocks. Then, the statistical distribution of impedance values ​​of all pixels in the entire map is calculated to obtain the average impedance value. and impedance standard deviation Set the impedance gradient threshold ,in Total number of pixels Set an impedance ratio threshold for the local gradient value. A value of 2.5 indicates that the impedance of the abnormal region must be at least 2.5 times the background impedance, setting a connected pixel threshold. 10 pixels are used to filter out isolated noise points; the depth index range is set. This involves layers 20 to 180 in the vertical direction of the atlas to focus on the effective exploration depth, followed by each pixel within the atlas. Perform the calculation, using the Sobel operator or the first-order difference formula to calculate the interlayer impedance difference. And the impedance difference in the horizontal direction, which in turn synthesizes the neighborhood impedance gradient. ,in coordinates The impedance value at that point is calculated, and the ratio of that impedance to the background impedance is also calculated. The gradient value is then calculated. Greater than And the impedance ratio is greater than Pixels with certain values ​​are marked as "1", and the rest are marked as "0", forming a binary mask image. A connected component search algorithm is then performed on the mask image to check the number of pixels contained in each connected region. ,reserve The region is identified and its geometric center and boundary range are defined as the difference region. Finally, based on the coordinate range of the locked region in the map, the original dataset is traced back to the original dataset. The corresponding original instantaneous amplitude sequence is extracted along the direction of the survey line number and the direction of the depth index. The data segments of all stations belonging to the difference region are spliced ​​together end to end in spatial order to construct a long one-dimensional time series data and obtain the local fluctuation signal.

[0028] Based on the local fluctuation signal, the extreme value detection step size, envelope interpolation accuracy, symmetry ratio threshold, and energy concentration threshold are set. First, the decision window size for local extreme value points, i.e., the extreme value detection step size, is defined. This means that the current point must be greater than or less than the three points before and after it to be considered an extreme value. The envelope interpolation precision is set to 0.1 milliseconds, which is the time interval between interpolation points. A symmetry ratio threshold is also set. It is used to assess the degree to which the mean of the upper and lower envelopes deviates from the zero line and to set the energy concentration threshold. It is used to determine whether a component contains the main energy, and to evaluate the local fluctuation signal of the input. Perform a scan to identify all sequences of local maxima. and local minimum point sequence By using a cubic spline interpolation function that passes through all maxima and minima, an upper envelope is constructed. and lower envelope Calculate the average curve of the upper and lower envelopes. The screening process is performed, and the average value curve is subtracted from the original signal to obtain the intermediate signal. Calculate the symmetry index of the intermediate signal. ,like Less than and The energy percentage exceeds Then determine It is an intrinsic mode function component; otherwise, it will be The envelope fitting and stripping steps are repeated with the new input signal. After each component is separated, the original signal is subtracted from the separated component to obtain the remaining residual signal. The residual signal is then decomposed in the next round until the residual signal is a monotonic function or its amplitude is lower than the preset noise level. All separated components are then decomposed. Arrange them from high to low frequency to generate a set of intrinsic mode function components.

[0029] The steps to obtain the instantaneous frequency amplitude vector are as follows: Based on the intrinsic mode function component set, set the orthogonal constraint threshold, phase expansion step size and sampling interval index, construct mutually orthogonal reference vector sets for the intrinsic mode function component set in sequence and perform projection to obtain the orthogonal component sequence. Calculate the phase increment of adjacent sampling points according to the phase expansion result of the orthogonal component sequence and convert the frequency according to the sampling interval index. Combine the corresponding amplitude label with the time index to assemble into a vector structure to generate the instantaneous frequency amplitude vector.

[0030] Specifically, based on the intrinsic mode function component set, the orthogonality constraint threshold, phase expansion step size, and sampling interval index are set. First, the orthogonality evaluation index is set, and then the orthogonality constraint threshold is set. To ensure that the constructed conjugate pairs meet the orthogonality requirement, the phase expansion step size is set to... Set the physical time interval corresponding to the sampling interval index. For example, 4 microseconds, for each extracted intrinsic mode function component The orthogonal reference vector is constructed using the Hilbert transform. The calculation formula is: ,in Let it be the integral variable, and verify it. Is the correlation coefficient less than To satisfy the orthogonality constraint, an analytic signal is constructed using the obtained orthogonal pairs. ,in Using the imaginary unit, calculate the instantaneous phase of the analytic signal. Perform an unfolding operation on the phase sequence to detect the phase difference between adjacent time points. ,like Then the subsequent phase addition ,like Then reduce , to obtain continuous phase The instantaneous frequency is calculated using the derivative of the phase with respect to time, as shown in the formula: Simultaneously calculate the instantaneous amplitude. Frequency values ​​at each time point and amplitude By component number and time index By combining these features, a multidimensional feature vector containing frequency and amplitude is constructed, generating an instantaneous frequency amplitude vector.

[0031] The steps for obtaining the judgment indicators are as follows: Based on the instantaneous frequency amplitude vector, the amplitude sequence and the corresponding frequency sequence are extracted according to the component number. The square of the amplitude of each component is calculated as the energy component. The sum of all energy components is obtained and the energy distribution probability is obtained by dividing each component's energy by the sum. The component number, amplitude, and frequency information are retained to generate the energy distribution probability. Based on the energy distribution probability and the frequency parameters in the instantaneous frequency amplitude vector, the judgment index is calculated using the following formula: ; in, The criterion for judging the current time sample. For the first The probability distribution of energy for each component is calculated by dividing the energy of each component by the total energy. For the first The instantaneous frequency of each component is obtained from the instantaneous frequency amplitude vector. This is the sum of the instantaneous frequencies of all components. These are frequency modulation coefficients, used to adjust the influence of high-frequency components in entropy calculations. The total number of components, and These are the component serial numbers.

[0032] Specifically, based on the instantaneous frequency amplitude vector, the amplitude sequence and corresponding frequency sequence are extracted according to the component index. First, the multi-dimensional data structure inside the vector is parsed. This structure is usually indexed by time step, and each time step contains several intrinsic mode components obtained from the decomposition. Each time sampling point is traversed, and for each component index at that point, the sequence is processed from the first component to the... For each component, read its stored instantaneous amplitude and instantaneous frequency values ​​one by one, and record the read amplitude data as follows: Using the principle of signal energy analysis, the instantaneous amplitude of each component is squared, and this value is taken as the instantaneous energy component of that component at the current moment. For example, if the amplitude of a component at the current moment is 4 microvolts, then its energy component is calculated to be 16. After completing the energy calculation of all components at that moment, an accumulator is built to collect all the data. The energy values ​​of each component are added together to obtain the total energy value of the signal at the current moment. Then, a normalization calculation is performed, dividing the energy value of each component by the total energy value at that moment to obtain the relative proportion of each component in the overall energy. This proportion is the energy distribution probability, which describes the weighted distribution of energy carried by different frequency components. It is ensured that the sum of the probabilities of all calculated components is strictly equal to 1. At the same time, for the convenience of subsequent calculations, the calculated probability value is associated and bound with the original component number, the corresponding instantaneous amplitude, and the instantaneous frequency to construct a feature entry containing complete attributes. Finally, the entries of all processed time sampling points are reorganized in time sequence to generate the energy distribution probability.

[0033] In the calculation formula of the judgment index, the traditional information entropy algorithm is improved by introducing a frequency weighting factor. While evaluating the signal complexity, it also takes into account the distribution characteristics of frequency components. It utilizes the physical law that the electromagnetic noise generated by the flow of fissure water usually has specific high-frequency characteristics. By weighting the frequency proportion, the influence of high-frequency components on the entropy value is amplified, so that it can more sensitively identify the weak high-frequency disturbance signals generated by the water dynamics in karst pipes and distinguish between static geological anomalies and dynamic water anomalies.

[0034] This parameter represents the total number of intrinsic mode function (IMF) components and is a dimensionless integer value. The steps to obtain this parameter are as follows: during the empirical mode decomposition (EMD) or variational mode decomposition (VMD) of the original local geoelectric fluctuation signal, count the number of IMF components that ultimately meet the iteration stopping condition. Typically, depending on the signal complexity, this value ranges from 3 to 10. For example, after decomposing the signal from a karst fissure development zone, if a total of 5 valid IMF components are extracted, then... The value is 5.

[0035] Indicates the first The energy distribution probability of the i-th component is a dimensionless value between 0 and 1. This parameter is obtained by using the energy distribution data calculated in the previous steps to read the i-th component. Energy value of each component Total energy of all components It is obtained through ratio calculation, and the calculation formula is as follows: For example, at a certain moment, the energy of the first component is 20, and the total energy is 100, then... The value is 0.2, which reflects the first... The degree of dominance of each component at the current moment.

[0036] Indicates the first The instantaneous frequency of the nth component, in Hertz (Hz), is obtained by directly extracting the nth component from the input instantaneous frequency amplitude vector. The frequency value of each component at the current time point is obtained by differentiating the phase using the Hilbert transform. For example, if the instantaneous frequency of the first component at the current moment is 50Hz, then... The value is 50.

[0037] This represents the sum of the instantaneous frequencies of all components, measured in Hertz (Hz). The steps to obtain this parameter are: iterate through each component at the current moment, read its instantaneous frequency value, and sum them. The calculation formula is as follows: For example, if there are 3 components with frequencies of 50Hz, 100Hz and 150Hz, the sum is 300Hz.

[0038] The frequency modulation coefficient is a dimensionless empirical constant used to adjust the weight of the frequency term in entropy calculation. The steps to obtain this parameter are as follows: in known water model experiments, electromagnetic response data for both static and dynamic water bodies are collected, and different... The entropy difference under various values ​​is determined by selecting the value that maximizes the distinction between the entropy values ​​of dynamic and static water bodies as the optimal coefficient. For example, through the analysis of 10 sets of simulation experimental data, it was found that when... When set to 1.2, the ability to distinguish between static and dynamic states reaches its highest level; therefore, it is set to... It is 1.2.

[0039] Calculations based on parameters: Decomposed at the current time to obtain Each component, frequency modulation coefficient .

[0040] Component 1: Instantaneous frequency Hz, energy distribution probability .

[0041] Component 2: Instantaneous frequency Hz, energy distribution probability .

[0042] First, calculate the sum of frequencies: Hz.

[0043] Calculate the weighted terms for component 1: .

[0044] Calculate the weighted terms for component 2: .

[0045] Calculate the entropy contribution of component 1: (Approximate calculation with base e, taking the natural logarithm).

[0046] Calculate the entropy contribution of component 2: .

[0047] Calculate the final judgment criteria: .

[0048] The result indicates that the signal at the current moment has a high frequency-weighted entropy value. The value of 0.994 reflects the distribution of signal energy among different frequency components and the activity of high-frequency components. A higher value usually corresponds to the complex flow state of water in karst fissures, which produces electromagnetic induction noise with high chaos and high frequency.

[0049] The steps for obtaining the dynamic stability determination results of the water body are as follows: Based on the judgment index, each sample is numerically compared with the fluid dynamics noise threshold. Samples with judgment indices greater than or equal to the fluid dynamics noise threshold are marked as dynamic flow states, and samples with judgment indices lower than the fluid dynamics noise threshold are marked as static flow states, thus forming the water body dynamic stability discrimination result.

[0050] Specifically, based on the judgment criteria, each sample is numerically compared with the fluid dynamics noise threshold. First, a fluid dynamics noise threshold for distinguishing between dynamic and static states needs to be established. This threshold is set based on statistical analysis of data from known waterless areas or background areas within the survey area. A background geoelectric signal with a length of no less than 500 sampling points is selected, and its judgment criterion sequence is calculated according to the aforementioned steps. The average value of this sequence is then calculated. and standard deviation Set the fluid dynamics noise threshold to Add 3 times For example, if the calculated average value of the background area index is 0.4 and the standard deviation is 0.1, then the threshold is calculated as 0.4 plus 0.3, which equals 0.7. This setting can effectively shield the interference of geological background noise. Next, the state label array is initialized, and the judgment index values ​​of the test samples are read one by one according to the time index order. Each index value is compared with the set threshold of 0.7. If the judgment index of the current sample is greater than or equal to 0.7, it is considered that the signal characteristics at that moment conform to the high entropy characteristics caused by hydrodynamic disturbance, and it is judged that there is groundwater flow and it is marked as a dynamic flow state. If the judgment index is lower than 0.7, it is considered that the signal is at a stable background level and there is no significant water activity, and it is marked as a static flow state. In order to improve the robustness of the judgment, a sliding window voting mechanism is used to smooth the labeling results. For example, the window length is set to 5. Only when more than 3 samples in the window are marked as dynamic flow state is it confirmed as a valid hydrodynamic anomaly. Finally, all the smoothed labels are summarized to form the hydrodynamic stability judgment result.

[0051] The above are merely preferred embodiments of the present invention and are not intended to limit the present invention in any other way. Any person skilled in the art may make changes or modifications to the above-disclosed technical content to create equivalent embodiments that can be applied to other fields. However, any simple modifications, equivalent changes, and modifications made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the protection scope of the present invention.

Claims

1. A method for karst cavity exploration based on three-dimensional geoelectricity and fracture water identification, characterized in that, Includes the following steps: Based on the three-dimensional geoelectric induction array, a pulsed magnetic field is excited, and the induction coil is controlled to collect the secondary field induced electromotive force fed back by the stratum, generating a transient electromagnetic voltage sequence. Based on the transient electromagnetic voltage sequence, the frequency components and amplitude information at different time points are extracted to generate a time-frequency complex feature matrix. The instantaneous amplitude and phase data of the target frequency slice are extracted based on the time-frequency complex feature matrix. The frequency and conductivity parameters are substituted into the skin depth calculation formula to generate an instantaneous depth conversion dataset. Based on the instantaneous depth conversion dataset, the energy clusters in the time-frequency domain are mapped to the spatial geometric coordinate system to generate a karst geoelectric layer impedance map. Based on the karst geoelectric layered impedance spectrum, the differential regions are identified and the corresponding local fluctuation signals are extracted. The local fluctuation signals are decomposed to generate a set of intrinsic mode function components. The intrinsic mode function components are subjected to orthogonal transformation and phase differentiation calculation to generate an instantaneous frequency amplitude vector. Based on the instantaneous frequency amplitude vector, a judgment index is calculated and generated. The judgment index is then compared with a preset fluid dynamics noise threshold to generate a water body dynamic stability discrimination result.

2. The karst cavity exploration method based on three-dimensional geoelectricity and fracture water identification according to claim 1, characterized in that, The steps for obtaining the time-frequency complex feature matrix are as follows: Based on the three-dimensional geoelectric induction array, the transmission sequence number and pulse time mark are set, the induction coil is controlled to synchronously collect the secondary field induced electromotive force according to the channel, the channel number, timestamp, sampling step and gain level are recorded, the induced electromotive force is zero drift subtraction and gain normalization is performed, samples exceeding the dynamic range threshold are removed, and transient electromagnetic voltage sequence is generated. Based on the transient electromagnetic voltage sequence, configure the center time sequence, window width sequence, sliding step and frequency scan list of the variable-width moving Gaussian window, slide the window on the time axis according to the center time sequence and calculate the spectral amplitude of each window according to the frequency scan list, record the time node, frequency index and spectral amplitude triplet, and generate a set of frequency components and amplitude information at different time nodes. Based on the frequency components and amplitude information set at different time points, the phase markers in the time-frequency joint calculation results are read, and the amplitude information and phase markers are paired one by one according to the time node index and the frequency index. Complex entries are filled with time as the row index and frequency as the column index to form a time-frequency complex feature matrix.

3. The karst cavity exploration method based on three-dimensional geoelectricity and fracture water identification according to claim 1, characterized in that, The steps for obtaining the instantaneous depth transformation dataset are as follows: Based on the time-frequency complex feature matrix, a target frequency index and a time node sequence are set. The corresponding column is read according to the target frequency index and the complex entries are parsed one by one according to the time node sequence. The magnitude of the complex entries is calculated as the instantaneous amplitude and the argument of the complex entries is parsed as the phase. The channel marker and frequency value are bound according to the time node sequence and the order is sorted to generate the instantaneous amplitude and phase data of the target frequency slice. Based on the instantaneous amplitude and phase data of the target frequency slice, frequency parameters and conductivity parameters are introduced, magnetic permeability parameters are set and expressed as angular frequency parameters, and the frequency parameters and conductivity parameters are substituted into the skin depth calculation formula to calculate the depth value for each time node and pair it with the instantaneous amplitude and phase to generate an instantaneous depth conversion dataset.

4. The karst cavity exploration method based on three-dimensional geoelectricity and fracture water identification according to claim 1, characterized in that, The steps for obtaining the karst geoelectric stratification impedance spectrum are as follows: Based on the instantaneous depth conversion dataset, the time node index and depth value of the energy cluster are read, the three-dimensional coordinates are converted according to the survey line coordinates and station numbers, and the spatial position of the energy cluster is assigned. The amplitude and phase values ​​are aggregated according to the spatial position to form a layered unit, and the impedance interval is calculated and filled with color marks to generate a karst geoelectric layered impedance map.

5. The karst cavity exploration method based on three-dimensional geoelectricity and fracture water identification according to claim 1, characterized in that, The steps for obtaining the intrinsic mode function component set are as follows: Based on the karst geoelectric layered impedance map, impedance gradient threshold, impedance ratio threshold, connected pixel threshold and depth index range are set. The interlayer impedance difference and neighborhood impedance gradient are calculated pixel by pixel in the karst geoelectric layered impedance map. Candidate pixels that meet the thresholds are marked and connected domain search is performed to lock the difference region. The amplitude sequence is extracted along the survey line number and depth index of the difference region and spliced ​​according to the station order to obtain the local fluctuation signal. Based on the local fluctuation signal, the extreme value detection step size, envelope interpolation accuracy, symmetry ratio threshold, and energy concentration threshold are set. Local maxima and local minima are located by time index and upper and lower envelopes are generated. The envelope average is calculated and approximate single components are successively stripped from the local fluctuation signal until the symmetry ratio threshold and energy concentration threshold are satisfied simultaneously, thereby generating an intrinsic mode function component set.

6. The karst cavity exploration method based on three-dimensional geoelectricity and fracture water identification according to claim 1, characterized in that, The steps for obtaining the instantaneous frequency amplitude vector are as follows: Based on the intrinsic mode function component group, set the orthogonal constraint threshold, phase expansion step size and sampling interval index, construct mutually orthogonal reference vector groups for the intrinsic mode function component group in sequence and perform projection to obtain the orthogonal component sequence, calculate the phase increment of adjacent sampling points according to the phase expansion result of the orthogonal component sequence and convert the frequency according to the sampling interval index, and assemble the corresponding amplitude label into a vector structure according to the time index to generate the instantaneous frequency amplitude vector.

7. The karst cavity exploration method based on three-dimensional geoelectricity and fracture water identification according to claim 1, characterized in that, The steps for obtaining the determination index are as follows: Based on the instantaneous frequency amplitude vector, the amplitude sequence and the corresponding frequency sequence are extracted according to the component number. The square of the amplitude of each component is calculated as the energy component. The sum of all energy components is obtained and the energy distribution probability is obtained by dividing each component's energy by the sum. The component number, amplitude, and frequency information are retained to generate the energy distribution probability. The judgment index is calculated based on the energy distribution probability and the frequency parameters in the instantaneous frequency amplitude vector.

8. The method for karst cavity exploration based on three-dimensional geoelectricity and fracture water identification according to claim 1, characterized in that, The steps for obtaining the hydrodynamic stability discrimination result are as follows: Based on the judgment index, each sample is numerically compared with the fluid dynamics noise threshold. Samples with judgment indices greater than or equal to the fluid dynamics noise threshold are marked as dynamic flow states, and samples with judgment indices lower than the fluid dynamics noise threshold are marked as static flow states, thus forming the water body dynamic stability discrimination result.

Images