A deep geostress prediction method and system based on seismic wave field information

Abstract

This invention relates to the fields of geophysical exploration and geomechanics, and discloses a method and system for predicting deep geostress based on seismic wavefield information, aiming to solve the problems of single information dimension and multiple solutions in existing deep geostress prediction methods. The method includes: acquiring multi-component full-wavefield seismic records and separating P-wave, fast S-wave, and slow S-wave wavefields; extracting velocity, amplitude gradient, attenuation coefficient, dominant frequency, and instantaneous phase from migration profiles to construct a multi-modal seismic wavefield attribute data volume; inputting the attribute vectors into a deep learning model containing a dual-convolutional autoencoder and a bidirectional gated recurrent unit to jointly predict the minimum and maximum horizontal principal stress gradients and the maximum principal stress azimuth; and employing a composite loss function with adaptive adjustment of multi-task weights during training. This invention enables high-precision three-dimensional prediction of deep geostress.

Description

Technical Field

[0001] This invention belongs to the fields of geophysical exploration and geomechanics, specifically relating to a method and system for predicting deep geostress based on seismic wavefield information. Background Technology

[0002] In-situ stress prediction is an indispensable fundamental technical link in oil and gas resource exploration and development, underground space engineering safety assessment, and geological hazard mechanism research. Its prediction accuracy directly affects key engineering decisions such as drilling trajectory design, reservoir fracturing effectiveness, and wellbore stability control. As oil and gas exploration and development continue to advance into deeper formations, deep and ultra-deep reservoirs are gradually becoming important areas for increasing reserves and production, placing higher demands on the coverage depth and reliability of in-situ stress prediction technology. Among various in-situ stress prediction methods, those based on seismic wavefield information have become an important approach for achieving regional-scale three-dimensional in-situ stress field prediction due to their wide coverage and ability to obtain a complete picture of the underground wavefield response.

[0003] Specifically, deep in-situ stress prediction based on seismic wavefield information typically relies on wavefield separation, parameter inversion, and stress state calculation of seismic records received at the surface or in wells. In practice, seismic data processing first extracts wavefield attributes sensitive to in-situ stress, such as P-wave and S-wave velocities, attenuation coefficients, and amplitude variations with offset. Then, a rock physics model is used to establish relationships between elastic parameters, anisotropic parameters, and the in-situ stress tensor, thereby calculating the formation stress distribution. This technical approach essentially uses the multimodal information carried by the seismic wavefield to constrain the subsurface stress state; the completeness of the information and the degree of complementary utilization between modes determine the reliability of the prediction.

[0004] However, in existing technical implementations, the utilization of deep seismic wavefield information generally remains at the level of extracting single or shallowly correlated modes, which is insufficient to support accurate solutions for the stress state of ultra-deep strata. When deep seismic waves travel through propagation paths of thousands to tens of thousands of meters, they are severely depleted of high-frequency components due to factors such as stratum absorption and attenuation, geometric diffusion, and scattering by inhomogeneous bodies. This leads to a significant decrease in the wavefield signal-to-noise ratio and resolution, resulting in the submergence of subtle wavefield disturbances sensitive to changes in ground stress. More importantly, the multi-modal information in the deep wavefield, such as P-wave velocity, S-wave velocity, attenuation factor, dispersion characteristics, and phase changes, does not exist independently; they exhibit highly coupled response characteristics under the influence of ground stress. Existing methods typically extract only one or a limited number of these modes for inversion, such as relying solely on P-wave velocity changes or pre-stack amplitude gradient characteristics to estimate stress, failing to perform synergistic fusion and correlation analysis of the aforementioned multi-modal information. This lack of information dimension results in insufficient constraint of the ground stress field on the observation data during the inversion process, and the problem of multiple solutions is drastically amplified under conditions of weak signals in deep strata. The results show that in ultra-deep strata with a burial depth of more than 3,000 meters, the predicted direction and magnitude of geostress often have great uncertainty, which cannot provide a reliable basis for drilling and completion risk warning and reservoir fracturing scheme optimization, and has become a key bottleneck restricting the efficient utilization of deep oil and gas resources. Summary of the Invention

[0005] The purpose of this invention is to overcome the shortcomings of the prior art and provide a method and system for predicting deep geostress based on seismic wavefield information, which can effectively solve the problems mentioned in the background art.

[0006] The objective of this invention is achieved through the following technical solution: a method for predicting deep geostress based on seismic wavefield information, comprising the following specific steps: S1: Deploying a multi-component seismic data acquisition and observation system on the surface or in wells of the target work area, exciting and recording a full-wavefield seismic record containing P-waves, S-waves, and converted waves through artificial seismic sources, and obtaining raw multi-component common shot gather data; S2: Preprocessing the raw multi-component common shot gather data and performing multi-domain wavefield separation to obtain time-domain common reflection gathers corresponding to the P-wave wavefield, fast S-wave wavefield, and slow S-wave wavefield, and then performing pre-stack time migration processing on the separated wavefield data to obtain P-wave migration profiles, fast S-wave migration profiles, and slow S-wave migration profiles; S3: From the P-wave... In migration profiles, fast S&W migration profiles, and slow S&W migration profiles, P-wave velocity, fast S&W velocity, slow S&W velocity, P-wave reflection amplitude gradient, fast S&W reflection amplitude gradient, slow S&W reflection amplitude gradient, P-wave attenuation coefficient, fast S&W attenuation coefficient, slow S&W attenuation coefficient, P-wave dominant frequency, fast S&W dominant frequency, slow S&W dominant frequency, P-wave instantaneous phase, fast S&W instantaneous phase, and slow S&W instantaneous phase are extracted along the target layer to construct a multimodal seismic wavefield attribute data volume oriented towards the target layer. Each spatial sampling point in the multimodal seismic wavefield attribute data volume corresponds to a multidimensional attribute vector, and the components of this attribute vector are, in order, P-wave velocity, fast S&W velocity, slow S&W velocity, P-wave reflection amplitude gradient, and fast S&W reflection amplitude gradient. S4: Construct a deep learning model for predicting deep geostress. The input layer of the deep learning model for predicting deep geostress receives the multidimensional attribute vectors. Its hidden layer contains a first convolutional autoencoder, a second convolutional autoencoder, a bidirectional gated recurrent unit network, and a fully connected regression layer, which are cascaded together. Its output layer outputs the minimum horizontal principal stress gradient value, the maximum horizontal principal stress gradient value, and the maximum horizontal principal stress azimuth angle. S5: Utilize the measured geostress data at known well points and the multidimensional attribute vectors at corresponding spatial sampling points. A training sample set is constructed using attribute vectors to supervise the training of the deep learning model for deep stress prediction. During training, a composite loss function with adaptive adjustment of multi-task weights is used to constrain the updating of model parameters until the composite loss function converges to a preset threshold, thus obtaining the trained deep learning model for deep stress prediction. S6: The multidimensional attribute vectors of all spatial sampling points in the multimodal seismic wavefield attribute data volume of the target layer to be predicted are sequentially input into the trained deep learning model for deep stress prediction. Through forward propagation calculation, the minimum horizontal principal stress gradient data volume, the maximum horizontal principal stress gradient data volume, and the maximum horizontal principal stress azimuth angle data volume of the entire target layer are obtained, realizing the three-dimensional prediction of deep stress.

[0007] Preferably, in step S1, the multi-component seismic data acquisition and observation system adopts a linear or grid-like layout with equal or non-equal spacing. The single-channel detector is a three-component digital detector, whose three components record the particle vibration velocity or acceleration of the vertical component, the east-west horizontal component, and the north-south horizontal component, respectively. The artificial seismic source is an explosive source or a controllable source. The frequency range of the elastic waves generated by the excitation covers the preset low-frequency to high-frequency range. The recording length is a preset recording duration, and the time sampling interval meets the preset high-precision acquisition requirements.

[0008] Preferably, in step S2, the preprocessing includes sequentially performing observation system definition, trace editing, static correction, spherical diffusion compensation, surface uniformity amplitude compensation, deconvolution, and noise suppression on the original multi-component common shot gather data; the multi-domain wavefield separation utilizes a combination of polarization filtering and rotational projection to project the original three-component seismic record from the Cartesian coordinate system composed of the vertical component, the east-west horizontal component, and the north-south horizontal component to the ray coordinate system, and separates the P-wave field, fast S-wave field, and slow S-wave field by analyzing the polarization ellipse parameters of the particle vibration trajectory, and performs up-and-downward wave separation and S-wave static correction on the separated S-wave field.

[0009] Preferably, in step S3, the P-wave velocity, fast S-wave velocity, and slow S-wave velocity are obtained through pre-stack time migration velocity analysis along the target layer. The P-wave reflection amplitude gradient, fast S-wave reflection amplitude gradient, and slow S-wave reflection amplitude gradient are obtained by fitting the amplitude variation curve of the target layer reflected wave in their respective migration profiles with the migration distance and calculating its slope. The P-wave attenuation coefficient, fast S-wave attenuation coefficient, and slow S-wave attenuation coefficient are obtained by calculating the dominant frequency migration of the top and bottom reflected waves of the target layer using the spectral ratio method and combining it with the interlayer travel time. The P-wave dominant frequency, fast S-wave dominant frequency, slow S-wave dominant frequency, P-wave instantaneous phase, fast S-wave instantaneous phase, and slow S-wave instantaneous phase are obtained by performing Fourier transform and Hilbert transform on the records within the time window of the target layer in their respective migration profiles.

[0010] Preferably, in step S4, the first convolutional autoencoder has a preset first number of convolutional layers and deconvolutional layers, used for feature extraction and reconstruction pre-training of the input multidimensional attribute vector, and its encoder outputs a first intermediate feature vector; the second convolutional autoencoder has a preset second number of convolutional layers and deconvolutional layers, which receives the first intermediate feature vector and further extracts deep abstract features, and its encoder outputs a second intermediate feature vector; the bidirectional gated recurrent unit network includes a forward propagation layer and a backward propagation layer, with the number of hidden units in each layer being a preset number of units, the bidirectional gated recurrent unit network receives the second intermediate feature vector, captures the spatial correlation and mutual constraint relationship between multimodal attributes through bidirectional temporal modeling, and outputs a temporal fusion feature vector; the fully connected regression layer includes a hidden layer with a preset first number of neurons and an output layer with a preset second number of neurons, which receives the temporal fusion feature vector and regresses to output the minimum horizontal principal stress gradient value, the maximum horizontal principal stress gradient value, and the maximum horizontal principal stress azimuth angle.

[0011] Preferably, in step S5, the composite loss function consists of a mean squared error loss term and a multi-task weight adaptive adjustment term. The mean squared error loss term is a weighted sum of the predicted mean squared error of the minimum horizontal principal stress gradient, the predicted mean squared error of the maximum horizontal principal stress gradient, and the predicted mean squared error of the maximum horizontal principal stress azimuth angle. The multi-task weight adaptive adjustment term dynamically adjusts the weight coefficients of different task losses according to the reciprocal of the loss decrease rate of each task on the validation set, so that the optimization step size of each task remains consistent during model training, avoiding the training of a certain task from dominating the model update direction due to differences in units and numerical scales.

[0012] Preferably, in step S5, the measured in-situ stress data at the known well points are obtained through any of the following methods: acoustic emission Kaiser effect test, used to obtain the memory information of the core under historical tectonic stress; differential strain analysis (DSA), to obtain triaxial in-situ stress values ​​by analyzing the strain recovery characteristics of directional cores after unloading; hydraulic fracturing method, to determine the closure pressure and re-tension pressure when the well wall ruptures through the fracturing curve, and to calculate the minimum horizontal principal stress. The training sample set is constructed as follows: for each known well, the well trajectory is spatially matched with the multimodal seismic wave field attribute data volume, and the multidimensional attribute vector of the target layer at the well point is extracted as the input sample. The measured minimum horizontal principal stress gradient value, the measured maximum horizontal principal stress gradient value, and the measured maximum horizontal principal stress azimuth angle of the corresponding depth segment are used as labels to form a training sample. All training samples corresponding to all known wells are summarized into the training sample set.

[0013] Preferably, in step S3, before inputting the multidimensional attribute vector of each spatial sampling point in the multimodal seismic wavefield attribute data volume into the deep learning model for deep geostress prediction, the multidimensional attribute vector is standardized using the statistical features of all attribute vectors in the training sample set, so that each dimension of the feature follows a preset statistical distribution. When the deep learning model for deep geostress prediction is trained and predicted, the input data are the standardized multidimensional attribute vectors.

[0014] The present invention also provides a deep geostress prediction system based on seismic wavefield information, the system performing all the steps included in the deep geostress prediction method based on seismic wavefield information.

[0015] The present invention has the following advantages:

[0016] This invention systematically acquires and processes full-wavefield multi-component seismic data, extracting wavefield information encompassing multiple dimensions and attributes, including velocity, amplitude gradient, attenuation coefficient, dominant frequency, and instantaneous phase, from migration profiles of P-waves and two types of S-waves. This constructs a multimodal seismic wavefield attribute data volume that comprehensively characterizes the deep geostress response, overcoming the limitations of existing technologies that rely on only a single or limited information mode. Based on this, a deep learning model for deep geostress prediction integrating a dual-convolutional autoencoder and a bidirectional gated recurrent unit network is constructed. The convolutional autoencoder abstracts and extracts nonlinear high-order features of multimodal attributes layer by layer, and the bidirectional gated recurrent unit network performs temporal modeling of the spatial continuity and mutual constraints between features. Ultimately, it achieves joint high-precision prediction of the minimum horizontal principal stress gradient, maximum horizontal principal stress gradient, and azimuth in ultra-deep strata. The multi-task weight adaptive adjustment composite loss function introduced during training intelligently balances the contribution of tasks with different physical dimensions in optimization, significantly improving the model's convergence stability and generalization ability under deep, weak signal conditions. This invention does not rely on the assumption of a specific reservoir type and is adaptable to various complex geological conditions. In particular, it can still maintain reliable prediction capabilities for ultra-deep strata with a burial depth exceeding the preset depth, thus solving the key technical bottlenecks of single information dimension and multiple solutions amplification in deep geostress prediction. Attached Figure Description

[0017] Figure 1 This is a schematic diagram of the overall technical architecture of the deep geostress prediction method based on seismic wavefield information proposed in this invention.

[0018] Figure 2 This is a schematic diagram of the core principle framework of the deep learning model for predicting deep geostress in this invention. Detailed Implementation

[0019] The present invention will be further described below with reference to the accompanying drawings, but the scope of protection of the present invention is not limited to the following description.

[0020] Example 1:

[0021] This embodiment is specifically applied to the scenario of geostress prediction in deep, tight sandstone oil and gas reservoirs on land. The target strata are buried at depths between 4500 and 6500 meters, and the structural background is a steep anticline at the front of a thrust-nappe zone. The geostress distribution exhibits strong heterogeneity, and the direction of the horizontal principal stress varies significantly with the structural location. In this embodiment, by deploying a high-density multi-component seismic data acquisition and observation system and constructing a deep learning model for deep geostress prediction that integrates a dual-convolutional autoencoder and a bidirectional gated recurrent unit network, a three-dimensional joint prediction of the minimum horizontal principal stress gradient, the maximum horizontal principal stress gradient, and the azimuth of the maximum horizontal principal stress in the ultra-deep target strata is systematically achieved.

[0022] like Figure 1 As shown, at the system architecture level, the deep geostress prediction system based on seismic wavefield information adopted in this embodiment consists of six core functional modules: a data acquisition module, a preprocessing and wavefield separation module, a wavefield attribute extraction module, a deep learning model and prediction module, a model training and deployment module, and a final output and visualization module. These modules are interconnected via a high-speed data transmission bus and a shared storage array, and are deployed as a whole within an integrated exploration geophysical high-performance computing cluster environment. The data acquisition module is located at the forefront of the entire system and undertakes the physical implementation task of step S1. Specifically, it consists of an array of field surface seismic acquisition equipment and a downhole three-component geophone array. The field surface seismic acquisition equipment array adopts a linear layout, with each receiving line containing 128 acquisition channels. The distance between receiving lines is 400 meters, and the distance between receiving channels is 50 meters, forming a combined regularized observation grid of 16 receiving lines and a total of 2048 channels. Each acquisition channel is equipped with a three-component digital detector, which integrates three orthogonally arranged MEMS accelerometer sensor elements. The vertical component MEMS accelerometer is mounted vertically to record the projected value of the particle vibration velocity in the vertical direction; the east-west horizontal component MEMS accelerometer is mounted east-west oriented to record the east-west horizontal component of the particle vibration velocity; and the north-south horizontal component MEMS accelerometer is mounted north-south oriented to record the north-south horizontal component of the particle vibration velocity. The MEMS accelerometers are factory-calibrated for sensitivity, with a voltage sensitivity of 2.0V / g, a dynamic range greater than or equal to 120dB, a frequency response flatness deviation of no more than ±3dB within the range of 2Hz to 800Hz, and a self-noise level below 1. Each three-component digital detector is equipped with a 24-bit Σ-Δ analog-to-digital converter, with a uniform sampling rate set to 1000Hz and a 1ms sampling interval. The quantized digital signal is transmitted to the acquisition station as a differential signal via armored twisted-pair cable, with a differential signal swing of ±2.5V and a common-mode rejection ratio greater than 100dB. The acquisition station integrates a programmable gain amplifier, digital filters, and a data packetizer. The gain level of the programmable gain amplifier is adjustable in 6dB steps from 0dB to 36dB, depending on the source energy and detector sensitivity. The digital filter uses a linear-phase FIR low-pass filter with a cutoff frequency set to 400Hz to filter out noise above the useful signal frequency band. The data packetizer packages the three component data from each sampling point into a 72-bit data frame, with each component occupying 24 bits. An 8-bit acquisition channel identifier and a 32-bit GPS timestamp are appended to the frame header, and the data is transmitted to the central recording system via fiber optic cable through a 100Mbps Ethernet physical layer interface. The central recording system is deployed in a vehicle-mounted data processing container. Its core is a rack-mount server equipped with two 20-core CPUs (2.3GHz), 256GB of DDR4 ECC memory with cache, and a RAID 6 array of twelve 4-terabyte SAS hard drives, providing approximately 40 terabytes of usable storage. This server connects to a seismic data acquisition interface unit via a PCIe 3.0 x8 bus interface. The interface unit aggregates and deframes data streams from multiple acquisition stations transmitted via fiber optic cables and writes them to the disk array in a common shot point gather format.

[0023] The artificial seismic source excitation module employs three articulated controllable seismic source vehicles, each with a peak output of 280 kN. The excitation scanning signal is a linear frequency-converted sinusoidal signal, with a scanning frequency range of 2 Hz to 120 Hz. The scanning length is 20 seconds, the recording length is 10 seconds, and a 5-second listening time is reserved after the scan, for a total recording length of 35 seconds. The hydraulic servo control system of the controllable seismic source acquires the acceleration signals of the flat plate and the reaction mass in real time. Through a PID control algorithm, the phase difference between the output force signal and the reference scanning signal is controlled within ±2 degrees, and the amplitude deviation is controlled within ±3%. The central recording system receives the GPS synchronization trigger signal and scanning status parameters from the controllable seismic source control system via a data transmission radio. The timestamp accuracy of the trigger signal is better than 10 microseconds, ensuring the time synchronization accuracy of all acquisition channels with the seismic source excitation moment.

[0024] The preprocessing and wavefield separation module is deployed on the same rack-mounted server as the central recording system. This server is also equipped with four additional NVIDIA Tesla V100 GPU accelerator cards, each with 32GB of HBM2 memory and a memory bandwidth of 900GB / s. The GPUs are connected via an NVLink high-speed interconnect bus with a bidirectional bandwidth of 300GB / s, used to accelerate wavefield separation and migration processing calculations for large datasets. At the software level, this module consists of an observation system definition submodule, a trace editing and static correction submodule, an amplitude compensation and deconvolution submodule, a polarization filtering and wavefield separation submodule, and a pre-stack time migration submodule. These components are integrated in series as seismic data processing middleware. The observation system definition submodule reads the GPS coordinate information and acquisition trace identification codes from the raw multi-component common-shot-point gather data. Combined with the survey line stationings and receiver elevation data obtained from the surveying project, it constructs a spatial location relationship database, assigns the receiver coordinates and shot point coordinates of each seismic trace to a unified Cartesian coordinate system, and calculates the shot-receiver distance, azimuth, and center point coordinates, generating an observation system spatial attribute table. The trace editing and static correction submodule performs bad trace removal, polarity reversal correction, and first-arrival refraction static correction on seismic traces. Bad trace removal involves statistically analyzing the root mean square amplitude value of each seismic trace, marking traces with amplitude values ​​exceeding three standard deviations of the average root mean square amplitude across the entire work area as bad traces and filling them with zeros. Polarity reversal correction aims to unify the waveform polarity inconsistencies caused by different detector and source combinations. First-arrival refraction static correction utilizes the first-arrival travel-time tomography method to invert the near-surface velocity model, calculating the static correction amount at the shot and receiver points. The correction time length is within the range of 0 to 120 ms, and a full-trace time shift operation is performed when applying static correction. The amplitude compensation and deconvolution submodule sequentially performs spherical diffusion compensation, surface consistency amplitude compensation, and predictive deconvolution. The gain function of spherical diffusion compensation is a first-order gain over time, and the gain constant is calculated based on the energy of the reflected wave at the target layer and the theoretical attenuation of spherical diffusion. After compensation, the energy of the effective reflected wave at deeper layers is increased. Surface uniform amplitude compensation decomposes the amplitude of each seismic trace into the product of four factors: shot point term, receiver term, shot-receiver offset term, and common center point term. Gauss-Seidel iteration is used to solve for each factor. After 20 iterations, each factor tends to stabilize, and the energy of each trace tends to be spatially balanced after compensation. Predictive deconvolution uses the autocorrelation method to estimate the seismic wavelet, with a prediction step size set to four times the sampling point (4 ms). The operator length is 120 ms, and Wiener filtering is used to eliminate seismic wavelet sidelobes and improve resolution.

[0025] The polarization filtering and wavefield separation submodule is the core component of the preprocessing and wavefield separation module. It is responsible for transforming the original three-component seismic record from a Cartesian coordinate system composed of the vertical, east-west, and north-south horizontal components to a ray coordinate system via rotational projection. First, a covariance matrix is ​​constructed using the particle vibration trajectories within the first arrival time window. This covariance matrix is ​​a 3×3 real symmetric matrix, with its elements corresponding to the zero-delay cross-correlation values ​​of the three vibration signal components. Eigenvalue decomposition is performed on this covariance matrix, yielding three eigenvalues ​​and three mutually orthogonal eigenvectors. The direction of the eigenvector corresponding to the largest eigenvalue is the principal polarization direction of the particle vibration, while the direction of the eigenvector corresponding to the smallest eigenvalue is perpendicular to the principal polarization plane. The wave type is determined based on the angle between the principal polarization direction and the assumed ray direction. If the angle is less than 15 degrees, the wavefield within that time window is determined to be a P-wave; if the angle is greater than 75 degrees, it is determined to be a S-wave. For the portion identified as a shear wave, it is further projected onto a plane perpendicular to the ray direction. A 2×2 covariance matrix is ​​then constructed in this plane, and eigenvectors are calculated. The direction with the larger eigenvalue is the fast shear wave polarization direction, and the direction with the smaller eigenvalue is the slow shear wave polarization direction. Through the above polarization filtering and rotation projection operations, the original three-component records are projected onto the P-wave polarization direction, the fast shear wave polarization direction, and the slow shear wave polarization direction, respectively, obtaining the time-domain common reflection point gathers for the P-wave field, the fast shear wave field, and the slow shear wave field. The separated shear wave field also needs to undergo uplink and downlink wave separation processing. This processing uses FK filtering in the frequency-wavenumber domain to distinguish between uplink and downlink waves based on the sign of the apparent velocity, retaining only the uplink reflected shear wave. Subsequently, shear wave static correction was performed on the shear wave field. The shear wave static correction amount was obtained by multiplying the P-wave static correction amount by an empirical coefficient, which was selected based on the near-surface P-wave velocity ratio in the range of 2.0 to 4.0, and was determined by cross-validation and fine adjustment with the shear wave first arrival traveltime tomography results.

[0026] The pre-stack time migration submodule performs Kirchhoff integral method pre-stack time migration processing on the time-domain common reflection point gathers of the separated P-wave, fast S-wave, and slow S-wave fields, respectively, to obtain P-wave migration profiles, fast S-wave migration profiles, and slow S-wave migration profiles. The core of Kirchhoff integral method pre-stack time migration is amplitude-weighted superposition along the diffraction travel time curve. The travel time calculation is based on the root mean square velocity model, which is obtained through iterative analysis of pre-stack time migration velocity along the target layer. Velocity analysis is performed in interactive velocity picking software, selecting velocity spectra at multiple control points. The velocity spectrum calculation adopts a velocity scanning method based on the Semblance coherence metric, with a scanning interval of 20 m / s. The root mean square velocity corresponding to the peak energy cluster is manually picked, converted into layer velocity using the Dix formula, and an initial velocity model is constructed. After three iterations of migration and residual velocity analysis, the velocity model converges, and the continuity of the phase axis of the reflected wave from the target layer on the migration profile reaches the optimal level.

[0027] After obtaining the three-dimensional spatially covered P-wave migration profile, fast S-wave migration profile, and slow S-wave migration profile, the wavefield attribute extraction module performs stratigraphic tracing along the interpretation stratigraphic surface of the target stratigraphic level, one by one, along the main survey lines and connecting survey lines, generating spatial surface meshes for the top and bottom interfaces of the target stratigraphic level. Within the stratigraphic segment between the top and bottom interfaces, spatial sampling points are defined with a transverse grid spacing of 25 m × 25 m, with each spatial sampling point corresponding to an attribute extraction location. For each spatial sampling point, 15 attribute values ​​are simultaneously extracted from the P-wave migration profile, fast S-wave migration profile, and slow S-wave migration profile. Among them, the P-wave layer velocity, fast S-wave layer velocity, and slow S-wave layer velocity are obtained by transforming the root mean square velocity model that finally converged in the pre-stack time migration velocity analysis process using the Dix formula. During the transformation, the velocity distortion caused by the dip angle of the layer interface is considered, and the Dix formula is dip-corrected. The P-wave reflection amplitude gradient, fast S-wave reflection amplitude gradient, and slow S-wave reflection amplitude gradient are obtained by fitting the amplitude variation curves of the reflected waves from the target layer in their respective migration profiles with the migration distance. Specifically, within the time window of the phase axis of the reflected waves from the top interface of the target layer, the root mean square amplitude values ​​corresponding to different migration distances are extracted. A linear fit is performed with the square of the migration distance as the abscissa and the root mean square amplitude value as the ordinate. The least squares method is used for fitting, and the slope of the fitted line is the reflection amplitude gradient. If the correlation coefficient is lower than 0.6, the gradient at that sampling point is marked as a null value and supplemented by subsequent interpolation. The P-wave attenuation coefficient, fast S-wave attenuation coefficient, and slow S-wave attenuation coefficient are calculated using the spectral ratio method. Short-time Fourier transforms were performed on the reflected waves from the top and bottom interfaces of the target layer, with a time window length of 80 ms. The dominant frequency and logarithmic amplitude spectrum within the effective frequency band were extracted. The ratio of the logarithmic amplitude spectrum to the frequency was linearly fitted, with the slope proportional to the reciprocal of the interlayer quality factor Q. The attenuation coefficient was then obtained by combining this with the interlayer travel time calculated from the layer velocity and layer thickness. The dominant frequencies of the P-wave, fast S-wave, and slow S-wave were obtained by performing Fourier transforms on the records within the time window of the target layer in their respective migration profiles, calculating the centroid frequencies of their amplitude spectra. The instantaneous phases of the P-wave, fast S-wave, and slow S-wave were obtained by performing Hilbert transforms on the records within the same time window, calculating the instantaneous phase of the analytic signal. The Hilbert transform was performed by multiplying the frequency domain by a sign function and then inversely transforming back to the time domain. To avoid endpoint effects, a 10 ms cosine attenuation border was applied at both ends of the time window. The 15 attribute values ​​extracted from each spatial sampling point are organized into a 15-dimensional attribute vector in sequence. The components of this attribute vector are, in order: P-wave velocity, fast S-wave velocity, slow S-wave velocity, P-wave reflection amplitude gradient, fast S-wave reflection amplitude gradient, slow S-wave reflection amplitude gradient, P-wave attenuation coefficient, fast S-wave attenuation coefficient, slow S-wave attenuation coefficient, P-wave dominant frequency, fast S-wave dominant frequency, slow S-wave dominant frequency, P-wave instantaneous phase, fast S-wave instantaneous phase, and slow S-wave instantaneous phase.The 15-dimensional attribute vectors of all spatial sampling points in the entire work area together constitute the multimodal seismic wavefield attribute data volume oriented towards the target layer. The attribute data volume is stored in HDF5 format and is internally organized as a three-dimensional regular grid with 800 points in the transverse survey line direction, 1200 points in the longitudinal survey line direction, and 1 layer in the time or depth direction.

[0028] The deep learning model and prediction module is the core unit for completing steps S4 and S6 of the entire system. The hardware platform for this module is a dedicated deep learning inference server, configured with two Intel Xeon Platinum 8268 processors, totaling 48 physical cores, with a clock speed of 2.9GHz, 384GB DDR4 ECC memory cache, and eight NVIDIA A100 GPU accelerator cards. Each A100 has 40GB HBM2e video memory, with a memory bandwidth of 1555GB / s. The GPUs are connected via a third-generation NVSwitch full interconnect structure, with a bidirectional total bandwidth of 600GB / s, supporting direct access from a single GPU to the video memory of other GPUs. This server connects to the shared storage array of the central recording system via a 100GB Ethernet card, reading the multimodal seismic wavefield attribute data files stored therein. Figure 2As shown, at the software level, the deep learning model and prediction module load a pre-trained deep learning model for predicting deep geostress. The model's network structure consists of an input layer, a first convolutional autoencoder, a second convolutional autoencoder, a bidirectional gated recurrent unit network, a fully connected regression layer, and an output layer, cascaded sequentially. The input layer receives a 15-dimensional attribute vector after Z-score normalization, with a data dimension of 15. The input layer does not transform the data and directly feeds it into the first convolutional autoencoder. The first convolutional autoencoder consists of an encoder part and a decoder part. The encoder part contains three convolutional blocks, each consisting of a one-dimensional convolutional layer, a batch normalization layer, and a leaky rectified linear unit activation function connected sequentially. The first convolutional block has a kernel size of 3, 32 kernels, a stride of 2, and uses "same" padding, resulting in an output feature map size of 8×32. The second convolutional block has a kernel size of 3, 64 kernels, a stride of 1, and uses "same" padding, resulting in an output feature map size of 8×64. The third convolutional block has a kernel size of 3, 128 kernels, a stride of 1, and uses "same" padding, resulting in an output feature map size of 8×128. The encoder's final output is flattened to obtain a 1024-dimensional first intermediate feature vector. The decoder part of the first convolutional autoencoder is strictly symmetrical to the encoder, containing three deconvolutional blocks that progressively reconstruct the 1024-dimensional first intermediate feature vector back into a 15-dimensional input attribute vector. During model training, the first convolutional autoencoder undergoes separate pre-training using mean squared error as the loss function, optimizing only the encoder and decoder parameters without affecting subsequent network layers. The initial learning rate during the pre-training phase was set to 0.001, using the Adam optimizer and a batch size of 64. After 200 training epochs, the loss value decreased to the order of 10⁻⁴ and remained stable. After pre-training, the parameters of the encoder part of the first convolutional autoencoder were fixed, and its output 1024-dimensional first intermediate feature vector was used as the input of the second convolutional autoencoder.

[0029] The second convolutional autoencoder also adopts a symmetrical encoding and decoding structure. The encoder part contains four convolutional blocks with a progressive structure. The first convolutional block has a kernel size of 3, 128 kernels, a stride of 1, and "same" padding, with an output feature map size of 1024×128. The second convolutional block has a kernel size of 3, 256 kernels, a stride of 2, and "same" padding, with an output feature map size of 512×256. The third convolutional block has a kernel size of 3, 256 kernels, a stride of 1, and "same" padding, with an output feature map size of 512×256. The fourth convolutional block has a kernel size of 3, 512 kernels, a stride of 1, and "same" padding, with an output feature map size of 512×512. The encoder output is compressed into a 512-dimensional second intermediate feature vector through global average pooling. This second intermediate feature vector highly abstractly encodes the nonlinear high-order coupling features and correlation patterns between multimodal attributes caused by geostress. The second convolutional autoencoder also employs mean squared error loss and the Adam optimizer during the pre-training phase, using the 1024-dimensional first intermediate feature vector as the reconstruction target. The training runs for 150 epochs. After pre-training, the encoder parameters are fixed, and the 512-dimensional second intermediate feature vector is fed into a bidirectional gated recurrent unit network.

[0030] The bidirectional gated recurrent unit (ROU) network consists of parallel forward propagation gated ROU layers and backward propagation gated ROU layers, each containing 128 hidden units. The internal structure of each gated ROU includes two gating mechanisms: a reset gate and an update gate. The reset gate controls the degree to which the hidden state information from the previous time step is forgotten, while the update gate controls the mixing ratio between the current candidate hidden state and the hidden state from the previous time step. The second intermediate feature vector received by the bidirectional gated ROU network is interpreted in sequence dimension as a feature sequence of spatially sampled points along the survey line direction or the main structural direction. The forward propagation layer processes the sequence along the increasing direction of spatially sampled points, and the backward propagation layer processes the sequence along the decreasing direction of spatially sampled points. The hidden states of the last time step in both the forward and backward directions are concatenated into a 256-dimensional temporal fusion feature vector. This bidirectional temporal modeling mechanism captures the continuous variation patterns and mutual constraints of multimodal properties between adjacent sampling points along the spatial direction. Examples include the covariant trend of increasing velocity while decreasing attenuation coefficient, and the discontinuous variation patterns of instantaneous phase in tectonic anomaly zones. These temporal patterns reflect the wavefield response to changes in geostress direction and the concentration or release of its magnitude. The spliced ​​temporal fusion feature vector is then input into the fully connected regression layer.

[0031] The fully connected regression layer consists of two fully connected sublayers. The first fully connected sublayer contains 128 neurons, employs a leaky rectified linear unit activation function with a leakage coefficient set to 0.1, and performs a nonlinear mapping on the input 256-dimensional temporal fusion feature vector. The second fully connected sublayer, the output layer, contains 3 neurons and does not use an activation function. It directly regresses the minimum horizontal principal stress gradient (in kPa / m), the maximum horizontal principal stress gradient (also in kPa / m), and the maximum horizontal principal stress azimuth (in degrees, measured clockwise from true north). The three neurons in the output layer correspond to these three prediction tasks.

[0032] In the model training and deployment module, step S5 uses measured in-situ stress data at known well points and 15-dimensional attribute vectors at corresponding spatial sampling points to construct a training sample set, performing end-to-end supervised training on the aforementioned deep learning model for deep in-situ stress prediction. In this embodiment, measured in-situ stress data from five exploration wells within the work area were collected. For three wells, the minimum horizontal principal stress was obtained using hydraulic fracturing, and the maximum horizontal principal stress azimuth was determined by combining the wellbore collapse azimuth and the drilling-induced fracture azimuth. For the other two wells, triaxial in-situ stress values ​​and azimuths were obtained using acoustic emission Kaiser effect tests and differential strain analysis. During sample construction, the well trajectory coordinates of each well were aligned with the spatial coordinate system of the multimodal seismic wavefield attribute data volume. The intersection points of the well trajectory and the target layer were matched to grid nodes according to the nearest neighbor principle. The 15-dimensional attribute vector at each node was extracted, and the average measured in-situ stress value within a 10-meter range above and below the corresponding depth was selected as the label. A total of 275 valid training samples were obtained from the five wells. To expand the effective sample size and enhance the model's generalization ability, a data augmentation strategy based on physical constraints was applied to the training samples: Gaussian noise was added to the 15-dimensional attribute vector of each existing sample, with the noise standard deviation being 5% of the standard deviation of each attribute in all samples, generating 5 augmented samples. Simultaneously, random perturbations of no more than ±2% were added to the label values. After augmentation, the total number of training samples reached 1650, of which 80% (1320 samples) were randomly selected as the training set, and the remaining 20% ​​(330 samples) were used as the validation set.

[0033] Before training, the 15-dimensional attribute vectors of all spatial sampling points in the multimodal seismic wavefield attribute data volume were Z-score standardized. The mean and standard deviation of each dimension of the attribute vectors of all 1650 samples in the training sample set were calculated. These statistics were used to standardize the entire attribute data volume, ensuring that the features input to the deep learning model follow a distribution with a mean of 0 and a standard deviation of 1. The standardized training set was then fed into the model for supervised training. A composite loss function with adaptive adjustment terms for multi-task weights was used during training. The complete expression of this composite loss function is as follows:

[0034]

[0035] Where MSE1, MSE2, and MSE3 are the prediction mean square errors of the minimum horizontal principal stress gradient, the maximum horizontal principal stress gradient, and the maximum horizontal principal stress azimuth angle, respectively; w1, w2, and w3 are static weight coefficients initially set to 1 / 3 during the initial training phase, and subsequently adjusted by the multi-task weight adaptive adjustment term L. adapt Dynamic correction. The multi-task weight adaptive adjustment term dynamically adjusts the weight coefficients of different task losses based on the reciprocal of the loss decrease rate of each task on the validation set. Every 10 training epochs, the decrease in loss of each task on the validation set relative to the previous 10 epochs is calculated. Tasks with larger decreases indicate faster learning rates, and their corresponding weight coefficients are moderately reduced; tasks with smaller decreases indicate learning lag, and their weight coefficients are moderately increased. The adjustment step size is controlled by the coefficient λ, which is initially set to 0.01. This mechanism allows the prediction tasks with three different physical dimensions and numerical scales to be optimized at a more balanced pace during model training, avoiding the situation where the minimum and maximum horizontal principal stress gradient values ​​typically vary greatly within the range of 15 to 25 kPa / m, and the maximum horizontal principal stress azimuth angle is within the range of 0 to 180 degrees and exhibits angular periodicity, causing one task to dominate the global gradient update direction.

[0036] The training process uses the AdamW optimizer with an initial learning rate of 0.0005, a weight decay coefficient of 0.0001, a batch size of 32, and a total of 800 training epochs. In each epoch, the training set loss and validation set loss are calculated. An early stopping mechanism is triggered when the validation set loss does not significantly decrease for 50 consecutive epochs, and the model parameters with the lowest validation set loss are saved. After 750 training epochs, the validation set loss converges to a preset threshold of 3.5 × 10⁻⁻⁻⁻⁶. 4 The root mean square error is approximately 0.0187, indicating that model training is complete. Training took approximately 4.2 hours, with peak GPU memory usage of approximately 28GB.

[0037] At the workflow level, the deep geostress prediction system in this embodiment strictly follows steps S1 to S6 in chronological order. Each step is triggered by the output data of the previous step, forming an automated data processing and prediction pipeline. Step S1 is executed by the data acquisition module. At the field construction site, surveyors use a real-time dynamic differential GPS positioning instrument to precisely set the acquisition track station numbers according to the previously designed observation grid coordinates. 2048 three-component digital geophones are buried according to the station numbers, with the geophone tail cones vertically inserted into the dense soil layer 30 cm below the surface, ensuring tight coupling between the geophones and the ground, and the verticality is corrected using a horizontal bubble. All geophones are connected to the acquisition station via armored twisted-pair cables in a tree topology. Every four receiving lines share a cross station, which aggregates data from its respective acquisition station and transmits the data stream to the central recording system through a fiber optic backbone loop. Three articulated controllable seismic source vehicles sequentially vibrate and excite at preset shot locations, with a source spacing of 50 meters, generating a total of 1250 shots, covering an area of ​​approximately 80 square kilometers. Before each firing of the controllable seismic source, its GPS receiver performs real-time dynamic positioning and calculates the deviation from the designed shot point, with a positioning accuracy better than 0.3 meters. The control software of the central recording system starts data acquisition 500ms before the scheduled firing time, and all acquisition channels simultaneously enter the recording state. The sampling clock is locked to the GPS second pulse to ensure that the time synchronization accuracy between acquisition channels reaches the microsecond level. The elastic wave frequency generated by the seismic source covers a wide frequency range from 2Hz to 120Hz. The recording time is 35 seconds, the time sampling interval is 1ms, the data volume of a single shot is approximately 215 megabytes, and the total data volume of all 1250 shots is approximately 260GB, which is written to the disk array in SEG-Y format.

[0038] Step S2 involves the preprocessing and wavefield separation module processing the original multi-component common shot gather data. First, the observation system definition submodule reads the shot coordinates, receiver coordinates, and trace head information from the SEG-Y file header, constructs a spatial attribute database, and completes the observation system definition. The trace editing and static correction submodule automatically detects bad traces. In this work area, the proportion of bad traces is approximately 1.8%, mostly concentrated in the hillside areas with dramatic topographic relief. After removing bad traces, adjacent valid traces are reconstructed using local radial basis function interpolation. Static correction employs the refracted first-arrival travel-time tomography method to invert the near-surface velocity model. The inversion grid has a horizontal dimension of 50m × 50m and a vertical layer thickness of 10m. After 8 iterations, the root mean square value of the travel-time residual converges to 3.2ms. Based on this, the static correction amounts for the shot and receiver points are calculated and applied. The amplitude compensation and deconvolution submodule sequentially performs spherical diffusion compensation, surface consistency amplitude compensation, and predictive deconvolution. After deconvolution, the seismic dominant frequency increases from 35Hz to approximately 55Hz, significantly improving the vertical resolution. The polarization filtering and wavefield separation submodule performs wavefield separation on the three-component data after amplitude compensation. For the common-shot gathers of each shot, a time window from 100ms before the first arrival to 300ms after arrival is selected to construct the polarization covariance matrix and perform eigenvalue decomposition. Based on the angle between the polarization direction and the ray direction, P and S waves are separated. Within the S wave field, fast S and slow S waves are further distinguished, resulting in three independent common-reflection gathers. After the separated S wave gathers are filtered by FK filtering to remove downflow waves and S wave static correction, the signal-to-noise ratio of the data is significantly improved. The pre-stack time migration submodule calls the GPU-accelerated Kirchhoff integral migration algorithm to migrate the P-wave gather, fast S-wave gather, and slow S-wave gather respectively. The migration aperture is set to 4000 meters, and the maximum tilt angle of the anti-spoofing filter is set to 70 degrees. The three migration profiles generated after migration are completely aligned in spatial coordinates, and can be directly used for layer comparison and joint analysis of multi-wave attributes.

[0039] Step S3 is performed by the wavefield attribute extraction module based on the migration profile. Interpreters use seismic interpretation software to trace the phase axes of the top and bottom interfaces of the target layer on the P-wave migration profile. The interpretation grid density is 50 m × 50 m, and the interpretation results are converted into a continuous spatial surface. The attribute extraction program samples uniformly within the layer segment using a 25 m × 25 m grid, calculating 15 attribute values ​​for each sampling point. Velocity attributes are extracted directly from the velocity model; reflection amplitude gradient attributes are obtained through AVO gradient fitting; attenuation coefficient attributes are calculated using the spectral ratio method; and dominant frequency and instantaneous phase attributes are calculated through short-time window spectral analysis and Hilbert transform. Approximately 3.8 million spatial sampling points are extracted across the entire work area, each corresponding to a 15-dimensional attribute vector. All attribute data is stored in an 800 × 1200 × 1 HDF5 dataset with 32-bit floating-point data, and the data file size is approximately 27 GB.

[0040] Step S5, model training, is conducted based on the attribute data volume generated in step S3 and the existing wellpoint stress data. Training sample attribute vectors are extracted from the HDF5 dataset according to wellpoint coordinates. An initial sample set is constructed by combining this vector with the measured stress label values ​​of the wellpoints, and then augmented to 1650 samples through data augmentation. The mean and standard deviation of 15 dimensions are calculated for all samples and saved as a standard parameter file. These statistics are then used to standardize the attribute data volume using Z-scores. The standardized training and validation sets are then fed into a deep learning inference server for model training. Training is performed using the PyTorch 1.12 framework, with automatic mixed-precision training enabled to accelerate computation and save GPU memory. The training runs on an A100 GPU with a batch size of 32 for 800 iterations, converging after 750 iterations. The resulting model weight parameter file, along with the standardized parameter file, is saved as a model deployment package.

[0041] Step S6 executes the final 3D geostress prediction. The standardized 15-dimensional attribute vector data from 3.8 million sampling points across the entire work area are sequentially input into the pre-trained deep learning model for deep geostress prediction. The model performs forward propagation calculations in batches, with a batch size of 256, resulting in approximately 14,844 batches across the entire work area. Each batch of data flows through a fixed pre-trained encoder convolutional layer and a bidirectional gated recurrent unit network, outputting the predicted minimum horizontal principal stress gradient, maximum horizontal principal stress gradient, and maximum horizontal principal stress azimuth at the fully connected regression layer. The entire prediction process takes approximately 23 minutes. The three prediction results are reconstructed into three data volumes using a spatial grid, with each data volume having the same size as the input attribute data volume. The final output and visualization module converts these three data volumes into standard SEG-Y format files and VTK structured mesh files, respectively, facilitating display and mapping in conventional seismic interpretation software and 3D visualization software. The prediction results show that the minimum horizontal principal stress gradient value of the target layer is between 16.8 kPa / m and 24.3 kPa / m, and the maximum horizontal principal stress gradient value is between 20.1 kPa / m and 29.5 kPa / m. The azimuth of the maximum horizontal principal stress gradually rotates from 45 degrees north of east in the core of the anticline to 110 degrees north of east in the flank. The difference in the direction of the maximum horizontal principal stress between the core and the flank is consistent with the background of the regional tectonic stress field. The average relative error between the prediction results and the measured values ​​of existing exploration wells is 8.3%, which is significantly reduced compared with the average relative error of 19.5% of the traditional method that only uses P-wave velocity and amplitude gradient information, thus verifying the technical effectiveness of the proposed scheme.

[0042] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for predicting deep geostress based on seismic wavefield information, characterized in that, Includes the following steps: S1: Deploy a multi-component seismic data acquisition and observation system on the surface or in wells of the target work area, and obtain raw multi-component common shot gather data by exciting and recording the full wavefield seismic record containing P-waves, S-waves and converted waves through artificial seismic sources. S2: The original multi-component common shot gather data is preprocessed and multi-domain wavefield separation is performed to obtain the time-domain common reflection gathers corresponding to the P-wave field, fast S-wave field and slow S-wave field respectively. Then, the separated wavefield data are subjected to pre-stack time migration processing to obtain the P-wave migration profile, fast S-wave migration profile and slow S-wave migration profile. S3: From the P-wave migration profile, fast S-wave migration profile, and slow S-wave migration profile, extract the P-wave layer velocity, fast S-wave layer velocity, slow S-wave layer velocity, P-wave reflection amplitude gradient, fast S-wave reflection amplitude gradient, slow S-wave reflection amplitude gradient, P-wave attenuation coefficient, fast S-wave attenuation coefficient, slow S-wave attenuation coefficient, P-wave dominant frequency, fast S-wave dominant frequency, slow S-wave dominant frequency, P-wave instantaneous phase, fast S-wave instantaneous phase, and slow S-wave instantaneous phase along the target layer to construct a multi-mode seismic wavefield oriented towards the target layer. The attribute data volume, wherein each spatial sampling point in the multimodal seismic wavefield attribute data volume corresponds to a multidimensional attribute vector, and the components of the attribute vector are, in order, P-wave layer velocity, fast S-wave layer velocity, slow S-wave layer velocity, P-wave reflection amplitude gradient, fast S-wave reflection amplitude gradient, slow S-wave reflection amplitude gradient, P-wave attenuation coefficient, fast S-wave attenuation coefficient, slow S-wave attenuation coefficient, P-wave dominant frequency, fast S-wave dominant frequency, slow S-wave dominant frequency, P-wave instantaneous phase, fast S-wave instantaneous phase, and slow S-wave instantaneous phase; S4: Construct a deep learning model for predicting deep geostress. The input layer of the deep learning model for predicting deep geostress receives the multidimensional attribute vector. Its hidden layer contains a first convolutional autoencoder, a second convolutional autoencoder, a bidirectional gated recurrent unit network, and a fully connected regression layer, which are cascaded in sequence. Its output layer outputs the minimum horizontal principal stress gradient value, the maximum horizontal principal stress gradient value, and the maximum horizontal principal stress azimuth angle. S5: Construct a training sample set using the measured geostress data at known well points and the multidimensional attribute vectors at corresponding spatial sampling points, and conduct supervised training on the deep learning model for deep geostress prediction. During training, a composite loss function containing adaptive adjustment terms for multi-task weights is used to constrain the update of model parameters until the composite loss function converges to a preset threshold, thereby obtaining the trained deep learning model for deep geostress prediction. S6: Input the multidimensional attribute vectors of all spatial sampling points in the multimodal seismic wavefield attribute data volume of the target layer to be predicted into the trained deep learning model for deep geostress prediction in sequence. After forward propagation calculation, the minimum horizontal principal stress gradient data volume, the maximum horizontal principal stress gradient data volume, and the maximum horizontal principal stress azimuth data volume of the entire target layer are obtained.

2. The method according to claim 1, characterized in that, The multi-component seismic data acquisition and observation system described in S1 adopts a linear or grid-like layout with equal or non-equal spacing. The single-channel detector is a three-component digital detector, whose three components record the particle vibration velocity or acceleration of the vertical component, the east-west horizontal component, and the north-south horizontal component, respectively.

3. The method according to claim 1, characterized in that, The preprocessing described in S2 includes sequentially performing observation system definition, trace editing, static correction, spherical diffusion compensation, surface uniformity amplitude compensation, deconvolution, and noise suppression on the original multi-component common shot point gather data. The multi-domain wavefield separation utilizes a combination of polarization filtering and rotational projection to project the original three-component seismic record from the Cartesian coordinate system composed of the vertical component, the east-west horizontal component, and the north-south horizontal component to the ray coordinate system. By analyzing the polarization ellipse parameters of the particle vibration trajectory, the P-wave field, fast S-wave field, and slow S-wave field are separated. The separated S-wave field is then subjected to up-and-downward wave separation and S-wave static correction processing.

4. The method according to claim 1, characterized in that, The P-wave velocity, fast S-wave velocity, and slow S-wave velocity mentioned in S3 are obtained through pre-stack time migration velocity analysis along the target layer. The P-wave reflection amplitude gradient, fast S-wave reflection amplitude gradient, and slow S-wave reflection amplitude gradient are obtained by fitting the amplitude variation curve of the reflected wave of the target layer with the migration distance in their respective migration profiles and calculating its slope. The P-wave attenuation coefficient, fast S-wave attenuation coefficient, and slow S-wave attenuation coefficient are obtained by calculating the dominant frequency migration of the top and bottom reflected waves of the target layer using the spectral ratio method and combining it with the interlayer travel time. The P-wave dominant frequency, fast S-wave dominant frequency, slow S-wave dominant frequency, P-wave instantaneous phase, fast S-wave instantaneous phase, and slow S-wave instantaneous phase are obtained by performing Fourier transform and Hilbert transform on the records within the time window of the target layer in their respective migration profiles.

5. The method according to claim 1, characterized in that, In S4, the first convolutional autoencoder has a preset first number of convolutional layers and deconvolutional layers, which are used to perform feature extraction and reconstruction pre-training on the input multidimensional attribute vector, and its encoder outputs a first intermediate feature vector; the second convolutional autoencoder has a preset second number of convolutional layers and deconvolutional layers, which receive the first intermediate feature vector and further extract deep abstract features, and its encoder outputs a second intermediate feature vector. The bidirectional gated recurrent unit network includes a forward propagation layer and a backward propagation layer, with a preset number of hidden units in each layer. The bidirectional gated recurrent unit network receives the second intermediate feature vector and captures the spatial correlation and mutual constraint relationship between multimodal attributes through bidirectional temporal modeling, outputting a temporal fusion feature vector. The fully connected regression layer includes a hidden layer with a preset first number of neurons and an output layer with a preset second number of neurons. It receives the temporal fusion feature vector and regresses to output the minimum horizontal principal stress gradient value, the maximum horizontal principal stress gradient value, and the maximum horizontal principal stress azimuth angle.

6. The method according to claim 1, characterized in that, The composite loss function described in S5 consists of a mean squared error loss term and a multi-task weight adaptive adjustment term. The mean squared error loss term is a weighted sum of the predicted mean squared error of the minimum horizontal principal stress gradient, the predicted mean squared error of the maximum horizontal principal stress gradient, and the predicted mean squared error of the maximum horizontal principal stress azimuth angle. The multi-task weight adaptive adjustment term dynamically adjusts the weight coefficients of different task losses based on the reciprocal of the loss decrease rate of each task on the validation set.

7. The method according to claim 1, characterized in that, The measured in-situ stress data at the known well points mentioned in S5 are obtained through any of the following methods: acoustic emission Kaiser effect test, used to obtain the memory information of the core under historical tectonic stress; differential strain analysis, to obtain the triaxial in-situ stress value by analyzing the strain recovery characteristics of the directional core after unloading; hydraulic fracturing method, to determine the closure pressure and re-tension pressure when the well wall ruptures through the fracturing curve, and to calculate the minimum horizontal principal stress. The training sample set is constructed as follows: for each known well, the well trajectory is spatially matched with the multimodal seismic wave field attribute data volume, and the multidimensional attribute vector of the target layer at the well point is extracted as the input sample. The measured minimum horizontal principal stress gradient value, the measured maximum horizontal principal stress gradient value, and the measured maximum horizontal principal stress azimuth angle of the corresponding depth segment are used as labels to form a training sample. All training samples corresponding to all known wells are summarized into the training sample set.

8. The method according to claim 1, characterized in that, Before being input into the deep learning model for deep geostress prediction, the multidimensional attribute vectors of each spatial sampling point in the multimodal seismic wavefield attribute data volume described in S3 are standardized using the statistical features of all attribute vectors in the training sample set, so that each dimension of the feature follows a preset statistical distribution. During training and prediction, the input data of the deep learning model for deep geostress prediction are the standardized multidimensional attribute vectors.

9. The method according to claim 1, characterized in that, The multidimensional attribute vector of each spatial sampling point in the multimodal seismic wavefield attribute data volume described in S3 is formed by simultaneously extracting the P-wave layer velocity, fast S-wave layer velocity, slow S-wave layer velocity, P-wave reflection amplitude gradient, fast S-wave reflection amplitude gradient, slow S-wave reflection amplitude gradient, P-wave attenuation coefficient, fast S-wave attenuation coefficient, slow S-wave attenuation coefficient, P-wave dominant frequency, fast S-wave dominant frequency, slow S-wave dominant frequency, P-wave instantaneous phase, fast S-wave instantaneous phase, and slow S-wave instantaneous phase from the P-wave migration profile, fast S-wave migration profile, and slow S-wave migration profile, and then organizing them according to the component order of the attribute vector.

10. A deep geostress prediction system based on seismic wavefield information, characterized in that, The system performs all the steps included in the deep geostress prediction method based on seismic wavefield information as described in any one of claims 1 to 9.

Images