An intelligent identification method for geological stress concentration areas
By laying multiple vector electromagnetic induction coils in the monitoring area and using multiple signal processing techniques, multi-dimensional data cubes are processed, geological stress distribution estimates are constructed, and a three-dimensional model of the geological stress concentration area is finally constructed, which solves the accuracy and accuracy of geological stress concentration area identification in the existing technology, and realizes high-precision stress field reconstruction under complex geological conditions.
Patent Information
- Application Number
- CN202510457940.1
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2025-04-14
- Publication Date
- 2025-06-13
- Estimated Expiration
- 2045-04-14
AI Technical Summary
The prior art faces problems of weak signal strength, severe noise interference, difficulty in positioning signal sources and complex data interpretation when identifying geological stress concentration areas, especially in complex geological environments, it is difficult to achieve high-precision stress field reconstruction.
By arranging multiple vector electromagnetic induction coils in the monitoring area, the original electromagnetic pulse signals containing the three direction components of X, Y, and Z are collected, and technical means such as denoising, time-frequency analysis, integer linear programming model and adaptive fuzzy correction algorithm are used to process multi-dimensional data cubes, construct geological stress distribution estimates, and finally build a three-dimensional model of the geological stress concentration area.
The accuracy of stress field boundary recognition is significantly enhanced, the blind area is reduced, high-precision stress field reconstruction under complex geological conditions is realized, the pseudo-anomaly recognition rate is reduced, and the identification accuracy of geological stress concentration areas is improved.
Smart Images

Figure CN119986821B_ABST
Abstract
Description
Technical Field
[0001] This application relates to the field of intelligent recognition technology, and particularly to an intelligent recognition method for geological stress concentration areas. Background Art
[0002] Geological stress concentration areas are regions with a high incidence of geological disasters. Their accurate identification is of great significance for mine safety, earthquake early warning, and geological disaster prevention and control. Traditional methods for identifying geological stress concentration areas mainly rely on borehole stress measurement, acoustic emission monitoring, or microseismic monitoring technologies. These methods not only require a large amount of on-site work and high equipment costs, but also have strong destructiveness and point limitations, making it difficult to achieve large-scale, real-time, and continuous monitoring of the geological stress field. With the development of geophysical exploration technologies, methods based on detecting natural fields have gradually received attention. Among them, natural electromagnetic pulse signals have become one of the important means for identifying geological stress concentration areas due to their advantages of non-destructiveness, real-time nature, and large-scale monitoring.
[0003] Natural electromagnetic pulse signals are electromagnetic phenomena generated by stress changes inside rocks. By analyzing these signals, the stress state of geological bodies can be inferred. However, there are many challenges in existing natural electromagnetic pulse signal monitoring technologies, including weak signal intensity, severe noise interference, difficult signal source positioning, and complex data interpretation. Especially in complex geological environments, traditional signal processing methods are difficult to effectively distinguish useful signals from interference signals, resulting in insufficient accuracy of discrimination results. At the same time, most existing technologies use empirical formulas or simple statistical models for data analysis, lacking a systematic method for accurately identifying geological stress concentration areas. Summary of the Invention
[0004] This application provides an intelligent recognition method for geological stress concentration areas. This application enhances the recognition accuracy of the stress field boundary and reduces the blind area, thereby achieving high-precision stress field reconstruction under complex geological conditions.
[0005] In the first aspect, this application provides an intelligent recognition method for geological stress concentration areas. The intelligent recognition method for geological stress concentration areas includes:
[0006] Arrange a plurality of vector electromagnetic induction coils in a grid pattern in the monitoring area to collect original electromagnetic pulse signals containing X, Y, and Z direction components;
[0007] Perform denoising and time-frequency analysis on the original electromagnetic pulse signals to obtain a multi-dimensional data cube;
[0008] Process the multi-dimensional data cube using an integer linear programming model and an adaptive fuzzy correction algorithm to obtain a corrected electromagnetic pulse signal data set;
[0009] Decompose the corrected electromagnetic pulse signal dataset and construct a mapping relationship to obtain an estimated value of the geological stress distribution;
[0010] Construct a three-dimensional model of the geological stress concentration area based on the estimated value of the geological stress distribution.
[0011] In the technical solution provided by this application, the spatial resolution is improved by a multi-vector electromagnetic induction coil monitoring network with a grid-like distribution. The equilateral triangle arrangement mode and the orthogonal three-component receiving method are adopted to achieve an all-round three-dimensional monitoring of natural electromagnetic pulse signals, significantly enhancing the stress field boundary recognition accuracy and reducing the blind area. The systematic signal preprocessing process includes power frequency filtering, wavelet denoising, and adaptive band-pass filtering. Combined with the time-frequency analysis of the Hilbert-Huang transform, the effective components in weak electromagnetic pulse signals are effectively extracted, and environmental interference and system noise are suppressed. The method of combining an integer linear programming model and an adaptive fuzzy correction algorithm is introduced to solve the problems of signal baseline drift and electromagnetic signal ambiguity that have long troubled this field. The time-reassigned multi-synchronous compression transform breaks through the limitations of traditional Fourier analysis and can effectively process non-stationary signals. Combined with the mapping relationship between signal features and geological stress parameters constructed by the extreme gradient boosting algorithm, high-precision stress field reconstruction under complex geological conditions is achieved. The stress concentration index evaluation system correlates the stress field characteristics with geological structures, and spatial constraint optimization is carried out through the Markov random field model, significantly reducing the false anomaly recognition rate. The three-dimensional model of the geological stress concentration area uses adaptive grid meshing technology to finely depict the spatial geometric characteristics and internal stress distribution of the stress concentration area. Brief Description of the Drawings
[0012] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings required for the description of the embodiments will be briefly introduced below. Obviously, the drawings in the following description are some embodiments of the present invention. For those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative efforts.
[0013] Figure 1 It is a schematic diagram of an embodiment of the intelligent identification method for the geological stress concentration area in the embodiments of this application. Detailed Embodiments
[0014] An embodiment of the present application provides an intelligent identification method for geological stress concentration areas. Terms such as "first", "second", "third", "fourth", etc. (if any) in the specification, claims and the above-mentioned drawings of the present application are used to distinguish similar objects, and do not have to be used to describe a specific order or sequence. It should be understood that the data used in this way can be interchanged under appropriate circumstances, so that the embodiments described here can be implemented in an order other than that illustrated or described here. In addition, the term "comprising" or "having" and any deformation thereof are intended to cover non-exclusive inclusion. For example, a process, method, system, product or device comprising a series of steps or units does not have to be limited to those steps or units clearly listed, but may include other steps or units not clearly listed or inherent to these processes, methods, products or devices.
[0015] For ease of understanding, the specific process of the embodiment of the present application will be described below. Please refer to Figure 1 , an embodiment of the intelligent identification method for geological stress concentration areas in the embodiment of the present application includes:
[0016] Step S101: Arrange a plurality of vector electromagnetic induction coils in a grid pattern in the monitoring area, and collect the original electromagnetic pulse signals including the X, Y, and Z direction components;
[0017] It can be understood that the execution subject of the present application can be an intelligent identification system for geological stress concentration areas, or a terminal or a server. Specifically, it is not limited here. An embodiment of the present application will be described by taking the server as the execution subject as an example.
[0018] Specifically, conduct a geological structure complexity analysis of the monitoring area to obtain key geological parameters such as the distribution characteristics of various structural units, fault density, joint development degree, and lithological differences within the monitoring area. The analysis content includes a comprehensive analysis of existing geological maps, remote sensing image data, and drilling results, supplemented by regional stress field simulation and geometric modeling of the fault system, to obtain an analysis result reflecting the geological structure complexity level. Based on this analysis result, determine which areas within the monitoring area have active geological activities and dense structural intersections, and which areas have simple structures and slow changes. According to the complexity levels of different areas, adopt an adaptive strategy to determine the layout spacing of the receiving devices, and uniformly use the equilateral triangle arrangement pattern as the basic layout form. For areas with higher geological structure complexity, the spacing of the receiving devices is set to 50 meters to improve the spatial resolution; while for areas with relatively simple structures, the spacing is widened to 200 meters to save resources and reduce the layout cost. Through the differential layout strategy, obtain the layout plan of the receiving devices. Install an electromagnetic signal receiving unit containing three orthogonal vector electromagnetic induction coils on each receiving device in the layout plan of the receiving devices. The coils arranged in three axes correspond to the three spatial directions of X, Y, and Z respectively, and simultaneously obtain the electromagnetic pulse vector components in the three directions to ensure the integrity and decoupling ability of the data in the spatial dimension. Each induction coil is connected to a high-sensitivity fluxgate sensor, and the sensor range is ±100 nT to meet the observation requirements of weak natural electromagnetic signals with a large frequency span. At the same time, each electromagnetic signal receiving unit integrates a high-precision GPS clock module, which realizes the unified time calibration between the receiving devices, enabling each station to collect data synchronously at the millisecond level, so as to ensure the time consistency of the responses of different points to the same electromagnetic event. To improve the network flexibility and deployment efficiency of the overall system, each receiving device is embedded with a low-power wireless transmission module, which establishes a communication link with the surrounding devices through this module to construct a mesh network with self-organizing and self-healing capabilities. This distributed mesh network can automatically find an alternative route when any node in the network fails, ensuring the continuity and robustness of data transmission. All electromagnetic signal receiving units form a distributed monitoring network with time synchronization and spatial coordination capabilities. While completing the layout of the main monitoring grid, select several reference sites outside the monitoring area. These sites are far from the main geological structure activity zones, and their main function is to record the electromagnetic interference signals outside the region, such as lightning activities, industrial electromagnetic waves, communication signals, etc. Since the waveforms, frequencies, and amplitudes of interference signals are sometimes similar to those of natural geological electromagnetic pulses, it is easy to cause misidentification. Therefore, independent recording and feature extraction are carried out through the reference sites. The obtained reference interference data will be used as an important control standard for subsequent signal preprocessing. The entire distributed monitoring network continuously collects electromagnetic pulse signals in real time for 24 hours and realizes multi-site synchronous recording through the GPS module.After the data is transmitted to the central processing system, the system automatically compares and analyzes the reference interference data with the electromagnetic signals received at each measuring point, uses algorithms such as waveform matching, frequency band characteristics, and phase difference discrimination to exclude interference signals from non-geological sources, judges the effectiveness of the remaining signals and marks them, and finally outputs the original electromagnetic pulse signals containing the X, Y, and Z direction components.
[0019] Step S102: Perform denoising and time-frequency analysis on the original electromagnetic pulse signals to obtain a multi-dimensional data cube;
[0020] Specifically, the original electromagnetic pulse signal is subjected to interference filtering. Since there is a large amount of background interference in the actual environment caused by power systems, radio communications, lightning activities, etc., which affects the authenticity and recognizability of natural electromagnetic signals, a notch filter is used to suppress the power frequency interference in the signal, especially the 50 Hz and its higher harmonic components. By designing the center frequency and bandwidth of the notch filter, the power system noise is effectively removed, and a preliminary filtered signal is obtained. Then, the wavelet threshold denoising method is applied to the preliminary filtered signal to eliminate random noise and low-frequency drift. The db4 wavelet function in the Daubechies wavelet family is selected as the mother wavelet, and the signal is decomposed into 5 layers by wavelet decomposition. The soft threshold processing method is used to compress the detail coefficients to effectively remove the noise and retain the characteristics of the useful signal. After this processing process, a denoised signal closer to the true response of the geological event is output. The denoised signal is processed by time-domain segmentation to enhance the resolution ability of time-frequency characteristics and adapt to the non-stationary characteristics of geological events. Each signal is divided into time-domain segments with a length of 60 seconds, and a 50% overlap rate is set between adjacent segments, so as to achieve a good balance between time resolution and frequency resolution. The Hilbert-Huang transform is used for time-frequency analysis of each segment during segmentation. This method decomposes the signal into several intrinsic mode functions through empirical mode decomposition, and then performs the Hilbert transform on each mode function to extract its instantaneous frequency, instantaneous amplitude and phase information. In this way, the dynamic characteristics of the frequency changing with time in the natural electromagnetic pulse are effectively captured. The analysis results of all segmented signals are summarized to form a high-dimensional signal feature matrix. Each row corresponds to the analysis result within a time window, and each column represents a specific signal parameter index, such as instantaneous frequency or instantaneous energy, etc. Polarization analysis is performed on the three-component vector electromagnetic signals collected at each receiving point in the obtained signal feature matrix. By calculating the polarization ellipse parameters of the signal within each time period, three key polarization parameters, namely ellipticity, dip angle and azimuth angle, are extracted. These parameters reflect the spatial propagation characteristics of electromagnetic waves and the source mechanism of earthquakes, and help to identify the geological origin of the signal. At the same time, in order to remove potential non-geological source interference signals, a group of adaptive notch filters are constructed in combination with the regional external electromagnetic interference data recorded by the previous reference sites. These filters dynamically attenuate the components in the signal spectrum that are consistent with the characteristic frequencies of the interference sources, and remove non-natural pulse interferences such as lightning, power, and industry to the greatest extent without damaging the structure of the geological effective signal, and obtain pure geological source electromagnetic signals. The pure geological source signals are subjected to spatial correlation analysis. By calculating the signal cross-correlation coefficients between each receiving device and the phase difference within the corresponding time window, the degree of propagation consistency and correlation of the signal in space is evaluated. High correlation and stable phase difference indicate that the signal originates from the same geological stress release event.All processed data are organized in a unified standard format to form a multi-dimensional data cube. Each data unit of this cube contains acquisition time, measuring point location, amplitudes, frequencies, and phase information in the X / Y / Z three directions, constructing a standardized input data volume that describes the spatio-temporal evolution characteristics of the geological electromagnetic field.
[0021] Step S103: Process the multi-dimensional data cube using an integer linear programming model and an adaptive fuzzy correction algorithm to obtain a corrected electromagnetic pulse signal data set;
[0022] Specifically, taking a multi-dimensional data cube as the input, an integer linear programming model with sparsity constraints is constructed and solved. The core objective of this model is to minimize the weighted sum of squared residuals between the estimated value of the natural electromagnetic pulse signal and the actual observed value, so as to restore the true geological origin characteristics of the signal as much as possible. Introducing an L1-norm sparsity constraint term into the optimization objective function helps to avoid overfitting, enhances the physical interpretability of the solution, and ensures that the finally extracted signal mainly focuses on the components highly responsive to geological stress changes. Through this process, a set of preliminarily optimized signal data is obtained. After obtaining the preliminarily optimized signal data, the singular value decomposition method is used to perform eigenanalysis on the system matrix involved in the integer linear programming model, and the condition number of the system matrix is calculated by decomposing its singular value spectrum. The condition number reflects the ill-conditioning degree of the system matrix. If its value is significantly greater than 1000, it indicates that there are problems such as numerical instability or overstretching of the solution space during the system solution process. Therefore, the Tikhonov regularization method is introduced to correct it. Tikhonov regularization effectively reduces the volatility of the solution and improves its robustness by adding a penalty term related to the signal norm to the optimization objective function. Aiming at the baseline drift problem commonly existing in the actual measurement process of electromagnetic signals, baseline correction is performed on the above-stabilized signal data. A piecewise polynomial function is used to model the signal baseline. The entire signal sequence is divided into several time periods, and a polynomial curve is fitted for each period to describe the background offset caused by non-geological factors. The polynomial coefficients are solved by minimizing the deviation between each segment of the signal and its fitted curve, so as to restore the signal baseline. After baseline correction, the signal waveform returns to its true undulation, eliminating the influence of low-frequency trend terms caused by non-structural factors such as environmental disturbances and equipment drift. An adaptive fuzzy correction algorithm is introduced to eliminate the ambiguity of the corrected signal. This algorithm adopts a Mamdani fuzzy inference model constructed based on expert knowledge and uses Gaussian membership functions to classify the input signal characteristics. The input variables of the fuzzy system include the time scale of the signal, spatial correlation, and environmental parameters during observation, while the output variable is the ambiguity correction coefficient for each segment of the signal. Through reasoning with the preset If-Then rules in the fuzzy rule base, the fuzzy control system adaptively generates the optimal correction strategy according to different input combinations, so as to classify, merge, and clarify the ambiguous information in the signal, and finally output the signal data with ambiguity eliminated. To solve the problem of electromagnetic signal positioning error caused by complex transmission paths or signal source ambiguity, time difference of arrival constraints are introduced based on the ambiguity-eliminated signal, an overdetermined observation equation set is constructed, and the iteratively reweighted least squares method is used to solve it. This method effectively deals with the non-linear characteristics in signal propagation by dynamically updating the weights of signal components in each round and gradually converging to the global optimal solution. During the iterative process, when the Euclidean distance change is less than the set threshold (such as Or reaching the maximum number of iterations (such as 100 times) as the stopping condition to ensure the best balance between the accuracy and the convergence speed of the calculation results. Obtain a set of electromagnetic signal samples with eliminated ambiguity and high spatio-temporal consistency. Conduct quality assessment on each segment of the signal. Quantitatively evaluate the integrity, effectiveness, and usability of the signal by calculating indicators such as its signal-to-noise ratio, amplitude stability, and phase coherence, and accordingly screen out the signal segments that meet the requirements of engineering applications. All the high-quality signal data passing the screening are uniformly archived to form a dataset of corrected electromagnetic pulse signals.
[0023] Step S104: Decompose the dataset of corrected electromagnetic pulse signals and construct a mapping relationship to obtain an estimated value of the geological stress distribution;
[0024] Specifically, high-quality electromagnetic pulse signal data is input into a time-reassigned multi-synchronous compression transform algorithm. This algorithm dynamically analyzes the local frequency characteristics of the signal and adaptively adjusts the window length to achieve fine decomposition of the signal. To ensure optimal time-frequency resolution in different frequency bands, the algorithm introduces an adaptive window length function, such that a longer window function is applied to signal segments with slow frequency changes to improve frequency resolution, while a shorter window function is used for segments with rapid frequency changes to enhance time resolution. After the transformation process, the corrected electromagnetic pulse signal is decomposed into several single-component functions with narrow frequency bands. Each function represents a specific frequency component and has a clearer physical feature directivity. From the narrowband single-component functions, a set of high-dimensional signal feature parameter sets related to geological stress activities is extracted. This parameter set includes the central frequency, bandwidth, distribution characteristics of signal energy in the time domain and frequency domain, duration, and the rate of change of instantaneous frequency of each single component, etc. These parameters comprehensively reflect the spectral structure, energy evolution, and dynamic propagation behavior of the electromagnetic pulse. Regarding the initial electromagnetic pulse wavefront information in the corrected electromagnetic pulse signal dataset, the signal responses at corresponding moments among multiple receiving points are selected as the basis. By using the differences in these response times to measure the time difference of arrival, a set of time difference of arrival data is formed, which reflects the path characteristics of the electromagnetic pulse propagating from the source point to each measuring point. Combining the known spatial distribution coordinate information of the receiving devices, the time difference of arrival data and the spatial distance matrix are jointly modeled to construct an overdetermined form of geological stress localization equation. Its basic expression form is V·T = D, where V represents the velocity tensor of electromagnetic wave propagation in the geological medium, T is the time difference of arrival matrix, and D is the spatial distance matrix between the measuring points. Since this equation belongs to a typical multi-variable overdetermined problem, the singular value decomposition method is used for robust solution to obtain a spatially anisotropic propagation velocity tensor, which reflects the propagation ability of electromagnetic waves in different directions and indirectly reveals the elastic properties and stress distribution characteristics of the underground medium. After completing the propagation modeling, to establish a quantitative mapping relationship between the signal characteristics and geological stress parameters, the extreme gradient boosting algorithm, i.e., the XGBoost model, is introduced to train and analyze the signal feature parameter set. As an enhanced tree model, XGBoost can efficiently capture non-linear mapping relationships while taking into account overfitting control and feature selection capabilities. In this step, using the known geological stress data in the historical monitoring data as labels, a supervised learning model is constructed. The input features are the previously extracted spectrum-based, energy-based, and time-frequency dynamic-based parameters, and the model output is the stress magnitude, principal stress direction, and stress gradient value of the target area. During model training, hyperparameters such as the number of decision trees, depth, and learning rate are initially set, and the complex relationships between the data are gradually fitted through the standard gradient boosting mechanism.To improve the generalization ability and prediction accuracy of the model, the importance of different input features in the model is sorted and screened. During the screening process, a feature importance analysis method is adopted to identify the top M key features with the greatest contribution in the stress prediction task. This process effectively eliminates redundant features, reduces the model complexity, and improves the training efficiency. The hyperparameters of the XGBoost model are automatically adjusted in combination with the Bayesian optimization method, including the regularization term coefficient, the maximum depth of the tree, the learning rate, the minimum sample weight, etc., to maximize the performance on the validation set. The optimization process searches for the optimal solution point in the high-dimensional parameter space through the Bayesian framework, avoiding the inefficiency of the traditional grid search. The geological stress prediction model constructed through the above steps can accurately map the electromagnetic signals in any monitoring area and output the estimated values of the structured geological stress distribution, including the direction vector of the principal stress, the absolute value of the stress intensity, and the stress gradient distribution between different regions.
[0025] Fast Fourier transform is performed on the corrected electromagnetic pulse signal data set. The fast Fourier transform method is used to process the signal in the entire frequency band. With its advantages in computational efficiency and frequency resolution, a series of initial spectral characteristic parameters such as energy distribution, main frequency peak, and frequency band extension of the electromagnetic signal in the frequency domain are quickly obtained. The local eigenvalues of the Hessian matrix are calculated for each frequency band based on the initial spectral characteristic parameters of the signal. The frequency change curvature model of the signal is constructed using the local Hessian matrix, and the main eigenvalues of the Hessian matrix corresponding to each frequency point are calculated to obtain the intensity of the signal change in the frequency band where the point is located. The more drastic the frequency change interval, the larger the eigenvalue of the Hessian matrix, reflecting that the signal in this frequency band has a strong time-varying property, and a higher time resolution should be given in subsequent analysis. Based on this analysis result, the window length scaling factor is constructed in combination with the frequency change rate information. By designing an adaptive window length function, the time window length of different frequency intervals of the signal is dynamically allocated, so that the low-frequency stable segment obtains a longer window function to enhance the frequency resolution, while the high-frequency mutation segment is allocated a shorter window function to improve the time resolution, forming an adaptive window length function system that automatically adjusts with frequency changes. The short-time Fourier transform is used as the basic framework of time-frequency analysis, and the adaptive window length function is embedded in the short-time Fourier transform process. For each frequency band of the signal, the corresponding window function length is used for piecewise convolution calculation to generate a time-frequency distribution map with high-resolution time-frequency structure. The time-frequency distribution map is subjected to time redistribution operation. According to the instantaneous frequency and group velocity information of the signal, the energy distributed in the fuzzy area is refocused to its real time-frequency position, so as to improve the energy concentration and structural clarity of the time-frequency map, and form an enhanced time-frequency distribution feature map with stronger readability. The enhanced time-frequency distribution feature is input into the synchronous compression kernel function to perform ridge extraction, and the synchronous depth control strategy is applied to stabilize the ridge extraction process. The synchronous compression kernel function is used to identify the energy peak coherent trajectory and extract the time-frequency ridge sequence representing the main path of the signal. These ridges correspond to the most stable and geophysical frequency components of the signal in the time-frequency map. The synchronous depth control adjusts the threshold parameters in the extraction process according to the local consistency index of the signal at multiple scales to ensure that each extracted ridge has good energy coherence and physical consistency. Through this process, multiple original single-component functions are obtained, each of which corresponds to a signal sub-component with clear physical meaning and stable frequency characteristics, reflecting the independent electromagnetic response path generated during a specific geological activity. Energy contribution evaluation is performed on the original single-component function. This evaluation calculates the proportion of each single component in the total signal energy, and selects the parts with significant energy proportion and stable structure as the final analysis object. Only signal components with energy contribution exceeding the preset threshold are retained, and low-energy, dispersed structure, and difficult-to-interpret secondary components are eliminated, and finally multiple narrow-band single-component functions with high energy integration, spectral stability and time domain resolution are obtained.
[0026] Step S105: Construct a three-dimensional model of the geological stress concentration area based on the estimated value of the geological stress distribution.
[0027] Specifically, taking the estimated value of the geological stress distribution as the modeling basis, and using the Kriging interpolation method in spatial statistics to reconstruct the continuous field of the stress parameters of all discrete sampling points. As the optimal unbiased estimation method in geological spatial interpolation, Kriging interpolation takes into account the spatial distance between measurement points and introduces the semi-variance function to describe the spatial correlation between measurement points, showing high fitting accuracy when dealing with physical quantities with spatial heterogeneity such as geological stress. Through this method, the stress estimated values distributed at different depths and positions are uniformly interpolated into a structurally coherent and continuously distributed stress field data, forming a spatial stress tensor field expressed in the form of a three-dimensional grid, where each grid cell contains the magnitude and direction information of the principal stress. Based on the obtained stress field data, tensor analysis is carried out to extract the principal stress direction and magnitude in each spatial unit, and calculate the three principal stress values , and and their corresponding direction cosines to construct a complete expression of the principal stress ellipsoid. By statistically analyzing the distribution characteristics of the principal stress in the whole region, the regions with high stress gradient and large principal stress deviation are identified, and such regions are candidates for the spatial distribution of potential stress concentration areas. In order to achieve intelligent discrimination, combined with the regional average stress state, a set of clear discrimination criteria for stress concentration areas are constructed. For example, it is stipulated that if the maximum principal stress in a certain region exceeds 1.5 times the average principal stress of the whole region, and the difference between the maximum and minimum principal stresses exceeds 2 times its regional average difference, then this region is marked as a preliminarily identified stress concentration area. The Markov random field model is introduced to perform spatial optimization on the preliminary stress concentration area. By setting the local spatial adjacency relationship and state transition probability, this model establishes a statistical correlation between the regional stress data and the neighborhood relationship, making the optimization result take into account both data-driven and spatial structure rules. The model regards each spatial grid cell as a state node and determines its final attribution by comparing the stress states of its adjacent nodes, thereby eliminating the isolated pseudo-abnormal points caused by abnormal sampling or interpolation errors in the identification results, and at the same time enhancing the geometric coherence and geological structure consistency of the concentration area in space. After optimization, a more reasonable, clearly bounded and structurally continuous stress concentration area in space is obtained. Based on the optimized stress concentration area, a Stress Concentration Index (SCI) evaluation system is established to achieve quantitative grading of the stress concentration degree in different regions. The SCI index comprehensively considers three factors: the relative strength of the maximum principal stress, the difference between the principal stresses, and the distance from the fault. Its expression form is defined as:
[0028] , where and respectively represent the average values of the maximum principal stress and the stress difference within the region, and D is the normalized distance index between the center point of the unit and the nearest active fault. The SCI value of each grid cell is calculated by this formula, and then the stress concentration area is divided into three levels, namely level I (SCI > 4.0), level II (2.0 ≤ SCI ≤ 4.0), and level III (1.0 ≤ SCI < 2.0) according to the set threshold, realizing the quantitative expression and classification processing of stress risk. Based on the classified stress concentration area data, 3D modeling work is carried out to integrate the spatial distribution of geological stress and geological background information, and a 3D model of the comprehensive geological stress concentration area containing various information elements is constructed. The model takes the stress tensor field as the backbone and is supplemented by lithology distribution, fracture development degree, and groundwater distribution data as geological structure constraints. To improve the accuracy and adaptability of the model, an adaptive grid meshing technology is adopted to automatically refine the grid in the area with a large stress gradient to more accurately describe the complex geological response characteristics. The model appears as a multi-layer nested structure in 3D space and is visualized through various forms such as isosurface maps, slice maps, and vector field maps, presenting the internal structure of the stress concentration area, the stress evolution trend, and its coupling relationship with the surrounding tectonic system in a complete manner.
[0029] In the embodiment of the present application, the spatial resolution is improved by a multi-vector electromagnetic induction coil monitoring network with a grid-like distribution. An equilateral triangle arrangement pattern and an orthogonal three-component receiving method are adopted to achieve an all-round three-dimensional monitoring of natural electromagnetic pulse signals, significantly enhancing the stress field boundary recognition accuracy and reducing the blind area. The systematic signal preprocessing process includes power frequency filtering, wavelet denoising, and adaptive band-pass filtering. Combining with the time-frequency analysis of the Hilbert-Huang transform, the effective components in the weak electromagnetic pulse signals are effectively extracted, and environmental interference and system noise are suppressed. A method combining an integer linear programming model and an adaptive fuzzy correction algorithm is introduced to solve the problems of signal baseline drift and electromagnetic signal ambiguity that have long troubled this field. The multi-synchronous compression transform with time reassignment breaks through the limitations of traditional Fourier analysis and can effectively process non-stationary signals. Combining with the signal feature and geological stress parameter mapping relationship constructed by the extreme gradient boosting algorithm, the high-precision stress field reconstruction under complex geological conditions is realized. The stress concentration index evaluation system correlates the stress field characteristics with geological structures and optimizes the spatial constraints through the Markov random field model, significantly reducing the false anomaly recognition rate. The 3D model of the geological stress concentration area adopts an adaptive grid meshing technology to finely depict the spatial geometric characteristics and internal stress distribution of the stress concentration area.
[0030] In a specific embodiment, the process of executing step S101 may specifically include the following steps:
[0031] Perform a geological structure complexity analysis on the monitoring area to obtain the geological structure complexity analysis results, and determine an equilateral triangle arrangement pattern based on the geological structure complexity analysis results to obtain a receiving device layout plan;
[0032] Install three orthogonal vector electromagnetic induction coils on each receiving device in the receiving device layout plan to obtain an electromagnetic signal receiving unit;
[0033] Integrate a GPS clock module for each electromagnetic signal receiving unit to obtain a time synchronization acquisition system, and connect the electromagnetic signal receiving units according to the time synchronization acquisition system and a low-power wireless transmission module to form a mesh network to obtain a distributed monitoring network;
[0034] Select a reference site in the monitoring area and record the external electromagnetic interference signals in the area based on the reference site to obtain reference interference data;
[0035] Monitor the electromagnetic pulse signals of the distributed monitoring network, and combine the reference interference data to mark the valid signals to obtain the original electromagnetic pulse signals containing the X, Y, and Z direction components.
[0036] Specifically, conduct a geological structure complexity analysis of the monitoring area to systematically identify the tectonic units, fault systems, lithological combinations and their distribution laws within the area. This analysis process comprehensively utilizes geological maps, remote sensing images, drilling results and historical earthquake data. By constructing a geological structure complexity evaluation model, quantitatively evaluate key elements such as joint density, fault intersection degree, number of tectonic intersection points, number of stratigraphic unconformity interfaces, etc. within the area, and conduct multi-factor weighted calculations in combination with the frequency of seismic activities, crustal movement speed and changes in geomorphic units to obtain the geological structure complexity analysis result expressed in the form of a structure complexity level. This analysis result is spatially divided into multiple regional units, each unit is assigned a structure complexity level label, and is output in the form of a layer. According to the geological structure complexity analysis result, formulate the layout strategy of the receiving device, adopt an equilateral triangle grid as the basic arrangement mode to ensure an isotropic signal acquisition ability in space. In areas with high structure complexity, such as areas with dense fault intersections, fold-intensive areas or lithological mutation areas, set a higher density of receiving nodes, and control the side length of the triangle between 50 meters and 100 meters to improve the local spatial resolution; while in areas with relatively stable or simple structures, appropriately reduce the layout density and expand the side length to 150 meters to 200 meters to reduce system resource consumption and simplify the layout engineering quantity. Through the layout strategy driven by structure complexity, form a layout plan of the receiving device that not only has macroscopic coverage ability but also takes into account local high-resolution acquisition. At each layout point, install an electromagnetic signal receiving unit according to a unified standard. In each receiving device, three orthogonally arranged vector electromagnetic induction coils are fixed, corresponding to the three spatial directions of X, Y, and Z respectively, for synchronously collecting the electromagnetic field change data in each direction. The inductance coils are all connected to high-sensitivity fluxgate sensors, whose range is set to ±100 nT and the sensitivity reaches 0.1 pT / √Hz to ensure that extremely low-amplitude electromagnetic signals generated during the release of natural geological stress can be captured. Introduce magnetic shielding and non-magnetic housing material design in the inductance acquisition structure to prevent the device structure itself from interfering with or shielding the signal. At the same time, environmental sensors such as temperature, humidity, and vibration are equipped to record the acquisition environment parameters in real time, which are important reference information for subsequent signal processing and interference identification. To achieve time synchronization between receiving units, a high-precision GPS clock module is integrated in each electromagnetic signal receiving unit. This module supports millisecond-level synchronization accuracy. By receiving the standard time signal of the global navigation satellite system, it provides a unified time reference for electromagnetic signal acquisition, ensuring that all devices have timestamp consistency when collecting the same natural electromagnetic event. After synchronization, the receiving units compare the time of arrival of electromagnetic signals from the same source point and invert the propagation path, supporting subsequent signal positioning and stress field modeling work. On the basis of achieving time synchronization, add a low-power wireless transmission module to each receiving device, and construct a mesh network architecture through a multi-hop communication mechanism to realize automatic interconnection and data collaborative forwarding between nodes.The mesh network supports dynamic path reconstruction and self-repair mechanisms between devices. When a node fails or communication is interrupted, other nodes automatically rebuild the transmission path to ensure the continuous and stable operation of the entire network. Through this distributed mesh network, all receiving units upload the collected data to the central data processing server in real time to realize data centralized processing, remote monitoring and fault diagnosis functions, and build a data collection system with full coverage, high reliability and strong coordination. After the distributed collection system is deployed, several reference sites are selected at the outer edge of the monitoring area. These reference points try to avoid the main geological tectonic activity belt and strong artificial interference sources. Their main task is to record non-geological electromagnetic interference signals from outside the area or background sources. The data collected by the reference site is used to build a background interference library, including typical interference modes such as lightning activity, electromagnetic leakage of power lines, and radio frequency radiation. Its time, frequency, and waveform characteristics will be used as comparison templates to assist in the subsequent identification of effective geological electromagnetic events. After the entire distributed monitoring network is put into operation, each receiving unit will conduct uninterrupted real-time monitoring of the electromagnetic pulse signals that continue to appear in the natural environment. All collected signals are calibrated with GPS timestamps and then uploaded to the central server. The system automatically compares and analyzes them with the background interference data of the reference site. Multi-dimensional feature extraction and comparison methods such as frequency domain matching, waveform similarity, polarization parameters and wave arrival time are used to accurately identify and mark valid signals with clear sources and geological characteristics. The selected signals are the original electromagnetic pulse signals containing the three directional components of X, Y, and Z.
[0037] In a specific embodiment, the process of executing step S102 may specifically include the following steps:
[0038] The original electromagnetic pulse signal is subjected to interference filtering to obtain a preliminary filtered signal, and the preliminary filtered signal is subjected to wavelet threshold denoising to obtain a denoised signal;
[0039] Perform time domain segmentation processing on the denoised signal to obtain a segmented time domain signal, and perform time-frequency analysis on the segmented time domain signal to obtain a signal feature matrix;
[0040] The polarization parameters including ellipticity, inclination and azimuth are calculated based on the three-component vector signal in the signal characteristic matrix, and adaptive notch filtering is performed in combination with the external interference signal recorded at the reference station to obtain a pure geological source signal.
[0041] The pure geological source signals are subjected to spatial correlation analysis, the mutual correlation coefficients and phase differences between adjacent measuring points are calculated, and they are organized in a unified format to form a multidimensional data cube containing time, space, amplitude, frequency, and phase.
[0042] Specifically, taking the three-component raw electromagnetic pulse signals collected at each measuring point in the monitoring area as input data, which are strongly affected by power system interference, lightning electromagnetic waves, communication emission signals, and man-made electromagnetic sources, and performing interference filtering on them. A notch filter is designed to specifically suppress power frequency interference such as 50 Hz and its harmonics. By using precise band-stop filter parameter settings, the noise brought by the power system is eliminated to the greatest extent without affecting the effective components in the target frequency band. At the same time, for interference sources with a relatively wide frequency band and a fast-changing time scale, a multi-channel filter bank is introduced, and band-pass filtering is performed according to the frequency domain analysis results to obtain a preliminary filtered signal. For the preliminary filtered signal, a wavelet threshold denoising method is introduced for multi-scale decomposition. The db4 in the Daubechies wavelet family is selected as the mother wavelet. By performing 5-layer discrete wavelet decomposition on the signal, it is divided into approximation coefficients and detail coefficients at multiple scales, and then the high-frequency noise in the detail coefficients is suppressed according to the soft threshold principle. The soft threshold function has good smoothing performance, can effectively eliminate high-frequency noise interference while retaining the mutation characteristics of the effective signal, and is suitable for the signal structure of natural electromagnetic signals that are non-stationary, non-linear, and contain mutation characteristics. Through the process of denoising layer by layer and then reconstruction, a denoised electromagnetic signal is obtained. This signal shows continuous waveform, clear mutation points, and stable baseline in the time domain, and has a concentrated and clearly structured main frequency distribution in the frequency domain. The denoised signal is processed by time-domain segmentation to achieve time-varying extraction of non-stationary signal characteristics. The sliding window mechanism is used to divide the entire signal sequence into multiple time periods with a length of 60 seconds, and adjacent windows are set with a 50% overlap rate to improve time resolution and enhance the ability to capture event continuity. Each time-domain segmented signal is input into the time-frequency analysis module, and the Hilbert-Huang transform is selected as the main analysis tool. This method decomposes each segment of the signal into a set of intrinsic mode functions through empirical mode decomposition, and then performs the Hilbert transform on each mode function to extract time-varying characteristic parameters such as instantaneous frequency, instantaneous amplitude, and instantaneous phase. All these information are summarized into a multi-dimensional signal feature matrix, where each row corresponds to a time period, and each column represents a characteristic index, such as average frequency, frequency change rate, amplitude envelope, number of modes, etc., forming a characteristic description result with time series and frequency structure. Polarization parameter analysis is performed on the three-component vector signals contained in the signal feature matrix to identify the polarization characteristics of electromagnetic waves during spatial propagation. By combining the signals in the X, Y, and Z directions to construct a spatial vector trajectory, the propagation ellipse is restored in three-dimensional space, and the ellipticity of the ellipse is calculated to judge its deviation from circular polarization. At the same time, the signal inclination angle is derived according to the main axis direction to measure the inclination degree of the propagation path relative to the ground surface, while the azimuth angle reflects the main propagation direction of the horizontal component of the signal. The calculation of polarization parameters helps to distinguish the source type and propagation mode of the signal. Especially in the multi-source or multi-path propagation scenario, the polarization characteristics are used to exclude secondary interference sources.Implement an adaptive notch filtering mechanism by combining the background interference data recorded at the reference site. The reference site is located in an area far from the main tectonic zones and engineering activities, and the signals collected there are regarded as a pure background interference model. By constructing the interference template spectrum and performing frequency-domain comparison and analysis with the target signal, the main frequency and its energy distribution of the interference components are extracted. Then, notch filters are designed specifically in the target signal to dynamically attenuate these components, so as to eliminate the influence of external disturbances while retaining the effective geological components, and output a pure geological source electromagnetic signal. Perform spatial correlation analysis on the pure geological source signal to evaluate the consistency and spatial coupling degree of the signals among various measurement points. By calculating the cross-correlation coefficient between adjacent receiving points, analyze the similarity of their signal waveforms in time; at the same time, extract the instantaneous phase of each signal segment and calculate the phase difference to reveal the propagation delay characteristics of the signal wavefront. Organize all the above processing results according to the unified data structure standard to construct a multi-dimensional data cube with temporal sequence, spatial distribution, and complete frequency characteristics. Each data unit of this cube contains parameters such as timestamp, spatial coordinates of the measurement point, normalized amplitudes in three directions, main frequency, and phase.
[0043] In a specific embodiment, the process of executing step S103 may specifically include the following steps:
[0044] Input the multi-dimensional data cube into an integer linear programming model, and by minimizing the weighted residual sum of squares between the signal estimated value and the observed value and introducing the L1 norm sparse constraint, obtain the preliminarily optimized signal data;
[0045] Perform singular value decomposition on the preliminarily optimized signal data to obtain the condition number, and apply the Tikhonov regularization method according to the condition number to obtain the stabilized signal data;
[0046] Perform baseline correction on the stabilized signal data using a piecewise polynomial function to obtain the signal after baseline correction;
[0047] Apply Mamdani fuzzy inference to the signal after baseline correction using an adaptive fuzzy correction algorithm, and classify the signal characteristics according to the Gaussian membership function to obtain the signal data with ambiguity eliminated;
[0048] Perform iterative reweighted least squares processing on the signal data with ambiguity eliminated using the time difference of arrival constraint condition to obtain the signal with ambiguity eliminated;
[0049] Calculate the signal-to-noise ratio, amplitude stability, and phase coherence parameters according to the signal with ambiguity eliminated, and screen the signal segments that meet the requirements according to the quality index to obtain the corrected electromagnetic pulse signal dataset.
[0050] Specifically, the multi-dimensional data cube is input into an integer linear programming model, in which an objective function is defined to minimize the weighted sum of squared residuals between the estimated value and the observed value of the electromagnetic signal as the optimization objective, and by introducing the L1-norm sparsity constraint term, the solution space becomes more compressible and physically interpretable. This optimization problem is expressed as minimizing , where A represents the system matrix, x is the signal vector to be solved, b is the actual observed vector, λ is the sparse constraint coefficient, and here the weights are adaptively allocated according to the noise level and signal credibility of the signal, thereby enhancing the fitting accuracy of the model in the high signal-to-noise ratio region and suppressing the interference contribution in the noise region. After solving this integer linear programming, a set of preliminarily optimized signal data is obtained. Perform singular value decomposition on the preliminarily optimized signal data, and decompose the system matrix A into UΣV TForm, and calculate the condition number of the matrix by observing the change of the singular value sequence. If the condition number is high, it indicates that there is ill-conditioning in the system matrix during the solution process, which is likely to cause instability of the solution or high sensitivity to initial perturbations. At this time, the Tikhonov regularization method is introduced for robust processing. Tikhonov regularization controls the oscillation behavior of the solution by adding an L2 norm penalty term of the solution vector to the objective function, making the solution result more robust in the face of redundant information, thus obtaining a set of stabilized signal data after regularization processing. Baseline correction is performed on the stabilized signal data to eliminate the low-frequency trend or baseline drift problems caused by factors such as environmental changes and instrument zero drift. This step adopts a piecewise polynomial function modeling strategy, divides the entire signal into several time periods, and fits a polynomial baseline curve for each segment separately. The common choices are quadratic or cubic polynomials to accurately capture the non-linear change characteristics of the baseline, and the polynomial coefficients are determined by the least squares method. Then, with this polynomial baseline curve as a reference, the fitting result is subtracted from the original signal point by point to obtain the signal after baseline correction, making the signal fluctuate around zero numerically and not contain additional drift components, significantly improving the accuracy of subsequent ambiguity elimination and physical parameter extraction. On the basis of baseline correction, to enhance the clarity of the signal structure and reduce interpretation ambiguity, an adaptive fuzzy correction algorithm is introduced to process the ambiguity of the signal, and the Mamdani fuzzy inference system is used as the core inference mechanism. This system represents the fuzzy sets of the input signal characteristic variables with Gaussian membership functions, sets the input dimensions including the instantaneous frequency change rate, amplitude fluctuation degree, phase continuity index, etc. of the signal, and uses the baseline correction factor as the output quantity. Fuzzy inference is carried out through the If-Then form rules in the rule base to output an adaptively adjusted ambiguity correction value, obtaining a set of signal data without ambiguity. The iterative reweighted least squares method is used to process the signal data without ambiguity using the time difference of arrival constraint condition. Calculate the time difference of arrival of the signal based on the GPS timestamp and the spatial measurement point position, construct a propagation time error function, use the signal estimated value as the optimization variable, minimize the propagation path error and the signal residual, and punish the outliers by assigning different weights to the results of each round of calculation, thereby suppressing the error amplification caused by inconsistent propagation models. During the iteration process, the weight matrix and the residual function are continuously updated until the Euclidean distance between the solutions of two consecutive iterations is less than the preset threshold or the maximum iteration number limit is reached, and finally a set of signals without ambiguity with optimal time consistency and spatial propagation characteristics is output. To ensure the effectiveness of the final result in data analysis and model training, multi-dimensional quality evaluation is performed on the signals without ambiguity, and quantitative scoring is carried out by calculating indicators such as the signal-to-noise ratio, amplitude stability, and phase coherence of each segment of the signal. The signal-to-noise ratio reflects the proportion of the effective components in the signal, the amplitude stability measures the fluctuation degree of the signal in a continuous time window, and the phase coherence is used to judge whether the signal shows the physical consistency of continuous propagation.By setting evaluation thresholds for these metrics, high-quality signal segments are screened out, abnormal data that does not meet the standards is eliminated, and a set of corrected electromagnetic pulse signal datasets with standardization, structurization, and strong physical interpretability is generated.
[0051] In a specific embodiment, the process of executing step S104 may specifically include the following steps:
[0052] Apply the time-reassigned multi-synchronous compression transform algorithm to set the adaptive window length function, and dynamically adjust the local frequency characteristics of the corrected electromagnetic pulse signal dataset to obtain multiple narrowband single-component functions;
[0053] Extract a set of signal characteristic parameters including center frequency, bandwidth, energy distribution, duration, and instantaneous frequency change rate from the multiple narrowband single-component functions;
[0054] For the initial electromagnetic pulse wavefront information in the corrected electromagnetic pulse signal dataset, measure the time difference of arrival between different measurement points to obtain the time difference of arrival data;
[0055] Based on the time difference of arrival data and combined with the spatial distribution coordinates of the receiving device, construct an overdetermined geostress localization equation, and solve it to obtain the electromagnetic wave propagation velocity tensor;
[0056] Adopt the extreme gradient boosting algorithm to establish the mapping relationship between the set of signal characteristic parameters and the geostress parameters to obtain the geostress prediction model;
[0057] Through feature importance analysis, screen out the top M features that contribute the most to stress prediction, and combine the Bayesian optimization method to automatically adjust the hyperparameters of the extreme gradient boosting algorithm to optimize the geostress prediction model, and output the estimated values of geostress distribution including the principal stress direction, stress magnitude, and stress gradient.
[0058] Specifically, the corrected electromagnetic pulse signal data set is used as input, and the time-reassigned multi-synchronous compression transform algorithm is used to perform high-precision decomposition of its non-stationary time-frequency characteristics. The algorithm introduces a window length adjustment mechanism based on the calculation of the local frequency change rate of the signal, that is, by analyzing the frequency intensity of the signal in different time windows, the length of the analysis window function is dynamically set, so that a longer window is used in the frequency stable area to enhance the frequency resolution, and a shorter window is used in the area of drastic frequency changes to enhance the time positioning ability, forming an adaptive window length function system, so that the signal shows a highly concentrated time-frequency distribution structure after transformation. This process effectively compresses the background noise and improves the accuracy of multi-component recognition while retaining the main physical characteristics of the signal. The signal processed by the multi-synchronous compression transform algorithm is decomposed into multiple narrow-band single-component functions with clear frequency structure and concentrated energy distribution. Each function corresponds to a frequency principal component with independent propagation path and geological significance. A set of characteristic parameters with strong representativeness and clear physical meaning are extracted from the above multiple narrow-band single-component functions to construct a signal characteristic parameter set. The parameter set includes the center frequency, which reflects the main frequency position where the signal energy is concentrated; the bandwidth, which indicates the frequency diffusion range and the stability of the propagation source; the energy distribution, which describes the change law of the signal strength in the frequency domain per unit time; the duration, which reflects the length and stability of the single component signal in the time dimension; and the instantaneous frequency change rate, which reveals the dynamic disturbance characteristics of the frequency component during the propagation process. The initial wavefront information in the electromagnetic signal is extracted, especially focusing on the first wave response time of each pulse event at different receiving points. By measuring the first wave arrival time difference between multiple measuring points, a wave arrival time difference data set is formed. This data contains the electromagnetic signal transmission delay between each pair of measuring points, which directly reflects the signal propagation speed and path differences. Combining these time difference data with the known three-dimensional spatial distribution coordinates of the receiving device, an electromagnetic signal propagation path model is established, and then an overdetermined geological stress location equation is constructed. Its basic form is V·T=D, where V is the unknown electromagnetic wave propagation velocity tensor, T is the wave arrival time difference matrix, and D is the spatial distance matrix between each measuring point. Since the number of sampling points is much larger than the tensor dimension, the equation is a typical overdetermined system. Robust solution algorithms such as singular value decomposition are used to infer the tensor characteristics of the propagation speed of electromagnetic waves in different directions while ensuring the stability of the solution. The propagation speed tensor represents the medium response capability of the electromagnetic signal, and also indirectly reflects the electrical anisotropy and stress state distribution of the underground rock mass. After obtaining the propagation tensor and signal characteristic parameter set, the extreme gradient boosting algorithm (XGBoost) is used to establish a nonlinear mapping relationship between the two and the actual geological stress parameters. In this supervised learning framework, the input of the model is the multi-dimensional variables in the signal characteristic parameter set, and the output is the geological stress index of the target area, including the principal stress magnitude, stress direction vector, and stress gradient value.As an ensemble learning method, XGBoost achieves strong prediction performance by constructing multiple weak classification tree models and using weighted accumulation, effectively identifying the non-linear dependence structure between variables and outputs in a high-dimensional complex feature space. During the model training process, the data with known stress labels in the historical monitoring data is used as the training set, and the generalization ability of the model is evaluated by cross-validation. Through this training mechanism, the XGBoost model can extract the potential relationship between multi-variable combinations and stress responses, and construct a geological stress prediction model with high robustness and high prediction accuracy. After the model construction is initially completed, a feature importance analysis method is introduced to rank the importance of the input signal feature parameter set to identify the top M features that have the most significant impact on the prediction results. Feature importance is statistically evaluated based on indicators such as split gain, frequency, or coverage inside the model, which can effectively eliminate redundant or noisy features, thus simplifying the model structure and enhancing interpretability. The Bayesian optimization method is combined to automatically adjust the core hyperparameters of the XGBoost model, including the maximum depth of the tree, learning rate, subsampling rate, regularization coefficient, etc. Bayesian optimization probabilistically models the parameter space by constructing a surrogate model and selects the optimal parameter combination according to the expected improvement function strategy, achieving global optimality of performance while controlling the training time cost. The geological stress prediction model after multiple rounds of training and optimization can quickly map any newly input electromagnetic signal feature parameter set into a stress distribution estimate value including the principal stress direction, stress magnitude, and stress gradient.
[0059] In a specific embodiment, the process of performing the multi-synchronous compression transform algorithm with time reassignment to set the adaptive window length function and dynamically adjust the local frequency characteristics of the corrected electromagnetic pulse signal dataset to obtain multiple narrowband single-component functions may specifically include the following steps:
[0060] Perform a fast Fourier transform on the corrected electromagnetic pulse signal dataset to obtain the initial spectral characteristic parameters of the signal;
[0061] Calculate the local eigenvalues of the Hessian matrix for each frequency band based on the initial spectral characteristic parameters of the signal, and apply the multi-synchronous compression transform algorithm with time reassignment to formulate a window length scaling factor in combination with the frequency change rate to obtain the adaptive window length function;
[0062] Apply the short-time Fourier transform to the corrected electromagnetic pulse signal dataset, and use the adaptive window length function to set the time window lengths of different frequency bands to obtain a time-frequency distribution map;
[0063] Perform a time reassignment operation on the time-frequency distribution map to obtain enhanced time-frequency distribution features;
[0064] Input the enhanced time-frequency distribution features into the synchrosqueezing kernel function to perform ridge extraction and synchrosqueezing depth control, obtaining multiple original single-component functions;
[0065] Evaluate the energy contribution degrees of the multiple original single-component functions to obtain multiple narrowband single-component functions.
[0066] Specifically, a fast Fourier transform is performed on the corrected electromagnetic pulse signal data set, and a full spectrum conversion operation is performed on each signal segment. With its efficient computing performance and global frequency perspective, the time domain signal is mapped to the frequency domain, and the energy distribution characteristics of the signal at each frequency are extracted, thereby obtaining the initial spectrum characteristic parameters containing key information such as the main frequency component, spectrum energy density, and spectrum width. The local eigenvalues of the Hessian matrix are calculated for each frequency band based on the initial spectrum characteristic parameters of the signal. A Hessian matrix based on the second-order derivative is constructed for the local area of each frequency band in the spectrum, and its eigenvalue is calculated to characterize the second-order curvature of the signal energy density change in the frequency band. The larger the eigenvalue, the more drastic the energy change and the more non-stationary the frequency band is, and a shorter time window function needs to be used to improve the time positioning capability; on the contrary, if the eigenvalue is small, it means that the energy change of the frequency band is slow and the frequency is stable, and a longer window function is used to improve the frequency resolution. According to the statistical results of the Hessian matrix eigenvalues and the frequency change rate index of the initial spectrum, the window length scaling factor is constructed, and the analysis window length parameters exclusive to each frequency band are defined accordingly to form a frequency-dependent and adaptively changing window length function. Based on the constructed adaptive window length function, the short-time Fourier transform is performed on the corrected electromagnetic pulse signal data set. The window length function is bound to the frequency region, and the corresponding window function length is selected for different frequency bands, so as to dynamically adjust the balance weight between time and frequency in the analysis. In the signal frequency mutation section, the window length is shortened to improve the time resolution; while in the frequency stable region, the window length is appropriately lengthened to improve the frequency resolution. This strategy improves the ability to analyze complex signal structures and generates a nonlinear time-frequency distribution diagram with a time-varying window structure. The time-frequency distribution diagram is time-redistributed. The fuzzy energy distributed in the non-central area is concentrated and compressed back to its real physical position to obtain a more structured time-frequency representation. By utilizing the instantaneous frequency and group delay information of the signal, each point of energy is remapped to its real time-frequency center, removing the frequency leakage effect and fuzzy diffusion caused by the window function, and obtaining an enhanced time-frequency distribution feature diagram. The enhanced time-frequency distribution features are input into the synchronous compression kernel function to perform ridge extraction and synchronous depth control operations. The essence of ridge extraction is to identify the main energy path in the time-frequency diagram, which corresponds to the main response trajectory generated by different physical components in the electromagnetic signal during the propagation process; and the synchronous compression kernel function combines local consistency evaluation with global curvature fitting to ensure that the extracted ridges have significant characteristics in energy intensity, frequency coherence and time consistency. At the same time, in order to prevent the misidentification of non-principal components during the extraction process, a synchronous depth control strategy is introduced. According to the embedding depth of the signal ridge in the time-frequency diagram, the ridge density and the relative weight relationship with the adjacent ridges, only the main ridge structure with stable structure and physical rationality is retained, and multiple original single-component functions are output. Each function is a structural unit extracted from the enhanced time-frequency diagram, with a clear time-frequency distribution range and stable physical propagation directionality.In order to compress the data scale and improve the analysis focus, it is necessary to evaluate the energy contribution of these original single-component functions. This evaluation process calculates the proportion of each single component in the total energy, and comprehensively scores its energy concentration, time-frequency stability, and propagation consistency. Components with high energy contribution, clear signals, and complete structures are selected and retained as the final results in the evaluation results, while redundant, low-intensity, or morphologically unstable components are removed. The multiple narrowband single-component functions finally obtained have high physical interpretability, good mathematical stability, and information decoupling ability.
[0067] In a specific embodiment, the process of executing step S105 may specifically include the following steps:
[0068] Based on the estimated value of the geological stress distribution, the stress parameters of discrete measurement points are interpolated into stress field data by using Kriging interpolation method;
[0069] Perform tensor analysis on the stress field data to obtain the principal stress distribution characteristics, and formulate a discrimination criterion for stress concentration areas according to the principal stress distribution characteristics to obtain the preliminarily identified stress concentration areas;
[0070] Apply the Markov random field model to the preliminarily identified stress concentration areas for spatial constraint optimization to obtain the optimized stress concentration areas;
[0071] Based on the optimized stress concentration areas, establish a stress concentration index evaluation system, and calculate the classified stress concentration areas according to the stress concentration index evaluation system;
[0072] Based on the classified stress concentration, conduct 3D modeling to generate a 3D model of the geological stress concentration area that includes stress tensor field, lithology distribution, fracture development degree, and groundwater distribution information.
[0073] Specifically, based on the estimated value of the geological stress distribution, the stress parameters of discrete measurement points are interpolated into stress field data by using Kriging interpolation method. Kriging interpolation fits the variogram between measurement points, fully considers the spatial autocorrelation and the weight relationship of the measurement point spacing, conducts interpolation calculation on the entire area, realizes the optimal estimation of stress values at unobserved points, keeps the estimation error minimized and the overall spatial structure reasonable, and generates a 3D continuous stress field data with a geological interpretation basis and mathematical stability. Based on this stress field data, tensor analysis is carried out to extract the principal values of the stress tensor and their principal axis directions in each grid cell, specifically including the maximum principal stress 、the intermediate principal stress 、the minimum principal stress And their corresponding direction cosines are used to construct the expression form of the stress ellipsoid at each point. By performing cluster analysis and gradient operation on the spatial distribution of the principal stress values in the regional stress tensor field, areas with highly concentrated and rapidly changing local stresses, that is, stress abnormal aggregation zones, are identified. To standardize the identification mechanism, clear criteria for judging stress concentration areas are formulated. For example, it is stipulated that the maximum principal stress in a certain area exceeds 1.5 times the average principal stress of the whole area, and the difference from exceeds more than 2 times the average stress difference, then this area is marked as a potentially stress-concentrated area initially identified. The Markov random field model is applied to optimize the spatial constraints for the initially identified stress concentration areas. By defining the adjacency relationship and state transition probability between grid cells, the model establishes the joint probability distribution between the stress state and the surrounding environment, so as to strengthen the spatial consistency constraint while retaining the local identification results, thereby eliminating pseudo-abnormal points, filling spatial fracture zones, and optimizing the structural marking of the entire stress field through the maximum a posteriori estimation method. Based on the optimized stress concentration areas, a stress concentration index evaluation system is established. This index comprehensively considers the maximum principal stress intensity, the principal stress difference, and the relative spatial relationship with geological structures (especially faults). Its formula is:
[0074] ;
[0075] where and represent the average values of the maximum principal stress and the stress difference in the area respectively, and D is the normalized distance index from the center point of this unit to the nearest active fault.
[0076] Objectively score and rank the stress concentration areas based on the SCI index, and set the grading thresholds according to the scores. For example: SCI > 4.0 is the first-level high-risk area, 2.0 ≤ SCI ≤ 4.0 is the second-level medium-risk area, and 1.0 ≤ SCI < 2.0 is the third-level area of concern. After grading, the stress concentration areas of different levels are processed separately, and different simulation weights and response priorities are assigned during subsequent modeling. Based on the above graded stress concentration areas, a three-dimensional geological stress concentration area model integrating multi-source information is constructed. In this model, the stress tensor field forms the main framework, finely expressing the spatial distribution law of the principal stresses; the lithology distribution serves as the structural background of the model, reflecting the mechanical response differences of different stratigraphic units; the degree of fracture development constructs a spatial weak surface system by superimposing geological structure feature data such as faults, fracture zones, and shear zones, providing a basis for mechanical boundaries and stress reconstruction; the groundwater distribution is embedded in the model as an important factor for fluid-induced stress field perturbation to evaluate the triple coupling relationship of water-rock-stress. To improve the model resolution, the adaptive grid meshing method is applied in areas with higher SCI values to locally refine the grid cells and enhance the modeling accuracy. The model is presented in a three-dimensional grid data structure, and the output forms include various expressions such as vector field diagrams, visualization of principal stress ellipsoids, cross-sectional diagrams, and isosurface diagrams, with spatial operability and analysis expandability.
[0077] Among them, the Markov random field model is applied to the preliminarily identified stress concentration areas for spatial constraint optimization to obtain the optimized stress concentration areas, including: dividing the preliminarily identified stress concentration areas into regular grid cells, establishing the connection relationships between each cell and its adjacent cells, constructing a Markov random field spatial network with a four-neighborhood or eight-neighborhood structure to obtain the spatial constraint network structure; defining the single-site potential function and the pairwise potential function based on the spatial constraint network structure, where the single-site potential function characterizes the conformity between the stress value of the grid cell and the observed data, and the pairwise potential function characterizes the strength of the spatial continuity constraint between adjacent cells to obtain the Markov random field energy function; using the maximum likelihood estimation or quasi-maximum likelihood estimation method, combined with the stress distribution characteristics of the known geological structure areas, to train and calculate the weight parameters in the Markov random field energy function to obtain the optimal model parameters; setting hard boundary constraints and soft transition conditions for the Markov random field model according to the geological structure unit boundaries, fault positions, and lithology boundaries to obtain the constraint conditions considering the geological structure characteristics; applying the iterative conditional mode algorithm or the graph cut algorithm to the Markov random field model with constraint conditions, and performing label optimization by minimizing the global energy function, and iteratively calculating until convergence or reaching the preset maximum number of iterations to obtain the optimized label result; post-processing the optimized label result, including removing isolated areas with an area smaller than the threshold, merging similar areas with too small spacing, smoothing the area boundaries, and verifying the result according to the principle of geomechanical equilibrium to obtain the optimized stress concentration areas that eliminate pseudo-anomalies and maintain spatial continuity.
[0078] Those skilled in the art can clearly understand that for the convenience and conciseness of description, the specific working processes of the systems, systems and units described above can refer to the corresponding processes in the foregoing method embodiments, and will not be elaborated herein.
[0079] If the integrated unit is implemented in the form of a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on such an understanding, the technical solution of the present invention, in essence, or the part that contributes to the prior art, or all or part of the technical solution, can be embodied in the form of a software product. The computer software product is stored in a storage medium and includes several instructions for causing an intelligent identification device for geological stress concentration areas (which can be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of the present invention. The foregoing storage medium includes: various media such as USB flash drives, mobile hard disks, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical discs that can store program codes.
[0080] The above is the case. The above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit them. Although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that they can still modify the technical solutions recorded in the foregoing embodiments or perform equivalent replacements for some of the technical features. These modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the various embodiments of the present invention.
Claims
1. A method for intelligently identifying geological stress concentration areas, characterized in that: include: Multiple vector electromagnetic induction coils are arranged in a grid shape in the monitoring area to collect original electromagnetic pulse signals containing three-directional components: X, Y, and Z. Performing denoising and time-frequency analysis on the original electromagnetic pulse signal to obtain a multidimensional data cube; Use the integer linear programming model to process the multidimensional data cube to obtain the preliminary optimized signal data; perform singular value decomposition and apply the Tikhonov regularization method to obtain the stabilized signal data; use the piecewise polynomial function to perform baseline correction; use the adaptive fuzzy correction algorithm to apply the Mamdani fuzzy reasoning to process the signal; use the time difference of arrival constraint condition to perform iterative reweighted least squares processing to obtain the disambiguated signal, and select the signal segments that meet the requirements according to the calculation parameters of the disambiguated signal to obtain the corrected electromagnetic pulse signal data set; Decomposing and mapping the corrected electromagnetic pulse signal data set to obtain an estimated value of geological stress distribution; A three-dimensional model of the geological stress concentration area is constructed according to the estimated value of the geological stress distribution.
2. The intelligent identification method of geological stress concentration area according to claim 1 is characterized in that: The method of arranging a plurality of vector electromagnetic induction coils in a grid-like manner in the monitoring area to collect original electromagnetic pulse signals in three directions of X, Y and Z includes: Performing a geological structure complexity analysis on the monitoring area to obtain a geological structure complexity analysis result, and determining an equilateral triangle arrangement mode according to the geological structure complexity analysis result to obtain a receiving device arrangement plan; Each receiving device in the receiving device arrangement scheme is equipped with three orthogonal vector electromagnetic induction coils to obtain an electromagnetic signal receiving unit; Integrate a GPS clock module into each electromagnetic signal receiving unit to obtain a time synchronization acquisition system, and connect the electromagnetic signal receiving units to form a mesh network according to the time synchronization acquisition system and the low-power wireless transmission module to obtain a distributed monitoring network; Selecting a reference site in the monitoring area, and recording the external electromagnetic interference signal of the area according to the reference site to obtain reference interference data; The electromagnetic pulse signal of the distributed monitoring network is monitored, and the effective signal is marked in combination with the reference interference data to obtain the original electromagnetic pulse signal containing the components in three directions of X, Y and Z.
3. The intelligent identification method of geological stress concentration area according to claim 1 is characterized in that: The performing of denoising and time-frequency analysis on the original electromagnetic pulse signal to obtain a multidimensional data cube includes: Performing interference filtering on the original electromagnetic pulse signal to obtain a preliminary filtered signal, and performing wavelet threshold denoising on the preliminary filtered signal to obtain a denoised signal; Performing time-domain segmentation processing on the noise-reduced signal to obtain a segmented time-domain signal, and performing time-frequency analysis on the segmented time-domain signal to obtain a signal feature matrix; Calculating polarization parameters including ellipticity, inclination and azimuth according to the three-component vector signal in the signal characteristic matrix, and performing adaptive notch filtering in combination with the external interference signal recorded at the reference station to obtain a pure geological source signal; The pure geological source signal is subjected to spatial correlation analysis, the mutual correlation coefficient and phase difference between adjacent measuring points are calculated, and the signals are organized in a unified format to form a multi-dimensional data cube including time, space, amplitude, frequency and phase.
4. The intelligent identification method of geological stress concentration area according to claim 1 is characterized in that: The method uses an integer linear programming model to process a multidimensional data cube to obtain preliminary optimized signal data; performs singular value decomposition and applies the Tikhonov regularization method to obtain stabilized signal data; performs baseline correction using a piecewise polynomial function; uses an adaptive fuzzy correction algorithm and applies Mamdani fuzzy reasoning to process the signal; uses a time difference of arrival constraint condition to perform iterative reweighted least squares processing to obtain an ambiguous signal, and selects signal segments that meet the requirements according to the calculation-related parameters of the ambiguous signal to obtain a corrected electromagnetic pulse signal data set, including: Inputting the multidimensional data cube into an integer linear programming model, obtaining preliminary optimized signal data by minimizing the weighted residual square sum between signal estimation values and observation values and introducing L1 norm sparsity constraints; Performing singular value decomposition on the initially optimized signal data to obtain a condition number, and applying a Tikhonov regularization method according to the condition number to obtain stabilized signal data; Performing baseline correction on the stabilized signal data using a piecewise polynomial function to obtain a baseline-corrected signal; Adopting an adaptive fuzzy correction algorithm to apply Mamdani fuzzy reasoning to the baseline-corrected signal, classifying the signal features according to the Gaussian membership function, and obtaining signal data with ambiguity eliminated; Using the time difference of arrival constraint condition, the signal data to be deambiguated is processed by iterative reweighted least square method to obtain a deambiguated signal; The signal-to-noise ratio, amplitude stability and phase coherence parameters are calculated according to the eliminated ambiguous signal, and the signal segments meeting the requirements are screened according to the quality index to obtain a corrected electromagnetic pulse signal data set.
5. The method for intelligently identifying geological stress concentration areas according to claim 1, characterized in that: Decomposing and mapping the corrected electromagnetic pulse signal data set to obtain an estimated value of geological stress distribution includes: Applying a time-redistributed multi-synchronous compression transformation algorithm to set an adaptive window length function, and dynamically adjusting the signal local frequency characteristics of the corrected electromagnetic pulse signal data set to obtain multiple narrowband single-component functions; Extracting a signal characteristic parameter set including center frequency, bandwidth, energy distribution, duration and instantaneous frequency change rate from the plurality of narrowband single-component functions; For the initial electromagnetic pulse wavefront information in the corrected electromagnetic pulse signal data set, measuring the time difference of arrival between different measuring points to obtain the time difference of arrival data; Based on the time difference of arrival data and the spatial distribution coordinates of the receiving device, an overdetermined geological stress location equation is constructed, and the electromagnetic wave propagation velocity tensor is obtained by solving the equation; An extreme gradient boosting algorithm is used to establish a mapping relationship between the signal characteristic parameter set and the geological stress parameter to obtain a geological stress prediction model; The top M features that contribute most to stress prediction are selected through feature importance analysis, and the extreme gradient boosting algorithm hyperparameters are automatically adjusted in combination with the Bayesian optimization method to optimize the geological stress prediction model and output an estimated value of geological stress distribution including principal stress direction, stress magnitude and stress gradient.
6. The intelligent identification method of geological stress concentration area according to claim 5, characterized in that: The multi-synchronous compression transformation algorithm using time redistribution is used to set an adaptive window length function, and the local frequency characteristics of the signal are dynamically adjusted for the corrected electromagnetic pulse signal data set to obtain multiple narrowband single-component functions, including: Performing a fast Fourier transform on the corrected electromagnetic pulse signal data set to obtain initial signal spectrum characteristic parameters; Based on the initial spectrum characteristic parameters of the signal, the local eigenvalues of the Hessian matrix are calculated for each frequency band, and a time-redistributed multi-synchronous compression transformation algorithm is applied in combination with the frequency change rate to formulate a window length scaling factor to obtain an adaptive window length function; Applying short-time Fourier transform to the corrected electromagnetic pulse signal data set, using the adaptive window length function to set the time window lengths of different frequency bands, and obtaining a time-frequency distribution diagram; Performing a time reallocation operation on the time-frequency distribution graph to obtain an enhanced time-frequency distribution feature; Inputting the enhanced time-frequency distribution feature into a synchronous compression kernel function to perform ridge extraction and synchronous depth control to obtain a plurality of original single-component functions; Energy contribution evaluation is performed on the multiple original single-component functions to obtain multiple narrow-band single-component functions.
7. The method for intelligently identifying geological stress concentration areas according to claim 1, characterized in that: The step of constructing a three-dimensional model of a geological stress concentration area according to the estimated value of geological stress distribution comprises: Based on the estimated value of geological stress distribution, the stress parameters of discrete measuring points are interpolated into stress field data using Kriging interpolation method; Performing tensor analysis on the stress field data to obtain principal stress distribution characteristics, and formulating a stress concentration area discrimination criterion based on the principal stress distribution characteristics to obtain a preliminarily identified stress concentration area; Applying a Markov random field model to perform spatial constraint optimization on the initially identified stress concentration area to obtain an optimized stress concentration area; Establishing a stress concentration index evaluation system based on the optimized stress concentration area, and calculating the graded stress concentration area according to the stress concentration index evaluation system; Based on the hierarchical stress set, three-dimensional modeling is performed to generate a three-dimensional model of the geological stress concentration area including stress tensor field, lithology distribution, fracture development degree and groundwater distribution information.
Citation Information
Patent Citations
Processing geophysical data
CN103038670A
Method and system for analyzing filling for karst reservoir based on spectrum decomposition and machine learning
US20230083651A1