Dual-phase medium elastic inversion system and method based on deep learning

Abstract

The present application relates to the field of deep learning, in particular to a dual-phase medium elastic inversion system and method based on deep learning. The seismic wave field data is received and time-frequency orthogonal transformation and adaptive band-pass filtering are performed to extract the characteristics of fast longitudinal wave, slow longitudinal wave and transverse wave field; the parameter regression network maps the above wave field characteristics to output solid phase skeleton elastic parameters and fluid phase fluid parameters; the output parameters are substituted into the dual-phase medium wave equation for forward calculation to obtain the theoretical synthetic wave field, and the waveform residual between the theoretical synthetic wave field and the seismic wave field is converted into a physical constraint loss term; the physical constraint loss term and the data-driven regression loss term are weighted and summed to obtain the total loss, and the network weight is updated based on the total loss. The present application overcomes the defect that the pure data-driven inversion result violates the physical law, reduces the multi-solution interference, improves the generalization ability and inversion stability.

Description

Technical Field

[0001] This invention relates to the field of deep learning, specifically to a deep learning-based two-phase medium elastic inversion system and method. Background Technology

[0002] In the field of geophysical exploration, elastic inversion of two-phase media aims to obtain the elastic parameters of the solid skeleton and the fluid parameters of the subsurface medium from seismic wavefield data. Existing conventional techniques typically employ purely data-driven deep neural networks to directly construct an end-to-end mapping between seismic wavefield observation data and the elastic parameters of the two-phase medium. In this approach, the network relies solely on pairs of wavefield data and parameter labels within the training set for forward propagation to calculate prediction errors, and updates network weights through backpropagation. The entire training process is controlled only by data-driven regression loss, lacking consideration of the physical mechanisms of interaction between the solid skeleton and pore fluids in the two-phase medium, thus constituting a black-box mapping mode.

[0003] The aforementioned pure data-driven black-box mapping mode results in a lack of constraints from the solid-fluid coupling physical mechanism between the solid skeleton parameters and the fluid parameters output by the deep learning network. This often leads to parameter combinations that violate the physical common sense of the wave equation in two-phase media, resulting in poor physical consistency of the predicted parameters. Summary of the Invention

[0004] The purpose of this invention is to provide a deep learning-based two-phase medium elastic inversion system and method, which can effectively solve the problems in the background art.

[0005] To achieve the above objectives, the technical solution adopted by the present invention is as follows: A deep learning-based two-phase medium elastic inversion system includes: a data input interface for receiving seismic wavefield data; The time-frequency decoupling feature extractor is configured to perform time-frequency orthogonal transformation and adaptive bandpass filtering on the seismic wavefield data to separate and extract fast P-wave wavefield features, slow P-wave wavefield features and shear wavefield features. The parameter regression network, coupled to the time-frequency decoupling feature extractor, receives the fast longitudinal wave field features, the slow longitudinal wave field features, and the transverse wave field features, and maps and outputs the solid-phase skeleton elastic parameters and the fluid parameters of the two-phase medium. A physical constraint loss builder, coupled with the parameter regression network, substitutes the solid skeleton elastic parameters and the fluid parameters into the Biot two-phase medium wave equation for finite difference forward modeling to obtain theoretical synthetic wave field data, calculates the waveform residual between the theoretical synthetic wave field data and the seismic wave field data, and transforms the waveform residual into a physical constraint loss term; A joint optimizer, coupled with the physical constraint loss builder and the parametric regression network, obtains the physical constraint loss term and the data-driven regression loss term based on labeled elastic parameters. The physical constraint loss term and the data-driven regression loss term are weighted and summed to obtain the total loss value. Gradient backpropagation is performed based on the total loss value to update the weight parameters of the time-frequency decoupled feature extractor and the parametric regression network.

[0006] Preferably, the time-frequency decoupling feature extractor includes a multi-scale time-frequency orthogonal decomposer and a bandgap energy focusing filter; The multi-scale time-frequency orthogonal decomposer performs multi-scale wavelet basis expansion on the seismic wavefield data along the time dimension to obtain a multi-scale time-frequency coefficient matrix. The frequency band energy focusing filter evaluates the energy concentration degree of the multi-scale time-frequency coefficient matrix along the frequency dimension, determines the main frequency band boundaries of fast P-wave, slow P-wave, and shear wave based on the energy concentration degree distribution, generates a frequency domain mask based on the main frequency band boundaries, applies the frequency domain mask to the multi-scale time-frequency coefficient matrix, separates the fast P-wave time-frequency coefficients, slow P-wave time-frequency coefficients, and shear wave time-frequency coefficients, and obtains the fast P-wave wave field characteristics, the slow P-wave wave field characteristics, and the shear wave field characteristics through inverse transformation, respectively.

[0007] Preferably, the parameter regression network includes a solid-fluid decoupling mapping layer and a physical boundary constraint output layer; The solid-fluid decoupling mapping layer receives the fast longitudinal wave field features, the slow longitudinal wave field features, and the transverse wave field features. It calculates the fluid influence weight matrix of the slow longitudinal wave field features on the fast longitudinal wave field features and the transverse wave field features through a cross-attention mechanism. It then fuses the fluid influence weight matrix with the fast longitudinal wave field features and the transverse wave field features to obtain solid-phase decoupling features and extracts the flow-phase decoupling features from the slow longitudinal wave field features. The physical boundary constraint output layer maps the solid phase decoupling features to the solid phase skeleton elastic parameter space and the fluid phase decoupling features to the fluid phase fluid parameter space. It also introduces a boundary threshold for the two-phase medium property parameters into the output activation function to truncate out-of-bounds predicted values.

[0008] Preferably, the physical constraint loss builder includes a Biot equation tensor quantization difference operator and a residual scaling aligner; The tensor difference operator of the Biot equation constructs the elastic parameters of the solid skeleton and the fluid parameters of the flow phase into a fourth-order tensor, which is then substituted into the partial differential operator of the Biot equation. The automatic differentiation mechanism is used to calculate the difference recursion relationship between the solid displacement field and the flow phase displacement field on the spatiotemporal grid to obtain the theoretical synthesized wave field data. The residual scaling aligner extracts the envelope difference between the theoretically synthesized wavefield data and the seismic wavefield data, performs amplitude normalization and phase alignment operations on the envelope difference, and uses the envelope difference after amplitude normalization and phase alignment operations as the waveform residual, thereby compressing the dynamic range of the waveform residual.

[0009] Preferably, the joint optimizer includes an adaptive loss weighter and a gradient conflict resolver; The adaptive loss weighter calculates a dynamic balance coefficient based on the variance of the physical constraint loss term and the variance of the data-driven regression loss term, and uses the dynamic balance coefficient to perform a weighted summation of the physical constraint loss term and the data-driven regression loss term to obtain the total loss value. The gradient conflict resolver extracts the physical gradient direction corresponding to the physical constraint loss term and the data gradient direction corresponding to the data-driven regression loss term during the gradient backpropagation process based on the total loss value. It calculates the cosine value of the angle between the physical gradient direction and the data gradient direction. When the cosine value is less than 0, it performs orthogonal projection decomposition of the physical gradient direction and the data gradient direction on the shared parameters.

[0010] Preferably, the physical constraint loss builder further includes a solid-fluid coupling energy conservation corrector; The solid-fluid coupling energy conservation corrector calculates the sum of solid-phase kinetic energy, fluid-phase kinetic energy, solid-phase potential energy, and fluid-phase potential energy based on the solid-phase skeleton elastic parameters, the fluid-phase parameters, and the theoretically synthesized wavefield data to obtain the predicted total energy. Based on the seismic wavefield data, it calculates the observed total energy, calculates the energy deviation integral between the predicted total energy and the observed total energy, multiplies the energy deviation integral by the energy conservation penalty coefficient, and then adds it to the waveform residual to obtain the corrected waveform residual. The corrected waveform residual is determined as the physical constraint loss term.

[0011] Preferably, the multi-scale time-frequency orthogonal decomposer includes an adaptive wavelet basis selector and a multi-resolution cascaded expander; The adaptive wavelet basis selector extracts the spectral centroid and kurtosis of the seismic wavefield data, and matches the main wavelet basis function and auxiliary wavelet basis function in a preset wavelet basis library based on the spectral centroid and the kurtosis. The multi-resolution cascaded expander uses the main wavelet basis function to perform coarse-scale low-frequency approximation decomposition on the seismic wavefield data, and uses the auxiliary wavelet basis function to perform fine-scale high-frequency detail decomposition on the low-frequency components obtained by the coarse-scale low-frequency approximation decomposition. The results of the coarse-scale low-frequency approximation decomposition and the fine-scale high-frequency detail decomposition are concatenated along the channel dimension to form the multi-scale time-frequency coefficient matrix, which is combined to form a cross-band asymmetric wavelet basis decomposition structure.

[0012] Preferably, the solid-fluid decoupling mapping layer includes a pore connectivity graph builder and a graph convolution aggregator; The pore connectivity graph builder extracts the fluid spatial distribution probability matrix based on the slow longitudinal wave field features, determines the pixel nodes in the fluid spatial distribution probability matrix that are greater than the probability threshold as graph nodes, calculates the edge weights between adjacent graph nodes based on wave field coherence, and constructs a non-Euclidean graph structure. The graph convolution aggregator performs message passing on the non-Euclidean graph structure, projects the fast longitudinal wave field features and the transverse wave field features onto the node feature space of the non-Euclidean graph structure, aggregates the skeleton elastic response features of neighboring nodes in the non-Euclidean graph structure, and outputs the solid-phase decoupling features that include pore topological associations.

[0013] Preferably, the residual scaling aligner includes a wavelet phase matcher and an amplitude envelope normalizer; The wavelet phase matcher performs Hilbert transform on the theoretical synthesized wavefield data and the seismic wavefield data respectively to obtain the instantaneous phase, calculates the phase difference between the instantaneous phase of the theoretical synthesized wavefield data and the instantaneous phase of the seismic wavefield data, and performs time shift correction on the theoretical synthesized wavefield data based on the phase difference to obtain the phase-aligned wavefield. The amplitude envelope normalizer extracts the amplitude envelopes of the phase-aligned wavefield and the seismic wavefield data, calculates the logarithmic ratio of the amplitude envelope of the phase-aligned wavefield to the amplitude envelope of the seismic wavefield data, scales the amplitude of the phase-aligned wavefield based on the logarithmic ratio, obtains the amplitude-normalized wavefield, and determines the difference between the amplitude-normalized wavefield and the seismic wavefield data as the envelope difference after amplitude normalization and phase alignment operations.

[0014] A deep learning-based method for elastic inversion of two-phase media includes: receiving seismic wavefield data; Perform time-frequency orthogonal transformation and adaptive bandpass filtering on the seismic wavefield data to separate and extract fast P-wave wavefield features, slow P-wave wavefield features and shear wavefield features; Receive the fast longitudinal wave field characteristics, the slow longitudinal wave field characteristics, and the transverse wave field characteristics, and map and output the solid-phase skeleton elastic parameters and the fluid parameters of the two-phase medium; Substitute the solid-phase skeleton elastic parameters and the fluid parameters into the Biot two-phase medium wave equation for finite difference forward modeling to obtain theoretical synthetic wave field data. Calculate the waveform residual between the theoretical synthetic wave field data and the seismic wave field data, and convert the waveform residual into a physical constraint loss term. Obtain the physical constraint loss term and the data-driven regression loss term based on the labeled elasticity parameters. Sum the physical constraint loss term and the data-driven regression loss term with weights to obtain the total loss value. Perform gradient backpropagation based on the total loss value to update the weight parameters for performing the time-frequency orthogonal transformation, adaptive bandpass filtering operation, and mapping operation.

[0015] Compared with the prior art, the beneficial effects of the present invention are as follows: 1. By introducing a physical constraint loss builder at the output of the parametric regression network, the solid-phase skeleton elastic parameters and fluid parameters output by the network are substituted into the two-phase medium wave equation for forward modeling. The waveform residuals of theoretical synthetic wavefield data and seismic wavefield data are transformed into physical constraint loss terms, which are jointly optimized with the data-driven regression loss terms to optimize the network weights. This enables the network to satisfy the physical conservation law of the two-phase medium wave equation while performing parametric regression, overcoming the defect that the pure data-driven inversion results violate the physical laws of solid-fluid coupling. It also reduces the multi-solution interference of multi-parameter joint inversion and enhances the model's generalization prediction ability for data not used in training.

[0016] 2. Independent wavefield features of fast P-waves, slow P-waves, and S-waves are extracted through time-frequency orthogonal transformation and adaptive bandpass filtering. A cross-attention mechanism is used to calculate the fluid influence weight matrix to obtain solid-phase decoupling and flow-phase decoupling features, matching the multi-waveform physical characteristics of the two-phase medium. A residual scaling aligner is used to normalize the amplitude and align the phase of the waveform residuals. Boundary thresholds for the two-phase medium property parameters are introduced into the output activation function to truncate out-of-bounds predicted values, improving the stability of parameter regression. A solid-fluid coupling energy conservation corrector is introduced to calculate the energy deviation integral superimposed on the waveform residuals, constraining the energy conservation of the predicted parameters. A dynamic equilibrium coefficient weighted loss term is calculated, and orthogonal projection decomposition is performed when the cosine of the angle between the physical gradient direction and the data gradient direction is less than 0, resolving the gradient conflict between physical constraints and data-driven approaches. Attached Figure Description

[0017] Figure 1 This is a flowchart illustrating the overall workflow of the deep learning-based two-phase medium elastic inversion system of the present invention. Figure 2 This is a flowchart of the time-frequency decoupling feature extractor of the present invention; Figure 3 This is a flowchart of the parameter regression network solid-fluid decoupling mapping and parameter output of the present invention; Figure 4 This is a flowchart of the waveform residual calculation process for the physical constraint loss builder of the present invention; Figure 5 This is a flowchart of the joint optimizer loss weighting and gradient conflict resolution process of the present invention; Figure 6This is a flowchart of the loss correction process for the solid-fluid coupling energy conservation corrector of the present invention. Detailed Implementation

[0018] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0019] Please refer to Figure 1 This embodiment provides a deep learning-based elastic inversion system and method for two-phase media. The deep learning-based elastic inversion system for two-phase media includes a data input interface, a time-frequency decoupling feature extractor, a parameter regression network, a physical constraint loss builder, and a joint optimizer. The data input interface receives seismic wavefield data, which is discrete sampling data on a two-dimensional or three-dimensional spatiotemporal grid. Each sampling point contains the particle vibration velocity or displacement amplitude at the corresponding spatiotemporal location. The time-frequency decoupling feature extractor performs time-frequency orthogonal transformation and adaptive bandpass filtering on the seismic wavefield data to separate and extract fast P-wave, slow P-wave, and S-wave features. The parameter regression network is coupled to the time-frequency decoupling feature extractor, receives the fast P-wave, slow P-wave, and S-wave features, and maps and outputs the solid-phase skeleton elastic parameters and the fluid parameters of the two-phase medium. A physical constraint loss builder is coupled to the parametric regression network. It substitutes the solid-phase skeleton elastic parameters and the fluid phase parameters into the Biot two-phase medium wave equation for finite-difference forward modeling to obtain theoretically synthesized wavefield data. It then calculates the waveform residual between the theoretically synthesized wavefield data and the seismic wavefield data, transforming this residual into a physical constraint loss term. A joint optimizer is coupled to the physical constraint loss builder and the parametric regression network. It obtains the physical constraint loss term and a data-driven regression loss term based on labeled elastic parameters. The physical constraint loss term and the data-driven regression loss term are weighted and summed to obtain a total loss value. Gradient backpropagation is then performed based on this total loss value to update the weight parameters of the time-frequency decoupled feature extractor and the parametric regression network.

[0020] Specifically, the seismic wavefield data received by the data input interface is represented as a tensor. ,in This indicates the number of spatial sampling points in the horizontal direction. This indicates the number of spatial sampling points in the vertical direction. Indicates the number of time sampling points. Each element Indicates spatial location ,time Seismic wavefield observations at the location. The time-frequency decoupling feature extractor... Performing a time-frequency orthogonal transform along the time dimension converts the time-domain signal into a time-frequency domain representation, yielding the time-frequency coefficient matrix. ,in This represents the number of frequency sampling points. The time-frequency orthogonal transform satisfies the orthogonality condition, meaning the product of the transform matrix and its transpose is the identity matrix, ensuring no information loss during the transformation process.

[0021] Time-frequency decoupling feature extractor for time-frequency coefficient matrix Adaptive bandpass filtering is performed to separate the time-frequency coefficients of the fast longitudinal wave based on the differences in energy distribution of different waveforms in the time-frequency domain. Slow longitudinal wave time-frequency coefficient With transverse wave time-frequency coefficient In the adaptive bandpass filtering process, the energy concentration degree of each frequency channel is first calculated. The energy concentration degree is defined as the ratio of the sum of squares of the time-frequency coefficients in that frequency channel to the total energy of all frequency channels. Based on the energy concentration degree distribution, the main frequency band ranges of fast P-waves, slow P-waves, and S-waves are determined, and corresponding frequency domain masks are generated. The frequency domain masks are then multiplied element-wise with the time-frequency coefficient matrix to obtain the time-frequency coefficients of each waveform.

[0022] For the time-frequency coefficients of fast longitudinal waves Slow longitudinal wave time-frequency coefficient With transverse wave time-frequency coefficient Perform inverse time-frequency transforms to obtain the characteristics of the fast P-wave field. Slow longitudinal wave field characteristics Characteristics of transverse wave field The inverse time-frequency transform is the inverse of the corresponding time-frequency orthogonal transform, satisfying the transform invertibility condition and ensuring that the time-domain wave field characteristics can be accurately recovered from the time-frequency coefficients.

[0023] Parameter regression network reception , and The solid-phase framework elastic parameters and the fluid parameters of the mobile phase in a two-phase medium are output through multilayer nonlinear mapping. The solid-phase framework elastic parameters include the solid bulk modulus. solid phase shear modulus With solid density The parameters of the fluid phase include the fluid bulk modulus. Fluid density With fluid viscosity The output of a parametric regression network is a parametric tensor. Each spatial location It corresponds to 6 elastic parameter values.

[0024] The physical constraint loss builder will convert the parameter tensor Substituting the values ​​into Biot's two-phase medium wave equation, finite-difference forward modeling is performed. Biot's two-phase medium wave equation describes the propagation of the solid-phase displacement field and the fluid-phase displacement field in a two-phase medium, considering the inertial, viscous, and elastic coupling between the solid skeleton and the porous fluid. Finite-difference forward modeling discretizes the continuous spatiotemporal domain into a uniform grid, and by approximating partial differential operations with finite differences, the recursive relationship between the solid-phase displacement field and the fluid-phase displacement field is obtained.

[0025] In finite difference forward modeling, the spatial domain is discretized into steps with a step size of The grid, discretized in the time domain with a step size of The grid. The solid-state displacement field is represented as... The flow phase displacement field is expressed as Using the second-order central difference approximation of spatial partial derivatives and the second-order central difference approximation of time partial derivatives, the solid-phase displacement field and the fluid-phase displacement field at the next time step are obtained: in Indicates the first Solid-state displacement field at each time step Indicates the first The flow phase displacement field at each time step , For the quality matrix, , , This is the stiffness matrix, whose elements are determined by the elastic parameters of the solid skeleton and the fluid parameters of the fluid phase.

[0026] The solid-phase displacement field and the fluid-phase displacement field at all time steps are obtained through recursive calculation. The solid-phase displacement field components at the surface observation locations are extracted to obtain the theoretically synthesized wavefield data. Computational theoretical synthesized wavefield data Input seismic wavefield data waveform residuals between : in This is the surface observation component of the input seismic wavefield data. The waveform residuals are then... The L2 norm is used as a physical constraint loss term. : The joint optimizer obtains the physical constraint loss term. With data-driven regression loss term The data-driven regression loss term is calculated based on labeled elasticity parameters, when the training data includes labeled biphasic medium elasticity parameters. At that time, the data-driven regression loss term is the output of the parametric regression network. With annotation parameters L2 norm between: The total loss value is obtained by weighting and summing the physical constraint loss term and the data-driven regression loss term. : in Data-driven loss weights The physical constraint loss weights are based on the total loss value. Gradient backpropagation is performed to calculate the gradient of the loss function with respect to all weight parameters of the time-frequency decoupled feature extractor and the parametric regression network. The weight parameters are then updated using the gradient descent algorithm. The gradient descent algorithm employs the Adam optimizer, and its parameter update rule is as follows: in For the first The gradient of the next iteration. This is a first-moment estimate of the gradient. For the second moment estimation of the gradient, , The exponential decay rate is estimated by moments. For learning rate, For numerically stable terms, For the first The weight parameters for the next iteration.

[0027] Table 1 shows the physical meaning and typical value range of the elastic parameters of a two-phase medium. These parameters together determine the propagation characteristics of seismic waves in a two-phase medium.

[0028] Table 1. Physical meaning and typical value range of elastic parameters of two-phase media. In this embodiment, the elastic parameters output by the parametric regression network are substituted into the Biot two-phase medium wave equation for forward modeling to construct a physical constraint loss term. This physical constraint loss term is then jointly optimized with the data-driven regression loss term to optimize the network weights, ensuring that the network simultaneously satisfies data fitting and physical constraints during the parametric regression process. The physical constraint loss term integrates the physical mechanism of the two-phase medium wave equation into the network training process, ensuring that the combination of elastic parameters output by the network conforms to the physical laws of solid-fluid coupling and avoiding prediction results that violate common physical sense.

[0029] In a preferred embodiment, reference Figure 2The time-frequency decoupling feature extractor includes a multi-scale time-frequency orthogonal decomposer and a bandgap energy focusing filter. The multi-scale time-frequency orthogonal decomposer performs multi-scale wavelet basis expansion on the seismic wavefield data along the time dimension to obtain a multi-scale time-frequency coefficient matrix. The bandgap energy focusing filter evaluates the energy concentration degree of the multi-scale time-frequency coefficient matrix along the frequency dimension. Based on the energy concentration degree distribution, it determines the main frequency band boundaries of fast P-waves, slow P-waves, and S-waves. A frequency domain mask is generated based on the main frequency band boundaries. The frequency domain mask is applied to the multi-scale time-frequency coefficient matrix to separate the time-frequency coefficients of fast P-waves, slow P-waves, and S-waves. The characteristics of the fast P-wave wavefield, slow P-wave wavefield, and S-wave wavefield are then obtained by inverse transformation.

[0030] Specifically, the multi-scale time-frequency orthogonal decomposer employs wavelet transform as the time-frequency orthogonal transform method. Wavelet transform has excellent time-frequency localization characteristics, enabling simultaneous analysis of signals in both the time and frequency domains. Multi-scale wavelet basis expansion decomposes seismic wavefield data into wavelet coefficients at different scales, with each scale corresponding to a different frequency range. The coarse scale corresponds to low-frequency components, and the fine scale corresponds to high-frequency components.

[0031] The multi-scale time-frequency orthogonal decomposer comprises an adaptive wavelet basis selector and a multi-resolution cascaded expander. The adaptive wavelet basis selector extracts the spectral centroid and kurtosis of the seismic wavefield data, and matches the main wavelet basis functions and auxiliary wavelet basis functions from a pre-defined wavelet basis library based on the spectral centroid and kurtosis. (Spectral centroid...) Defined as the weighted average frequency of the signal spectrum: in Fourier transform of seismic wavefield data, The Nyquist frequency. Spectral kurtosis. Defined as the ratio of the fourth central moment of the signal spectrum to the square of the second central moment: The centroid of the spectrum reflects the main concentrated frequency of the signal energy, while the kurtosis of the spectrum reflects the sharpness of the signal spectrum. The adaptive wavelet basis selector determines the center frequency of the main wavelet basis function based on the centroid and the support length of the main wavelet basis function based on the kurtosis. When the kurtosis is large, a wavelet basis function with a shorter support length is selected to obtain better frequency resolution; when the kurtosis is small, a wavelet basis function with a longer support length is selected to obtain better time resolution.

[0032] The pre-defined wavelet basis library includes various orthogonal wavelet basis functions such as Daubechies wavelet, Symlet wavelet, Coiflet wavelet, and Morlet wavelet. The adaptive wavelet basis selector calculates the cross-correlation coefficient between the seismic wavefield data and each wavelet basis function, selects the wavelet basis function with the largest cross-correlation coefficient as the primary wavelet basis function, and selects the wavelet basis function with the second largest cross-correlation coefficient as the auxiliary wavelet basis function.

[0033] The multi-resolution cascaded expander performs coarse-scale low-frequency approximation decomposition on the seismic wavefield data using the main wavelet basis function, and performs fine-scale high-frequency detail decomposition on the low-frequency components obtained from the coarse-scale low-frequency approximation decomposition using auxiliary wavelet basis functions. The coarse-scale low-frequency approximation decomposition employs the Mallat algorithm of discrete wavelet transform to decompose the signal into low-frequency approximation components and high-frequency detail components. For the... Low-frequency approximation components of the scale The decomposition process is as follows: in These are the coefficients of the low-pass filter. These are the coefficients of the high-pass filter, satisfying the conditions for an orthogonal mirror filter. For the first The low-frequency approximation component of the scale, For the first High-frequency detail components at scale.

[0034] Multi-resolution cascaded expander execution Level coarse-scale decomposition yields High-frequency detail components With one low-frequency approximation component For low-frequency approximation components Execution using auxiliary wavelet basis functions Fine-scale decomposition yields Fine-scale high-frequency detail components The high-frequency detail components obtained from coarse-scale decomposition are concatenated with those obtained from fine-scale decomposition along the channel dimension to form a multi-scale time-frequency coefficient matrix. .

[0035] The bandgap energy focusing filter evaluates the energy focusing degree of the multi-scale time-frequency coefficient matrix along the frequency dimension. For the first... The frequency channel (corresponding to the first frequency channel) Wavelet coefficients at various scales), energy concentration Defined as the ratio of the sum of squares of the time-frequency coefficients within that channel to the total energy of all channels: in This represents the total number of frequency channels. The energy concentration of all frequency channels is calculated to obtain the energy concentration distribution curve. Peak detection is performed on the energy concentration distribution curve to identify three main peaks, corresponding to the main frequency band centers of the fast P-wave, slow P-wave, and S-wave, respectively.

[0036] Using the center of each main frequency band as a reference, extend outwards to both sides until the energy concentration drops to 10% of the peak value to determine the main frequency band boundary of each waveform. The main frequency band boundary of the fast P-wave is... The dominant frequency band boundary of the slow longitudinal wave is The boundary of the main frequency band of the transverse wave is Generate a frequency domain mask based on the main frequency band boundary. , and The mask has 1s for the elements in the corresponding main frequency band and 0s for the rest.

[0037] Multiplying the frequency domain mask element-wise with the multi-scale time-frequency coefficient matrix yields the time-frequency coefficients for each waveform: Perform inverse wavelet transform on the time-frequency coefficients of each wave type to obtain the characteristics of the fast P-wave field. Slow longitudinal wave field characteristics Characteristics of transverse wave field The inverse wavelet transform is the inverse process of the wavelet transform. It uses the reconstruction filter coefficients to reconstruct wavelet coefficients at different scales into time-domain signals. Table 2 shows the characteristics and applicable scenarios of different types of wavelet basis functions. The adaptive wavelet basis selector selects the appropriate wavelet basis function according to the spectral characteristics of the seismic wavefield data to obtain the optimal time-frequency decomposition effect.

[0038] Table 2 Comparison of Wavelet Basis Function Characteristics and Applicable Scenarios refer to Figure 3The parameter regression network comprises a solid-fluid decoupling mapping layer and a physical boundary constraint output layer. The solid-fluid decoupling mapping layer receives fast P-wave, slow P-wave, and S-wave features. It calculates the fluid influence weight matrix of the slow P-wave features on the fast P-wave and S-wave features using a cross-attention mechanism. This fluid influence weight matrix is ​​then fused with the fast P-wave and S-wave features to obtain solid-phase decoupling features, and the slow P-wave features are extracted to obtain fluid-phase decoupling features. The physical boundary constraint output layer maps the solid-phase decoupling features to the solid-phase skeleton elastic parameter space and the fluid-phase decoupling features to the fluid-phase parameter space. A boundary threshold for the two-phase medium properties is introduced into the output activation function to truncate out-of-bounds predicted values.

[0039] Specifically, the solid-fluid decoupling mapping layer first maps the fast longitudinal wave field characteristics Slow longitudinal wave field characteristics Characteristics of transverse wave field Mapped to feature tensors of the same dimension through convolutional layers , and All dimensions ,in The number of feature channels. The cross-attention mechanism uses slow longitudinal wavefield features. As a query, the fast longitudinal wave field characteristics Characteristics of transverse wave field The fluid influence weight matrix is ​​calculated using the key and value as keys.

[0040] For the characteristics of fast longitudinal wavefields, the fluid influence weighting matrix The calculation process is as follows: in The bond tensor is obtained by linear transformation of the characteristics of the fast longitudinal wave field. This represents the dot product of two eigenvectors. This is a scaling factor used to prevent the gradient of the softmax function from vanishing due to excessively large dot product values.

[0041] Similarly, for the characteristics of the transverse wave field, the fluid influence weighting matrix The calculation process is as follows: in The bond tensor is obtained by linear transformation of the transverse wave field characteristics.

[0042] Fluid Influence Weight Matrix and This characterizes the degree to which the fluid characteristics at each spatial location influence the skeletal characteristics at other spatial locations. Multiplying the fluid influence weight matrix by the corresponding value tensor yields the fluid influence features: in and These are the value tensors obtained by linear transformation of the characteristics of the fast longitudinal wave field and the transverse wave field, respectively.

[0043] Subtracting the corresponding fluid influence characteristics from the fast longitudinal wave field characteristics and the transverse wave field characteristics yields the solid-phase decoupling characteristics: Will and By splicing along the channel dimension, the final solid-state decoupling characteristics are obtained. Characteristics of slow longitudinal wave field Feature extraction is performed using convolutional layers to obtain flow phase decoupling features. .

[0044] In a preferred embodiment, the solid-fluid decoupling mapping layer includes a pore connectivity graph builder and a graph convolution aggregator. The pore connectivity graph builder extracts the fluid spatial distribution probability matrix based on slow P-wave wavefield features, identifies pixel nodes in the fluid spatial distribution probability matrix that are greater than a probability threshold as graph nodes, calculates edge weights between adjacent graph nodes based on wavefield coherence, and constructs a non-Euclidean graph structure. The graph convolution aggregator performs message passing on the non-Euclidean graph structure, projects fast P-wave and S-wave wavefield features onto the node feature space of the non-Euclidean graph structure, aggregates the skeleton elastic response features of neighboring nodes in the non-Euclidean graph structure, and outputs solid-phase decoupling features containing pore topological associations.

[0045] Specifically, the pore connectivity graph constructor first analyzes the slow longitudinal wave field characteristics. Performing convolution operations and sigmoid activation yields the probability matrix of the fluid spatial distribution. Each element Indicates spatial location The probability of pore fluid existing at a given location. The probability matrix containing values ​​greater than a probability threshold. The pixel nodes are determined as graph nodes, and the set of graph nodes is represented as follows: .

[0046] For any two adjacent graph nodes and Calculate the edge weights between them. The edge weights are calculated based on the coherence of the slow longitudinal wave fields at the two nodes. Coherence is defined as the cross-correlation coefficient between the two wave field signals. edge weight This reflects the connectivity of the pore fluid between two nodes; higher coherence indicates better connectivity. Constructing a non-Euclidean graph structure. ,in Let the set be an edge set, and each edge has a corresponding edge weight. .

[0047] The graph convolution aggregator extracts fast longitudinal wave field features. Characteristics of transverse wave field Projecting onto the feature space of the graph nodes yields the initial feature vector for each graph node. Perform on a graph structure Layered graph convolution operations aggregate feature information from neighboring nodes at each layer: in For nodes The set of neighboring nodes, For nodes The degree, and For the first Learnable parameters of layer graph convolution. This is the activation function.

[0048] go through After layered graph convolutional aggregation, the features of each graph node contain skeletal elastic response information related to the pore topology in its neighborhood. Mapping the graph node features back to the original mesh space yields solid-state decoupling features containing pore topology associations. .

[0049] The physical boundary constraint output layer comprises two parallel convolutional neural network branches, processing solid-phase decoupling features and flow-phase decoupling features respectively. The solid-phase branch processes the solid-phase decoupling features... Mapped to the elastic parameter space of the solid skeleton, the solid bulk modulus is output. solid phase shear modulus With solid density The flow phase branching decouples the flow phases. Mapped to the fluid phase parameter space, output fluid bulk modulus. Fluid density With fluid viscosity .

[0050] At the output layer of each branch, the output parameters are truncated using an activation function with boundary constraints. The activation function with boundary constraints is defined as follows: in The sigmoid function maps the input to... interval, and These represent the physical boundary minimum and maximum values ​​of the corresponding parameters. The activation function with boundary constraints ensures that the elasticity parameters output by the network always remain within a physically reasonable range, avoiding out-of-bounds predictions.

[0051] In this embodiment, wavefield features of different wave types are separated through multi-scale time-frequency orthogonal decomposition and adaptive bandpass filtering. A cross-attention mechanism or graph convolution aggregator is used to extract solid-fluid decoupling features, and a physical boundary constraint output layer ensures the physical rationality of parameter predictions. Multi-scale time-frequency decomposition can capture the characteristics of seismic wavefields in different frequency ranges, and adaptive bandpass filtering can automatically determine the dominant frequency band of each wave type based on data characteristics, improving the accuracy of wavefield separation. The solid-fluid decoupling mapping layer can separate fluid influences from the skeleton features, improving the degree of decoupling between solid-phase parameter and fluid-phase parameter predictions.

[0052] In a preferred embodiment, reference Figure 4 The physical constraint loss builder includes a Biot equation tensor-based difference operator and a residual scaling aligner. The Biot equation tensor-based difference operator constructs fourth-order tensors from the solid-state elastic parameters and the fluid parameters of the flow phase. These tensors are then substituted into the partial differential operators of the Biot equation, and an automatic differentiation mechanism is used to calculate the difference recursion between the solid-state displacement field and the flow-state displacement field on the spatiotemporal grid, obtaining the theoretically synthesized wavefield data. The residual scaling aligner extracts the envelope difference between the theoretically synthesized wavefield data and the seismic wavefield data, performs amplitude normalization and phase alignment operations on the envelope difference, and uses the envelope difference after amplitude normalization and phase alignment as the waveform residual, compressing the dynamic range of the waveform residual.

[0053] Specifically, the tensor difference operator of the Biot equation constructs fourth-order tensors from the elastic parameters of the solid skeleton and the fluid parameters of the flow phase. The first three channels correspond to solid-phase parameters, and the last three channels correspond to fluid-phase parameters. Substituting the parameter tensors into the partial differential operators of Biot's two-phase medium wave equation, the tensor form of Biot's equation is: in The solid displacement vector, The flow phase displacement vector, To fix the corresponding force tensor, Pore ​​fluid pressure, The total density of the medium, Porosity For the effective mass of the flow phase, The viscous coupling coefficient is... This refers to penetration rate.

[0054] The relationship between the stress tensor and the displacement field is as follows: The relationship between pore fluid pressure and displacement field is as follows: in For strain tensor, For Biot coefficient, The effective bulk modulus of the fluid. This represents the bulk modulus of dry rock.

[0055] The tensor-based difference operator for the Biot equations utilizes an automatic differentiation mechanism to compute the difference approximation of the partial differential operators. This automatic differentiation mechanism constructs a computational graph, decomposing complex partial differential operations into a combination of basic operations, and can automatically compute derivatives of arbitrary order. In finite-difference forward modeling, the displacement field on the spatiotemporal grid is used as the nodes of the computational graph, and the parameter tensor is used as the input to the graph. The automatic differentiation mechanism then computes the partial derivatives of the displacement field with respect to time and space.

[0056] The partial derivatives calculated using the automatic differentiation mechanism are substituted into the difference recursion relation of the Biot equation to obtain the solid-phase displacement field and the fluid-phase displacement field at the next time step. The displacement fields at all time steps are recursively calculated, and the displacement components at the surface observation locations are extracted to obtain the theoretically synthesized wavefield data. .

[0057] The residual scaling aligner comprises a wavelet phase matcher and an amplitude envelope normalizer. The wavelet phase matcher performs Hilbert transforms on both the theoretical synthesized wavefield data and the seismic wavefield data to obtain the instantaneous phase. It calculates the phase difference between the instantaneous phases of the theoretical synthesized wavefield data and the seismic wavefield data, and performs time-shift correction on the theoretical synthesized wavefield data based on this phase difference to obtain the phase-aligned wavefield. The amplitude envelope normalizer extracts the amplitude envelopes of the phase-aligned wavefield and the seismic wavefield data, calculates the logarithmic ratio of the amplitude envelopes of the phase-aligned wavefield and the seismic wavefield data, scales the amplitude of the phase-aligned wavefield based on this logarithmic ratio, obtains the amplitude-normalized wavefield, and determines the difference between the amplitude-normalized wavefield and the seismic wavefield data as the envelope difference after amplitude normalization and phase alignment operations.

[0058] Specifically, the Hilbert transform is defined as: For real-valued signals Its analytical signal is: in The amplitude envelope of the signal. This represents the instantaneous phase of the signal.

[0059] Wavelet phase matcher for theoretically synthesized wavefield data With seismic wavefield data Perform Hilbert transforms on each to obtain their instantaneous phases. and Calculate each spatial location. Average phase difference at : Calculate time shift based on average phase difference : in For seismic wavefield data in spatial location The dominant frequency at that location. Time shift correction is performed on the theoretically synthesized wavefield data based on the time shift, resulting in a phase-aligned wavefield. : Amplitude envelope normalizer extracts phase-aligned wave field With seismic wavefield data amplitude envelope and Calculate each spatial location. Amplitude scaling factor at the location : The phase-aligned wavefield is scaled using an amplitude scaling factor to obtain the amplitude-normalized wavefield. : The difference between the amplitude-normalized wavefield and the seismic wavefield data is used as the waveform residual. : In a preferred embodiment, reference Figure 6The physical constraint loss builder also includes a solid-fluid coupling energy conservation corrector. Based on the solid skeleton elastic parameters, fluid parameters, and theoretically synthesized wavefield data, the solid-fluid coupling energy conservation corrector calculates the sum of solid kinetic energy, fluid kinetic energy, solid potential energy, and fluid potential energy to obtain the predicted total energy. Based on the seismic wavefield data, it calculates the observed total energy, calculates the energy deviation integral between the predicted total energy and the observed total energy, multiplies the energy deviation integral by the energy conservation penalty coefficient, and then superimposes it onto the waveform residual to obtain the corrected waveform residual. This corrected waveform residual is then determined as the physical constraint loss term.

[0060] Specifically, the solid-fluid coupled energy conservation corrector calculates the energy density at each spatiotemporal grid point. (Solid phase kinetic energy density) for: kinetic energy density of the flow phase for: Solid potential energy density for: Flow phase potential energy density for: in For fluid volume strain.

[0061] Total energy density The sum of all energy densities: Integrating the total energy density over the entire spatial domain yields the predicted total energy. : in and These represent the horizontal and vertical spatial lengths, respectively. The total observed energy is calculated based on seismic wavefield data. : in The average density of the medium.

[0062] Calculate the energy deviation integral between the predicted total energy and the observed total energy. : in The total recording time. The energy deviation integral is multiplied by the energy conservation penalty coefficient. We obtain the energy conservation correction term. : The energy conservation correction term is superimposed onto the L2 norm of the waveform residual to obtain the corrected physical constraint loss term: refer to Figure 5 The joint optimizer comprises an adaptive loss weighter and a gradient conflict resolver. The adaptive loss weighter calculates dynamic equilibrium coefficients based on the variances of the physical constraint loss term and the data-driven regression loss term. These coefficients are then used to weight and sum the physical constraint loss term and the data-driven regression loss term to obtain the total loss value. During gradient backpropagation based on the total loss value, the gradient conflict resolver extracts the physical gradient direction corresponding to the physical constraint loss term and the data gradient direction corresponding to the data-driven regression loss term. It calculates the cosine of the angle between the physical and data gradient directions. When the cosine is less than 0, the physical and data gradient directions are orthogonally projected onto shared parameters.

[0063] Specifically, the adaptive loss weighter maintains two sliding windows, each recording the most recent... The physical constraint loss value and the data-driven regression loss value are calculated for each iteration. The variance of the physical constraint loss value is also calculated. variance of data-driven regression loss value : in and These are the moving averages of the physical constraint loss value and the data-driven regression loss value, respectively.

[0064] Dynamic balance coefficient and Calculation based on the reciprocal of variance: This weighting method ensures that the two loss terms contribute approximately equally to the total loss, preventing one loss term from dominating the training process. The total loss value is: During gradient backpropagation, the gradient conflict resolver calculates the gradient of the physical constraint loss term with respect to the shared parameters. Gradient of the loss term with respect to shared parameters in data-driven regression Calculate the cosine of the angle between the two gradient directions. : when When the two gradients are in opposite directions, it indicates a gradient conflict. In this case, orthogonal projection decomposition of the two gradients is performed to eliminate the conflicting components. The correction value for the physical gradient is: The correction value for the data gradient is: The corrected gradients are summed to obtain the final gradient update direction: Table 3 shows the characteristics and effects of different gradient conflict resolution strategies. The orthogonal projection decomposition strategy can effectively eliminate gradient conflicts between physical constraints and data-driven approaches, thereby improving training stability.

[0065] Table 3 Comparison of Gradient Conflict Resolution Strategies In this embodiment, efficient finite-difference forward modeling is achieved through the Biot equation tensor quantization difference operator. A residual scaling aligner is used to normalize the amplitude and align the phase of the waveform residuals. A solid-fluid coupling energy conservation corrector is introduced to constrain the energy conservation of the predicted parameters. Joint optimization of physical constraints and data-driven approaches is achieved through an adaptive loss weighter and a gradient conflict resolver. The tensor quantization difference operator simplifies the solution process of partial differential equations using an automatic differentiation mechanism, improving the efficiency and differentiability of the forward modeling. The residual scaling aligner eliminates amplitude and phase differences between the theoretically synthesized wavefield and the observed wavefield, improving the effectiveness of the physical constraint loss. The energy conservation corrector further enhances the strength of the physical constraints, ensuring that the predicted parameters satisfy the energy conservation law. Adaptive loss weighting and gradient conflict resolution balance the contributions of physical constraints and data-driven approaches, avoiding gradient vanishing or exploding problems during training.

[0066] The deep learning-based elastic inversion method for two-phase media includes the following steps: receiving seismic wavefield data; performing time-frequency orthogonal transformation and adaptive bandpass filtering on the seismic wavefield data to separate and extract fast P-wave, slow P-wave, and S-wave features; receiving the fast P-wave, slow P-wave, and S-wave features, and mapping them to output the solid-phase skeleton elastic parameters and fluid parameters of the two-phase medium; substituting the solid-phase skeleton elastic parameters and fluid parameters into the Biot two-phase medium wave equation for finite-difference forward modeling to obtain theoretically synthesized wavefield data; calculating the waveform residual between the theoretically synthesized wavefield data and the seismic wavefield data, and converting the waveform residual into a physical constraint loss term; obtaining the physical constraint loss term and the data-driven regression loss term based on the labeled elastic parameters; weighted summing the physical constraint loss term and the data-driven regression loss term to obtain the total loss value; and performing gradient backpropagation based on the total loss value to update the weight parameters of the time-frequency orthogonal transformation, adaptive bandpass filtering, and mapping operations.

[0067] In a preferred embodiment, the steps of performing time-frequency orthogonal transformation and adaptive bandpass filtering on seismic wavefield data include: performing multi-scale wavelet basis expansion on the seismic wavefield data along the time dimension to obtain a multi-scale time-frequency coefficient matrix; evaluating the energy concentration of the multi-scale time-frequency coefficient matrix along the frequency dimension, and determining the main frequency band boundaries of fast P-wave, slow P-wave, and S-wave based on the energy concentration distribution; generating a frequency domain mask based on the main frequency band boundaries, applying the frequency domain mask to the multi-scale time-frequency coefficient matrix to separate the fast P-wave time-frequency coefficients, slow P-wave time-frequency coefficients, and S-wave time-frequency coefficients; and performing inverse transformations on the fast P-wave time-frequency coefficients, slow P-wave time-frequency coefficients, and S-wave time-frequency coefficients respectively to obtain the fast P-wave wavefield characteristics, slow P-wave wavefield characteristics, and S-wave wavefield characteristics.

[0068] In a preferred embodiment, the step of mapping the solid-phase skeleton elastic parameters and the fluid parameters of the output two-phase medium includes: calculating the fluid influence weight matrix of the slow longitudinal wave field characteristics on the fast longitudinal wave field characteristics and the transverse wave field characteristics through a cross-attention mechanism; fusing the fluid influence weight matrix with the fast longitudinal wave field characteristics and the transverse wave field characteristics to obtain solid-phase decoupling features, and extracting the fluid-phase decoupling features from the slow longitudinal wave field characteristics; mapping the solid-phase decoupling features to the solid-phase skeleton elastic parameter space, and mapping the fluid-phase decoupling features to the fluid parameter space; and introducing boundary thresholds for the two-phase medium physical property parameters into the output activation function to truncate out-of-bounds predicted values.

[0069] In a preferred embodiment, the step of substituting the solid-phase skeleton elastic parameters and the fluid phase parameters into the Biot two-phase medium wave equation for finite difference forward modeling includes: constructing the solid-phase skeleton elastic parameters and the fluid phase parameters as fourth-order tensors; substituting them into the partial differential operator of the Biot equation, and using an automatic differentiation mechanism to calculate the difference recursive relationship between the solid-phase displacement field and the fluid phase displacement field on the spatiotemporal grid; and obtaining theoretically synthesized wave field data.

[0070] In a preferred embodiment, the step of calculating the waveform residual between the theoretically synthesized wavefield data and the seismic wavefield data includes: performing Hilbert transform on the theoretically synthesized wavefield data and the seismic wavefield data respectively to obtain the instantaneous phase; calculating the phase difference between the instantaneous phase of the theoretically synthesized wavefield data and the instantaneous phase of the seismic wavefield data; performing time-shift correction on the theoretically synthesized wavefield data based on the phase difference to obtain the phase-aligned wavefield; extracting the amplitude envelopes of the phase-aligned wavefield and the seismic wavefield data; calculating the logarithmic ratio of the amplitude envelope of the phase-aligned wavefield to the amplitude envelope of the seismic wavefield data; scaling the amplitude of the phase-aligned wavefield based on the logarithmic ratio to obtain the amplitude-normalized wavefield; and determining the difference between the amplitude-normalized wavefield and the seismic wavefield data as the waveform residual.

[0071] In a preferred embodiment, the step of obtaining the total loss value by weighted summation of the physical constraint loss term and the data-driven regression loss term includes: calculating a dynamic balance coefficient based on the variance of the physical constraint loss term and the variance of the data-driven regression loss term; and using the dynamic balance coefficient to perform a weighted summation of the physical constraint loss term and the data-driven regression loss term to obtain the total loss value.

[0072] In a preferred embodiment, the step of performing gradient backpropagation based on the total loss value includes: extracting the physical gradient direction corresponding to the physical constraint loss term and the data gradient direction corresponding to the data-driven regression loss term respectively; calculating the cosine value of the angle between the physical gradient direction and the data gradient direction; when the cosine value of the angle is less than 0, performing orthogonal projection decomposition of the physical gradient direction and the data gradient direction on the shared parameters; and updating the weight parameters based on the decomposed gradient.

Claims

1. A deep learning-based two-phase medium elastic inversion system, characterized in that, include: Data input interface, used to receive seismic wavefield data; The time-frequency decoupling feature extractor is configured to perform time-frequency orthogonal transformation and adaptive bandpass filtering on the seismic wavefield data to separate and extract fast P-wave wavefield features, slow P-wave wavefield features and shear wavefield features. The parameter regression network, coupled to the time-frequency decoupling feature extractor, receives the fast longitudinal wave field features, the slow longitudinal wave field features, and the transverse wave field features, and maps and outputs the solid-phase skeleton elastic parameters and the fluid parameters of the two-phase medium. A physical constraint loss builder, coupled with the parameter regression network, substitutes the solid skeleton elastic parameters and the fluid parameters into the Biot two-phase medium wave equation for finite difference forward modeling to obtain theoretical synthetic wave field data, calculates the waveform residual between the theoretical synthetic wave field data and the seismic wave field data, and transforms the waveform residual into a physical constraint loss term; A joint optimizer, coupled with the physical constraint loss builder and the parametric regression network, obtains the physical constraint loss term and the data-driven regression loss term based on labeled elastic parameters. The physical constraint loss term and the data-driven regression loss term are weighted and summed to obtain the total loss value. Gradient backpropagation is performed based on the total loss value to update the weight parameters of the time-frequency decoupled feature extractor and the parametric regression network.

2. The system according to claim 1, characterized in that, The time-frequency decoupling feature extractor includes a multi-scale time-frequency orthogonal decomposer and a frequency band energy focusing filter; The multi-scale time-frequency orthogonal decomposer performs multi-scale wavelet basis expansion on the seismic wavefield data along the time dimension to obtain a multi-scale time-frequency coefficient matrix. The frequency band energy focusing filter evaluates the energy concentration degree of the multi-scale time-frequency coefficient matrix along the frequency dimension, determines the main frequency band boundaries of fast P-wave, slow P-wave, and shear wave based on the energy concentration degree distribution, generates a frequency domain mask based on the main frequency band boundaries, applies the frequency domain mask to the multi-scale time-frequency coefficient matrix, separates the fast P-wave time-frequency coefficients, slow P-wave time-frequency coefficients, and shear wave time-frequency coefficients, and obtains the fast P-wave wave field characteristics, the slow P-wave wave field characteristics, and the shear wave field characteristics through inverse transformation, respectively.

3. The system according to claim 2, characterized in that, The parameter regression network includes a solid-fluid decoupling mapping layer and a physical boundary constraint output layer; The solid-fluid decoupling mapping layer receives the fast longitudinal wave field features, the slow longitudinal wave field features, and the transverse wave field features. It calculates the fluid influence weight matrix of the slow longitudinal wave field features on the fast longitudinal wave field features and the transverse wave field features through a cross-attention mechanism. It then fuses the fluid influence weight matrix with the fast longitudinal wave field features and the transverse wave field features to obtain solid-phase decoupling features and extracts the flow-phase decoupling features from the slow longitudinal wave field features. The physical boundary constraint output layer maps the solid phase decoupling features to the solid phase skeleton elastic parameter space and the fluid phase decoupling features to the fluid phase fluid parameter space. It also introduces a boundary threshold for the two-phase medium property parameters into the output activation function to truncate out-of-bounds predicted values.

4. The system according to claim 1, characterized in that, The physical constraint loss builder includes a Biot equation tensor quantization difference operator and a residual scaling aligner; The tensor difference operator of the Biot equation constructs the elastic parameters of the solid skeleton and the fluid parameters of the flow phase into a fourth-order tensor, which is then substituted into the partial differential operator of the Biot equation. The automatic differentiation mechanism is used to calculate the difference recursion relationship between the solid displacement field and the flow phase displacement field on the spatiotemporal grid to obtain the theoretical synthesized wave field data. The residual scaling aligner extracts the envelope difference between the theoretically synthesized wavefield data and the seismic wavefield data, performs amplitude normalization and phase alignment operations on the envelope difference, and uses the envelope difference after amplitude normalization and phase alignment operations as the waveform residual, thereby compressing the dynamic range of the waveform residual.

5. The system according to claim 1, characterized in that, The joint optimizer includes an adaptive loss weighter and a gradient conflict resolver; The adaptive loss weighter calculates a dynamic balance coefficient based on the variance of the physical constraint loss term and the variance of the data-driven regression loss term, and uses the dynamic balance coefficient to perform a weighted summation of the physical constraint loss term and the data-driven regression loss term to obtain the total loss value. The gradient conflict resolver extracts the physical gradient direction corresponding to the physical constraint loss term and the data gradient direction corresponding to the data-driven regression loss term during the gradient backpropagation process based on the total loss value. It calculates the cosine value of the angle between the physical gradient direction and the data gradient direction. When the cosine value is less than 0, it performs orthogonal projection decomposition of the physical gradient direction and the data gradient direction on the shared parameters.

6. The system according to claim 4, characterized in that, The physical constraint loss builder also includes a solid-fluid coupling energy conservation corrector. The solid-fluid coupling energy conservation corrector calculates the sum of solid-phase kinetic energy, fluid-phase kinetic energy, solid-phase potential energy, and fluid-phase potential energy based on the solid-phase skeleton elastic parameters, the fluid-phase parameters, and the theoretically synthesized wavefield data to obtain the predicted total energy. Based on the seismic wavefield data, it calculates the observed total energy, calculates the energy deviation integral between the predicted total energy and the observed total energy, multiplies the energy deviation integral by the energy conservation penalty coefficient, and then adds it to the waveform residual to obtain the corrected waveform residual. The corrected waveform residual is determined as the physical constraint loss term.

7. The system according to claim 2, characterized in that, The multi-scale time-frequency orthogonal decomposer includes an adaptive wavelet basis selector and a multi-resolution cascaded expander. The adaptive wavelet basis selector extracts the spectral centroid and kurtosis of the seismic wavefield data, and matches the main wavelet basis function and auxiliary wavelet basis function in a preset wavelet basis library based on the spectral centroid and the kurtosis. The multi-resolution cascaded expander uses the main wavelet basis function to perform coarse-scale low-frequency approximation decomposition on the seismic wavefield data, and uses the auxiliary wavelet basis function to perform fine-scale high-frequency detail decomposition on the low-frequency components obtained by the coarse-scale low-frequency approximation decomposition. The results of the coarse-scale low-frequency approximation decomposition and the fine-scale high-frequency detail decomposition are concatenated along the channel dimension to form the multi-scale time-frequency coefficient matrix, which is combined to form a cross-band asymmetric wavelet basis decomposition structure.

8. The system according to claim 3, characterized in that, The solid-fluid decoupling mapping layer includes a pore connectivity graph builder and a graph convolution aggregator; The pore connectivity graph builder extracts the fluid spatial distribution probability matrix based on the slow longitudinal wave field features, determines the pixel nodes in the fluid spatial distribution probability matrix that are greater than the probability threshold as graph nodes, calculates the edge weights between adjacent graph nodes based on wave field coherence, and constructs a non-Euclidean graph structure. The graph convolution aggregator performs message passing on the non-Euclidean graph structure, projects the fast longitudinal wave field features and the transverse wave field features onto the node feature space of the non-Euclidean graph structure, aggregates the skeleton elastic response features of neighboring nodes in the non-Euclidean graph structure, and outputs the solid-phase decoupling features that include pore topological associations.

9. The system according to claim 4, characterized in that, The residual scaling aligner includes a wavelet phase matcher and an amplitude envelope normalizer. The wavelet phase matcher performs Hilbert transform on the theoretical synthesized wavefield data and the seismic wavefield data respectively to obtain the instantaneous phase, calculates the phase difference between the instantaneous phase of the theoretical synthesized wavefield data and the instantaneous phase of the seismic wavefield data, and performs time shift correction on the theoretical synthesized wavefield data based on the phase difference to obtain the phase-aligned wavefield. The amplitude envelope normalizer extracts the amplitude envelopes of the phase-aligned wavefield and the seismic wavefield data, calculates the logarithmic ratio of the amplitude envelope of the phase-aligned wavefield to the amplitude envelope of the seismic wavefield data, scales the amplitude of the phase-aligned wavefield based on the logarithmic ratio, obtains the amplitude-normalized wavefield, and determines the difference between the amplitude-normalized wavefield and the seismic wavefield data as the envelope difference after amplitude normalization and phase alignment operations.

10. A deep learning-based method for elastic inversion of two-phase media, characterized in that, include: Receive seismic wavefield data; Perform time-frequency orthogonal transformation and adaptive bandpass filtering on the seismic wavefield data to separate and extract fast P-wave wavefield features, slow P-wave wavefield features and shear wavefield features; Receive the fast longitudinal wave field characteristics, the slow longitudinal wave field characteristics, and the transverse wave field characteristics, and map and output the solid-phase skeleton elastic parameters and the fluid parameters of the two-phase medium; Substitute the solid-phase skeleton elastic parameters and the fluid parameters into the Biot two-phase medium wave equation for finite difference forward modeling to obtain theoretical synthetic wave field data. Calculate the waveform residual between the theoretical synthetic wave field data and the seismic wave field data, and convert the waveform residual into a physical constraint loss term. Obtain the physical constraint loss term and the data-driven regression loss term based on the labeled elasticity parameters. Sum the physical constraint loss term and the data-driven regression loss term with weights to obtain the total loss value. Perform gradient backpropagation based on the total loss value to update the weight parameters for performing the time-frequency orthogonal transformation, adaptive bandpass filtering operation, and mapping operation.

Images