An elastic wave multi-parameter synchronous inversion method based on a physical driving recurrent neural network and multi-observation system joint constraint

Abstract

The application provides an elastic wave multi-parameter synchronous inversion method and system based on a physical driving recurrent neural network and a multi-observation system joint constraint, belongs to the technical field of earth science and seismic exploration services, and comprises the following steps: 1, multi-parameter model construction; 2, observation system design and observation seismic record acquisition; 3, RNN elastic wave multi-parameter forward network design; 4, RNN elastic wave multi-parameter reverse gradient propagation design; 5, loss function design; and 6, multi-observation system joint inversion effect test. The method comprehensively utilizes the horizontal coverage advantage of the surface observation system and the constraint ability of the well observation system on the underground complex region and the deep region, enhances the constraint of the synchronous inversion of the underground shear wave velocity, longitudinal wave velocity and density model on the seismic data from different propagation paths and observation angles, and introduces multi-observation data weight distribution and smoothing regularization into the loss function to inhibit parameter crosstalk.

Description

Technical Field

[0001] This invention belongs to the field of Earth Science and Seismic Exploration Service Technology, specifically relating to a method for synchronous inversion of elastic waves based on physical-driven recurrent neural networks and joint constraints of multiple observation systems. Background Technology

[0002] In oil and gas exploration, reservoir characterization, and subsurface structural research, P-wave velocity (Vp), S-wave velocity (Vs), and density (rho) are crucial parameters describing the elastic properties of subsurface media. Compared to recovering only velocity parameters, the joint inversion of these three elastic parameters provides more complete physical property information for lithology identification, fluid prediction, and quantitative reservoir interpretation. Full-waveform inversion (FWI), using the wave equation as a physical constraint, minimizes the differences between simulated and observed seismic records, comprehensively utilizing the travel time, amplitude, phase, and different wave types of seismic waves. It has become one of the important methods for high-resolution subsurface parameter reconstruction.

[0003] However, elastic wave multiparameter FWI is an inverse problem with strong nonlinearity and parameter coupling characteristics. Seismic records are jointly controlled by Vp, Vs, and rho; perturbations of different parameters can easily produce similar or partially overlapping wavefield responses, leading to misallocation of data residuals during parameter updates, and consequently, inter-parameter error compensation and crosstalk. Density is generally less sensitive to the overall seismic data than velocity parameters and is easily affected by velocity parameter errors. Therefore, stably and synchronously recovering Vp, Vs, and rho under unified data residual constraints, especially improving the recovery capability of density background and interfaces, is a crucial issue in elastic wave multiparameter FWI.

[0004] To address the aforementioned issues, traditional elastic FWI research has improved the stability of multi-parameter simultaneous inversion through gradient correction, rock physics priors, and objective function design. For example, parameter contamination kernels are constructed based on multi-parameter Hessian off-diagonal blocks to quantify and suppress cross-contamination among the three parameters, and density recovery in three-parameter simultaneous inversion is improved through contamination term correction. Bayesian lithofacies estimation and rock physics prior constraints are introduced into FWI, and lithofacies estimation weights are adjusted according to the data sensitivity of different parameters to reduce the solution space of multi-parameter inversion and mitigate error compensation between parameters. However, traditional FWI typically requires explicit construction of the adjoint wave field and derivation of corresponding gradient expressions for different subsurface media parameters, making its gradient calculation and model parameter update processes complex.

[0005] Besides traditional FWI, data-driven deep learning methods have also been used for simultaneous multi-parameter prediction of elastic waves. These methods typically employ neural networks to learn the complex nonlinear mapping between seismic observation data and subsurface medium parameters, and utilize multi-task, multi-network collaboration, or multi-branch structures to mitigate parameter confusion and crosstalk, achieving simultaneous prediction of multiple parameters. However, these methods often lack explicit wave equation constraints and iterative matching of data residuals; their performance is highly dependent on the scale and distribution of the training data, and their generalization ability and interpretability remain limited.

[0006] To address the limitations of data-driven methods in terms of physical constraints and generalization capabilities, researchers have further proposed the FWI method based on physical-driven neural networks. Embedding the wave equation into a recurrent neural network (RNN) is a representative implementation. Since the time-recursive mechanism of the RNN is highly consistent with the numerical calculation process of time-step extrapolation of the seismic wavefield, each RNN unit can be associated with a wavefield time step, and the entire wavefield forward and backward gradient update process can be completed through recursive connections between units. Existing methods have verified the feasibility of this framework for simultaneous inversion of elastic wave three parameters and analyzed the impact of parameterization methods, objective functions, and optimization strategies on inversion performance and parameter crosstalk. However, further research and improvement are still needed in areas such as joint constraints of multiple observation systems, stable recovery of density parameters, and suppression of parameter crosstalk.

[0007] In summary, while existing multi-parameter full-waveform inversion methods for elastic waves can improve parameter coupling issues through gradient correction, prior constraints, and objective function design, traditional FWI methods are complex in terms of gradient derivation and model updating. Purely data-driven methods suffer from insufficient physical constraints, reliance on large-scale training samples, and limited generalization ability. The physics-driven RNN inversion framework can directly embed the elastic wave equation into the network, complete wavefield forward modeling through time recursion, and utilize automatic differentiation to calculate model parameter gradients, providing a more flexible approach for simultaneous multi-parameter inversion of elastic waves. However, relying solely on a single surface observation system is limited by the wavefield illumination range and propagation path, resulting in insufficient effective information in mid-deep and complex tectonic regions. This can easily lead to error compensation and parameter crosstalk between P-wave velocity, S-wave velocity, and density. The density parameter, due to its relatively weak data sensitivity, is particularly difficult to stabilize and recover.

[0008] Therefore, it is necessary to introduce joint constraints from surface and well observation data into the physical-driven RNN framework. By leveraging the better lateral coverage of surface observations and the stronger constraint capability of well observations on mid-deep and complex structures, wavefield information required for three-parameter inversion can be supplemented from different propagation paths. This improves the integrity of gradient information and parameter identification capability, thereby reducing parameter coupling and crosstalk, and achieving stable synchronous inversion of P-wave velocity, S-wave velocity, and density. Summary of the Invention

[0009] The purpose of this invention is to provide a multi-parameter synchronous inversion method. The main contributions include:

[0010] (1) The discrete elastic wave equation is embedded in the RNN framework. Multi-component seismic records are collected on the surface and in wells during the same wave field simulation process. The lateral coverage advantage of the surface observation system and the constraint ability of the well observation system on complex and deep underground areas are comprehensively utilized to enhance the constraint of seismic data on the synchronous inversion of underground Vp, Vs and rho from different propagation paths and observation angles.

[0011] (2) To address the differences in amplitude and energy between different observation systems and data components, this application constructs a grouped normalized weighted L1 data mismatch function, enabling various types of observation data to participate in gradient updates with controllable weights. At the same time, a model smoothing constraint is introduced into the loss function to suppress unstable local updates and high wavenumber artifacts.

[0012] (3) A two-stage inversion strategy from low frequency to high frequency is adopted. First, the large-scale background of the three parameters is restored, and then the interlayer interface and local structural details are gradually supplemented.

[0013] This invention provides the following technical solution: A method for simultaneous inversion of multiple parameters of elastic waves based on physical-driven recurrent neural networks and joint constraints of multiple observation systems includes the following steps:

[0014] Step (1): Multi-parameter model construction:

[0015] The construction of the multi-parameter model mainly includes the construction of the real model and the initial model: the real model is a complex two-dimensional isotropic elastic medium model containing various irregular strata, including P-wave velocity Vp, S-wave velocity Vs and density rho models; first, the three-parameter real model is downsampled to unify it to the same model size; then, smoothing is performed on the downsampled model to weaken the high wavenumber structure information of rock interface and local details, while retaining the overall velocity and density variation trend of the subsurface medium, thereby constructing the corresponding three-parameter real model and the initial model;

[0016] Step (2): Design of the observation system and acquisition of seismic records:

[0017] Based on the three-parameter real model constructed in step (1), a surface-multi-well joint observation system is designed:

[0018] First, multiple seismic sources are evenly distributed along the model surface, and geophones are placed on the surface at fixed intervals to obtain surface seismic records with a large lateral coverage.

[0019] Meanwhile, an observation well is set up on the left, middle and right sides of the model, and geophones are arranged at equal intervals along the depth direction of the well to record the wave field response of seismic waves at different depths and propagation paths underground. By combining the deployment of surface geophones and multi-well geophones, the underground medium is observed from different spatial locations and propagation directions to obtain various wave field information including direct waves, reflected waves, transmitted waves and converted waves.

[0020] Subsequently, the source was excited shot by shot according to the set source parameters. The actual P-wave velocity, S-wave velocity and density model in step (1) were used as the forward modeling medium model. The wave field propagation process under the excitation conditions of each source was calculated by elastic wave numerical forward modeling. The corresponding seismic waveforms were extracted synchronously at the locations of the geophones on the surface and in the well. Finally, the surface seismic record and multi-well observation seismic record were obtained for subsequent three-parameter inversion.

[0021] Step (3): Design of RNN elastic wave multi-parameter forward modeling network:

[0022] Based on step (2), RNN is used to generate predicted earthquake data; the first-order velocity-stress equation is embedded in the RNN framework, and the time recursion process of the elastic wave field is mapped to the cyclic propagation process of the RNN; each RNN unit takes the velocity field and stress field of the previous time step as input, and uses the velocity field and stress field as hidden states to recursively deduce along the time direction, and completes the wave field update by combining medium parameters, spatial derivatives and source information; and passes it to the next time step; elastic wave forward modeling is realized through continuous time step recursion, and the wave field response is extracted at the detector position to obtain the predicted earthquake record for subsequent model update;

[0023] Step (4): Design of multi-parameter back gradient propagation for RNN elastic waves:

[0024] Based on step (3), the wavefield extrapolation calculation is performed for each RNN unit corresponding to a time step, and all sampling time steps together constitute the time expansion depth of the RNN. As the RNN units are progressively extrapolated along the time direction, the complete elastic wave propagation process is represented as a continuous physical recursion network, and the corresponding wavefield state and its gradient information are saved at each time step for collecting seismic records and forming backpropagation links, providing a basis for seismic record extraction and backpropagation. Since the wavefield update operator, medium parameter mapping relationship and data acquisition process are differentiable, the gradient of the objective function with respect to the parameters of the model to be inverted can be obtained by the chain rule and automatic differentiation back calculation. Subsequently, the Adam optimization algorithm is used to iteratively update the model parameters based on the gradient, so that the predicted seismic record gradually approaches the observed record, thereby realizing the full waveform inversion of elastic waves based on RNN.

[0025] Step (5): Loss function design:

[0026] The loss function consists of a normalized multi-component data mismatch term and a smoothing regularization term. The data mismatch term measures the difference between the synthetic seismic record and the observed seismic record, and reduces the misallocation of data residuals among different parameters by balancing the contributions of different observation systems and data components, thereby reducing parameter crosstalk. The smoothing regularization term is used to suppress local anomalous updates and high wavenumber noise caused by crosstalk during multi-parameter synchronous inversion.

[0027] The specific steps in step (1) above are as follows: multi-parameter model downsampling and smoothing:

[0028] Based on the original three-parameter complex model, the original three parameters are scaled down to the same size by using the same downsampling factor down. Based on the scaled-down three-parameter model, a Gaussian blur with kernel x is used to create an initial three-parameter model. The prepared initial three-parameter model contains a shallow low-velocity, low-density region, a middle transition region, and a deep high-velocity, high-density region.

[0029] The specific steps in step (2) above are as follows: setting the seismic source, detector, and forward modeling parameters:

[0030] The n-shot seismic sources are evenly distributed on the ground surface at intervals of m meters, and the geophones are placed on the ground surface and in the wells at intervals of d meters, respectively. The wells are located at {x1, x2, x3} in the horizontal direction of the model. Then, forward modeling is performed using the two-dimensional isotropic first-order velocity-stress formula and the staggered grid C-order finite difference.

[0031] The specific steps for designing the RNN elastic wave multi-parameter forward modeling network in step (3) above are as follows:

[0032] The input to each RNN element is the velocity and stress components output by the previous element. These wavefield variables are recursively processed as hidden states along the time direction in the network. The spatial derivative term serves as an intermediate variable, which participates in the wavefield update after being calculated by the finite difference operator. The source term is injected into the normal stress field as an external excitation. The P-wave velocity Vp, S-wave velocity Vs, and density rho are parameters of the model to be inverted. These three are converted into Lamé coefficients and combined with the independent density term to participate in the recursive update of the velocity and stress fields. After completing the update of this time step, the element outputs new velocity and stress components and passes them to the next element, thereby realizing the continuous advancement of the seismic wavefield over time.

[0033] The total loss function in step (5) above is in the form of:

[0034] ;

[0035] in, For the total loss function, Represents weight, The data loss representing the Vz component of a certain observation system. Data loss representing the Vx component of a certain observation system, To smooth out the regularization loss;

[0036] The specific formula for the mismatch term in normalized multi-component data is:

[0037] ;

[0038] in, Represents the cannon number, Represents the number of cannons. Represents the sampling time. Represents the total sampling time. Represents earthquake records, It is a non-zero minimum constant; Represents data components, Representing the The first time network prediction Earthquake records of a certain velocity component from the No. 1 epicenter. Representing the The first time received The actual earthquake record of a certain velocity component of the source. They represent surface observation systems and well observation systems, respectively. Represents the well number, Representing the observation system, Vz and Vx represent the velocity components in the z and x directions, respectively. A mismatch function representing a velocity component of a certain observation system;

[0039] The specific form of the smoothing regularization term is:

[0040] ;

[0041] ;

[0042] in, and For spatial location index, and The grid point sizes are the parameters for the x and z directions, respectively. It is a first-order smoothing operator. The result of the first-order smoothing operator calculation for the model. To smooth out the regularization weights, For model parameters, Here are the parameter values ​​of the model at a certain point, where Vp, Vs, and ρ represent the transverse and longitudinal wave velocity models, the wave velocity model, and the density model, respectively. For the normalized result of the parameters, This represents the spatial average value of the corresponding parameter. , and These are the spatial average values ​​of P-wave velocity, S-wave velocity, and density, respectively. To smooth out the regularization loss.

[0043] The beneficial effects of this invention are as follows: This application provides a method for synchronous inversion of three parameters of elastic waves by combining the constraints of the physical equations of elastic waves, recurrent neural networks, and multiple observation systems. The aim is to construct a unified inversion framework that combines physical constraints and neural networks, and to reduce crosstalk problems in the synchronous inversion process by combining data collected from different observation systems and regularization constraints.

[0044] First, based on the RNN forward modeling framework driven by the discrete elastic wave equation, multi-component seismic records from the surface and wells are acquired during the same wavefield simulation process. By comprehensively utilizing the lateral coverage advantage of the surface observation system and the constraint capability of the well observation system for complex and deep underground regions, the constraints of seismic data on the synchronous inversion of subsurface Vp, Vs, and rho are enhanced from different propagation paths and observation angles.

[0045] Then, to address the amplitude and energy differences between different observation systems and data components, this paper constructs a grouped normalized weighted L1 data mismatch function, enabling various types of observation data to participate in gradient updates with controllable weights. At the same time, a model smoothing constraint is introduced into the loss function to suppress unstable local updates and high wavenumber artifacts.

[0046] Finally, a multi-parameter synchronous inversion experiment of elastic waves was conducted on the fuzzed Marmousi initial model to verify the effectiveness of the proposed method. Compared with inversion methods using only a single observation system, this invention effectively improves the accuracy and stability of multi-parameter synchronous inversion by fusing surface and well observation information. Furthermore, the RNN framework simplifies gradient calculation and model parameter update processes, reducing the implementation complexity of multi-parameter synchronous inversion and providing reliable technical support for further application in inversion of complex subsurface media and processing of actual seismic data. Attached Figure Description

[0047] Figure 1 This is the inversion flowchart of Example 1;

[0048] Figure 2 The diagram shows the Marmousi three-parameter real model and the three-parameter initial model prepared in Example 1, as well as the location distribution of the seismic source and the geophone. In the diagram: (a) is the real Vp model and the distribution of the well, seismic source and geophone; (b) is the initial Vp model; (c) is the real Vs model; (d) is the initial Vs model; (e) is the real rho model; (f) is the initial rho model.

[0049] Figure 3 These are surface and well seismic data records generated based on a real model in Example 1; in the figure: (a) surface seismic data record; (b) well seismic data record.

[0050] Figure 4 This is the RNN unit structure diagram of Example 1; in the diagram: the gray background represents the velocity components and their partial derivatives at the given time, the white background represents the stresses and their partial derivatives at the given time, and the dashed lines represent the calculation relationships between the physical quantities. The green background represents the trainable parameters (Vp, Vs, and rho), and the yellow background represents the Lamé coefficients;

[0051] Figure 5 This is Example 1, which is a synchronous inversion process of elastic wave three parameters based on RNN; in the figure: the upper part of the RNN unit generates surface seismic data records, and the lower part of the RNN unit generates well seismic data records.

[0052] Figure 6 This is a comparison chart of the inversion results of Example 2 and the inversion results using only the surface observation system; in the figure: (a): Vp inversion result of surface data, (b): Vp inversion result of the present invention (surface-multi-well joint), (c): Vs inversion result of surface data, (d): Vs inversion result of the present invention (surface-multi-well joint), (e): rho inversion result of surface data, (f): rho inversion result of the present invention (surface-multi-well joint).

[0053] Figure 7This is a comparison chart of the inversion results of Example 2 and the single-well observation system alone; in the figure: (a): Vp inversion result of single-well data, (b): Vp inversion result of the present invention (surface-single-well combined), (c): Vs inversion result of single-well data, (d): Vs inversion result of the present invention (surface-single-well combined), (e): rho inversion result of single-well data, (f): rho inversion result of the present invention (surface-single-well combined). Detailed Implementation

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

[0055] Example 1

[0056] As attached Figure 1 As shown in the flowchart, a method for simultaneous inversion of multiple parameters of elastic waves based on physical-driven recurrent neural networks and joint constraints of multi-observation systems includes the following steps:

[0057] Step (1): Multi-parameter model construction:

[0058] The multi-parameter model construction includes a three-parameter real model and a three-parameter initial model diagram. The three-parameter real model uses a two-dimensional complex isotropic elastic wave Marmousi model, including a transverse wave velocity model, a longitudinal wave velocity model, and a density model. The original three parameters are proportionally reduced to a 122x280 grid size by setting a downsampling factor down=20. Based on the reduced three-parameter model, a Gaussian blur with a kernel of x=12 is used to create a three-parameter initial model. The prepared three-parameter initial model contains a shallow low-velocity, low-density region, a middle transition region, and a deep high-velocity, high-density region. The resulting three-parameter real model and three-parameter initial model diagram are shown below. Figure 2 As shown, where, Figure 2 (a), (c) and (d) are three-parameter real models used for subsequent steps to acquire surface and well seismic data. Figure 2 (b), (d) and (e) are three-parameter initial models, which are used as initial values ​​in the subsequent inversion process.

[0059] Step (2): Design of the observation system and acquisition of seismic records:

[0060] Based on the real three-parameter model prepared in step (1), n=20 shot sources are evenly distributed on the ground surface at m=14 meters intervals, and geophones are placed on the ground surface and in wells at d=10 meters intervals, respectively. The wells are located at grid points {1, 140, 280} in the horizontal direction of the model. The distribution of the sources, geophones and wells is as follows. Figure 2As shown in (a). Subsequently, forward modeling was performed using the two-dimensional isotropic first-order velocity-stress formula and an 8th-order finite difference with staggered meshes to obtain a clear seismic record. The seismic record obtained in this step is as follows: Figure 3 As shown, it is used as an observation data record in the subsequent inversion process.

[0061] Step (3): Design of RNN elastic wave multi-parameter forward modeling network:

[0062] Based on step (2), to ensure that the wave field propagation process satisfies the elastic wave law, the first-order velocity-stress equation is embedded into the RNN inversion framework to guide the generation of predicted seismic records and the iterative update of model parameters during the forward modeling process. The input of each RNN unit is the velocity and stress components output by the previous unit. These wave field variables are recursively applied along the time direction as hidden states in the network; the spatial derivative term is used as an intermediate variable and participates in the wave field update after being calculated by the finite difference operator; the source term is injected into the normal stress field as an external excitation. The P-wave velocity Vp, S-wave velocity Vs, and density rho are used as the model parameters to be inverted. By converting the three into Lamé coefficients and combining them with the independent density term, they participate in the recursive update of the velocity and stress fields. The RNN forward modeling network unit structure is as follows: Figure 4 As shown, each RNN forward modeling network unit is equivalent to a time step. Within the unit, all parameters for the current time step are calculated, and the newly output velocity and stress components are passed to the next unit, thereby realizing the continuous advancement of the seismic wave field over time.

[0063] Step (4): Design of multi-parameter back gradient propagation for RNN elastic waves:

[0064] Based on step (3), one RNN unit corresponds to the wavefield extrapolation calculation within one time step, and the entire sampling time step can be regarded as the RNN network depth. Therefore, as the RNN unit gradually expands along the network depth direction, the complete elastic wave propagation process is represented as a physical recursive network. The network saves the wavefield state and its gradient information at each time step for collecting seismic records and forming backpropagation links, thus providing a basis for updating model parameters. Since the elastic wave forward modeling process consists of a differentiable wavefield update operator, medium parameter mapping relationship, and data acquisition system, the gradient of the objective function with respect to the parameters of the model to be inverted can be automatically calculated using the chain rule. After obtaining the gradient, the model parameters are iteratively updated using optimization algorithms such as Adam, so that the predicted seismic record gradually approaches the observed record, thereby realizing the elastic wave FWI backpropagation based on RNN. The RNN implementation process for elastic wave multi-parameter inversion described in this step is as follows: Figure 5 As shown;

[0065] Step (5): Loss function design:

[0066] The loss function consists of a normalized multi-component data mismatch term and a smoothing regularization term. The data mismatch term measures the difference between synthetic and observed seismic records, and reduces the misallocation of data residuals among different parameters by balancing the contributions of different observation systems and data components, thereby mitigating parameter crosstalk. The smoothing regularization term suppresses local anomalous updates and high wavenumber noise caused by crosstalk during multi-parameter synchronous inversion. Its overall form is:

[0067] .

[0068] in, For the total loss function, Represents weight, The data loss representing the Vz component of a certain observation system. Data loss representing the Vx component of a certain observation system, To smooth out the regularization loss, the weights of the surface-to-well VSP data mismatch function in this paper are: This means that the full value of the surface data is used in the loss function calculation, while the VSP data obtained from the three wells is compressed to 0.4 times before being used in the loss function calculation.

[0069] The specific formula for the mismatch term in normalized multi-component data is:

[0070] .

[0071] in, Represents the cannon number, Represents the number of cannons. Represents the sampling time. Represents the total sampling time. Represents earthquake records, It is a non-zero minimum constant. Represents data components, Representing the The first time network prediction Earthquake records of a certain velocity component from the No. 1 epicenter. Representing the The first time received The actual earthquake record of a certain velocity component of the source. They represent surface observation systems and well observation systems, respectively. Represents the well number, Representing the observation system, Vz and Vx represent the velocity components in the z and x directions, respectively. A mismatch function representing a velocity component of a certain observation system.

[0072] The specific form of the smoothing regularization term is:

[0073] ;

[0074] ;

[0075] ;

[0076] in, and For spatial location index, and The grid point sizes are the parameters for the x and z directions, respectively. It is a first-order smoothing operator. The result of the first-order smoothing operator calculation for the model. To smooth out the regularization weights, For model parameters, Here are the parameter values ​​of the model at a certain point, where Vp, Vs, and ρ represent the transverse and longitudinal wave velocity models, the wave velocity model, and the density model, respectively. For the normalized result of the parameters, This represents the spatial average value of the corresponding parameter. , and These are the spatial average values ​​of P-wave velocity, S-wave velocity, and density, respectively. To smooth the regularization loss, the smoothing regularization weights are: .

[0077] Example 2:

[0078] This embodiment utilizes the method of Embodiment 1 to compare the inversion quality of the proposed method with that of a single observation system when using only surface-single-well combined and surface-multi-well combined methods.

[0079] The experimental platform configuration for this embodiment is as follows: AMD EPYC 7773X, NVIDIA A100, Ubuntu 22.04 operating system, CUDA 11.7 acceleration, Python 3.9 experimental environment, and PyTorch 2.0.1 deep learning framework. The specific process of the example is as follows:

[0080] 2.1 Multi-parameter model construction:

[0081] The original Marmousi model's three-parameter .mat file was read into the system, and a consistency check was performed on each parameter file to ensure that the mesh size and spatial correspondence remained consistent. Based on this, the original Marmousi model was scaled down proportionally using the same scaling factor down=20, resulting in a three-parameter model with dimensions of 122x280. Then, a Gaussian blur with kernel x=12 was applied to blur it, creating an initial three-parameter model containing a shallow low-velocity, low-density region, a middle transition region, and a deep high-velocity, high-density region. The three-parameter model preparation process is as follows:

[0082] Step (1): Read the three-parameter file of the original Marmousi model. fields <- Load(longitudinal wave file, transverse wave file, density file).

[0083] Step (2): Check that all parameter field sizes are consistent. Check_Size_Consistency(fields).

[0084] Step (3): Construct a unified coordinate mapping based on the scaling factor. coords <- Build_Coordinates(original_size, down).

[0085] Step (4): Check if the correspondence between the multiple parameters is consistent. Check_Size_Consistency(resampled_field).

[0086] Step (5): Gaussian Blur. Initial_fields <- Gaussian Blur(resampled_field).

[0087] Step (6): Save the scaled-down Marmousi model and the initial model file. Save_To_Output(resampled_field, Initial_fields).

[0088] 2.2 Observation System Design and Seismic Record Acquisition:

[0089] Based on the realistic three-parameter model prepared in step 2.1, n=20 shot sources were uniformly distributed on the ground surface at intervals of m=14 meters, and geophones were placed at intervals of d=10 meters on the ground surface and in wells, with the wells located at grid points. Subsequently, forward modeling was performed using the two-dimensional isotropic first-order velocity-stress formula and an 8th-order finite difference method with staggered grids to obtain a clear seismic record. The process flow is as follows:

[0090] Step (1): Read the files of the three-parameter real model. fields <- Load(longitudinal wave file, transverse wave file, density file).

[0091] Step (2): Detector and true source settings. well_location <- (model horizontal width, well spacing): source[] <- allocation(number of shots, source spacing, detector spacing).

[0092] Step (3): Observation data acquisition. data_record<- forward(shot number, source, well_location, finite difference order).

[0093] 2.3 Design of RNN elastic wave multi-parameter forward modeling network:

[0094] Building upon step 2.2, to ensure the wavefield propagation process conforms to elastic wave laws, the first-order velocity-stress equation is embedded into the RNN inversion framework to guide the generation of predicted seismic records and the iterative updating of model parameters during forward modeling. The input to each RNN element is the velocity and stress components output from the previous element; these wavefield variables are recursively applied as hidden states along the time direction within the network. The spatial derivative term, as an intermediate variable, participates in the wavefield update after being calculated using a finite difference operator. The source term is injected into the normal stress field as an external excitation. The P-wave velocity Vp, S-wave velocity Vs, and density rho are used as parameters to be inverted in the model. By converting these three parameters into Lamé coefficients and combining them with an independent density term, they participate in the recursive updating of the velocity and stress fields. After completing the update at this time step, the element outputs new velocity and stress components and passes them to the next element, thus achieving continuous advancement of the seismic wavefield over time.

[0095] 2.4 RNN Elastic Wave Multi-parameter Backward Gradient Propagation Design:

[0096] One RNN unit corresponds to the wavefield extrapolation calculation within one time step, and the entire sampling time step can be regarded as the depth of the RNN network. Therefore, as the RNN units gradually expand along the network depth direction, the complete elastic wave propagation process is represented as a physical recursive network. The network saves the wavefield state and its gradient information at each time step to collect seismic records and form backpropagation links, thus providing a basis for updating model parameters. Since the elastic wave forward modeling process consists of a differentiable wavefield update operator, medium parameter mapping relationship, and data acquisition system, the gradient of the objective function with respect to the parameters of the model to be inverted can be automatically calculated using the chain rule. After obtaining the gradient, the model parameters are iteratively updated using optimization algorithms such as Adam, so that the predicted seismic record gradually approaches the observed record, thereby realizing the elastic wave FWI backpropagation based on RNN. The process is as follows:

[0097] Step (1): Set the trainable parameters: Vp, Vs, rho.

[0098] Step (2): Set forward modeling parameters, learning rate, etc.

[0099] Step (3): Multiple rounds of iteration.

[0100] Iteration begins (current round < maximum round):

[0101] Predicted data = RNN forward modeling (P-wave velocity, S-wave velocity, density);

[0102] Loss = Residual (predicted data, observed data) + Smoothing regularization;

[0103] Three-parameter model update (Adam, loss);

[0104] Iteration complete.

[0105] 2.5 Loss Function Design:

[0106] The loss function consists of a normalized multi-component data mismatch term and a smoothing regularization term. The data mismatch term measures the difference between synthetic and observed seismic records, and reduces the misallocation of data residuals among different parameters by balancing the contributions of different observation systems and data components, thereby mitigating parameter crosstalk. The smoothing regularization term suppresses local anomalous updates and high wavenumber noise caused by crosstalk during multi-parameter synchronous inversion. Its overall form is:

[0107]

[0108] in, For the total loss function, Represents weight, The data loss representing the Vz component of a certain observation system. Data loss representing the Vx component of a certain observation system, To smooth out the regularization loss, the weights of the surface-to-well VSP data mismatch function in this paper are: This means that the full value of the surface data is used in the loss function calculation, while the VSP data obtained from the three wells is compressed to 0.4 times before being used in the loss function calculation.

[0109] The specific formula for the mismatch term in normalized multi-component data is:

[0110]

[0111] in, Represents the cannon number, Represents the number of cannons. Represents the sampling time. Represents the total sampling time. Represents earthquake records, It is a non-zero minimum constant. Represents data components, Representing the The first time network prediction Earthquake records of a certain velocity component from the No. 1 epicenter. Representing the The first time received The actual earthquake record of a certain velocity component of the source. They represent surface observation systems and well observation systems, respectively. Represents the well number, Representing the observation system, Vz and Vx represent the velocity components in the z and x directions, respectively. A mismatch function representing a velocity component of a certain observation system.

[0112] The specific form of the smoothing regularization term is:

[0113] ;

[0114] ;

[0115] ;

[0116] in, and For spatial location index, and The grid point sizes are the parameters for the x and z directions, respectively. It is a first-order smoothing operator. The result of the first-order smoothing operator calculation for the model. To smooth out the regularization weights, For model parameters, Here are the parameter values ​​of the model at a certain point, where Vp, Vs, and ρ represent the transverse and longitudinal wave velocity models, the wave velocity model, and the density model, respectively. For the normalized result of the parameters, This represents the spatial average value of the corresponding parameter. To smooth the regularization loss, the smoothing regularization weights are: .

[0117] 2.6 Test of the joint inversion effect of multiple observation systems:

[0118] Under the same phase configuration and the same initial three-parameter model, two sets of comparative experiments were set up: one was to compare the inversion results of surface-multi-well joint observation data with the inversion results of surface observation data alone; the other was to compare the inversion results of surface-single-well joint observation data with the inversion results of intermediate well observation data alone, so as to verify the improvement effect of different joint observation methods on the synchronous inversion effect of three parameters.

[0119] 2.7 Implementation Results:

[0120] Figure 6 The results of inversion using only surface observation data and using surface-multi-well joint observation data proposed in this application are presented. When using only surface data, the three-parameter inversion results show significant differences due to insufficient wavefield illumination in the mid-deep and complex tectonic regions. Among them, Vp ( Figure 6 (a) in the middle can recover some of the main rock strata structure, but it is still accompanied by relatively obvious crosstalk noise; Vs ( Figure 6 (c) in the model can mainly distinguish the left side and shallow structure, while the recovery of the middle and deep complex structures is relatively blurry; rho ( Figure 6 The overall layered structure of (e) in the model is difficult to depict clearly. In contrast, after inversion using surface-multi-well joint observation data from this application, the three-parameter results ( Figure 6 The overall stratigraphic positions of (b), (d), and (f) are clearer, more continuous, and smoother. The parameter distribution in the shallow and deep regions is more stable, and the recovery of complex structures and inter-layer details is significantly improved. This indicates that the joint constraints of surface and multi-well observation data can effectively improve the overall quality of the three-parameter synchronous inversion.

[0121] Figure 7 The results of inversion using only intermediate well data and using surface-single-well joint observation data proposed in this invention are shown. (When using only intermediate well data...) Figure 7 In (a), (c), and (e) of the data, the Vp, Vs, and rho layered structures in the well-side region are relatively clearly recovered, but the inversion results in areas far from the wellbore are still poor, accompanied by significant crosstalk noise. In contrast, after using surface-single-well joint observation ( Figure 7 In (b), (d), and (f) of the model, the geological structures on both sides are clearer and more continuous, indicating that the good lateral coverage of surface data and the constraint of well data on the mid-to-deep region have a significant complementary effect. Joint observations further improved the recovery of layer interfaces, thin layers, and complex structures, with significant improvements in Vp and Vs, and the structural details of rho were also enhanced to some extent.

[0122] Overall, the surface-single-well joint observation data inversion method proposed in this invention can expand the effective wavefield illumination range, enhance the constraint ability of different regions on the three parameters, reduce the error compensation and crosstalk between parameters, thereby improving the stability and structural resolution of the three-parameter synchronous inversion.

Claims

1. A method for simultaneous inversion of multiple parameters of elastic waves based on physical-driven recurrent neural networks and joint constraints of multiple observation systems, characterized in that: Includes the following steps: Step (1): Multi-parameter model construction: The construction of the multi-parameter model mainly includes the construction of the real model and the initial model: the real model is a complex two-dimensional isotropic elastic medium model containing various irregular strata, including P-wave velocity Vp, S-wave velocity Vs and density rho models; first, the three-parameter real model is downsampled to unify it to the same model size; then, smoothing is performed on the downsampled model to weaken the high wavenumber structure information of rock interface and local details, while retaining the overall velocity and density variation trend of the subsurface medium, thereby constructing the corresponding three-parameter real model and the initial model; Step (2): Design of the observation system and acquisition of seismic records: Based on the three-parameter real model constructed in step (1), a surface-multi-well joint observation system is designed: First, multiple seismic sources are evenly distributed along the model surface, and geophones are placed on the surface at fixed intervals to obtain surface seismic records with a large lateral coverage. Meanwhile, an observation well is set up on the left, middle and right sides of the model, and geophones are arranged at equal intervals along the depth direction of the well to record the wave field response of seismic waves at different depths and propagation paths underground. By combining the deployment of surface geophones and multi-well geophones, the underground medium is observed from different spatial locations and propagation directions to obtain various wave field information including direct waves, reflected waves, transmitted waves and converted waves. Subsequently, the source was excited shot by shot according to the set source parameters. The actual P-wave velocity, S-wave velocity and density model in step (1) were used as the forward modeling medium model. The wave field propagation process under the excitation conditions of each source was calculated by elastic wave numerical forward modeling. The corresponding seismic waveforms were extracted synchronously at the locations of the geophones on the surface and in the well. Finally, the surface seismic record and multi-well observation seismic record were obtained for subsequent three-parameter inversion. Step (3): Design of RNN elastic wave multi-parameter forward modeling network: Based on step (2), RNN is used to generate predicted earthquake data; the first-order velocity-stress equation is embedded in the RNN framework, and the time recursion process of the elastic wave field is mapped to the cyclic propagation process of the RNN; each RNN unit takes the velocity field and stress field of the previous time step as input, and uses the velocity field and stress field as hidden states to recursively deduce along the time direction, and completes the wave field update by combining medium parameters, spatial derivatives and source information; and passes it to the next time step; Elastic wave forward modeling is achieved by recursively calculating continuous time steps, and wavefield response is extracted at the detector location to obtain the predicted earthquake record for subsequent model updates. Step (4): Design of multi-parameter back gradient propagation for RNN elastic waves: Based on step (3), the wavefield extrapolation calculation is performed for each RNN unit corresponding to a time step, and all sampling time steps together constitute the time expansion depth of the RNN. As the RNN units are progressively extrapolated along the time direction, the complete elastic wave propagation process is represented as a continuous physical recursion network, and the corresponding wavefield state and its gradient information are saved at each time step for collecting seismic records and forming backpropagation links, providing a basis for seismic record extraction and backpropagation. Since the wavefield update operator, medium parameter mapping relationship and data acquisition process are differentiable, the gradient of the objective function with respect to the parameters of the model to be inverted can be obtained by the chain rule and automatic differentiation back calculation. Subsequently, the Adam optimization algorithm is used to iteratively update the model parameters based on the gradient, so that the predicted seismic record gradually approaches the observed record, thereby realizing the full waveform inversion of elastic waves based on RNN. Step (5): Loss function design: The loss function consists of a normalized multi-component data mismatch term and a smoothing regularization term. The data mismatch term measures the difference between the synthetic seismic record and the observed seismic record, and reduces the misallocation of data residuals among different parameters by balancing the contributions of different observation systems and data components, thereby reducing parameter crosstalk. The smoothing regularization term is used to suppress local anomalous updates and high wavenumber noise caused by crosstalk during multi-parameter synchronous inversion.

2. The method for simultaneous inversion of multiple parameters of elastic waves based on physical-driven recurrent neural networks and multi-observation system constraints as described in claim 1, characterized in that: The specific steps in step (1) are as follows: multi-parameter model downsampling and smoothing: Based on the original three-parameter complex model, the original three parameters are scaled down to the same size by using the same downsampling factor down. Based on the scaled-down three-parameter model, a Gaussian blur with kernel x is used to create an initial three-parameter model. The prepared initial three-parameter model contains a shallow low-velocity, low-density region, a middle transition region, and a deep high-velocity, high-density region.

3. The method for simultaneous inversion of multiple parameters of elastic waves based on physical-driven recurrent neural networks and multi-observation systems according to claim 1, characterized in that: The specific steps in step (2) are as follows: setting the seismic source, detector, and forward modeling parameters: The n-shot seismic sources are evenly distributed on the ground surface at intervals of m meters, and the geophones are placed on the ground surface and in the wells at intervals of d meters, respectively. The wells are located at {x1, x2, x3} in the horizontal direction of the model. Then, forward modeling is performed using the two-dimensional isotropic first-order velocity-stress formula and the staggered grid C-order finite difference.

4. The method and system for simultaneous inversion of multiple parameters of elastic waves based on physical-driven recurrent neural networks and multi-observation systems according to claim 1, characterized in that: The specific steps for designing the RNN elastic wave multi-parameter forward modeling network in step (3) are as follows: The input to each RNN element is the velocity and stress components output by the previous element. These wavefield variables are recursively processed as hidden states along the time direction in the network. The spatial derivative term serves as an intermediate variable, which participates in the wavefield update after being calculated by the finite difference operator. The source term is injected into the normal stress field as an external excitation. The P-wave velocity Vp, S-wave velocity Vs, and density rho are parameters of the model to be inverted. These three are converted into Lamé coefficients and combined with the independent density term to participate in the recursive update of the velocity and stress fields. After completing the update of this time step, the element outputs new velocity and stress components and passes them to the next element, thereby realizing the continuous advancement of the seismic wavefield over time.

5. The method and system for simultaneous inversion of elastic waves based on physical-driven recurrent neural networks and multi-observation system constraints as described in claim 1, characterized in that: The total loss function in step (5) is in the form of: ； in, For the total loss function, Represents weight, The data loss representing the Vz component of a certain observation system. Data loss representing the Vx component of a certain observation system, To smooth out the regularization loss; The specific formula for the mismatch term in normalized multi-component data is: ； in, Represents the cannon number, Represents the number of cannons. Represents the sampling time. Represents the total sampling time. Represents earthquake records, It is a non-zero minimum constant; Represents data components, Representing the The first time network prediction Earthquake records of a certain velocity component from the No. 1 epicenter. Representing the The first time received The actual earthquake record of a certain velocity component of the source. They represent surface observation systems and well observation systems, respectively. Represents the well number, Representing the observation system, Vz and Vx represent the velocity components in the z and x directions, respectively. A mismatch function representing a velocity component of a certain observation system; The specific form of the smoothing regularization term is: ； ; ； in, and For spatial location index, and The grid point sizes are the parameters for the x and z directions, respectively. It is a first-order smoothing operator. The result of the first-order smoothing operator calculation for the model. To smooth out the regularization weights, For model parameters, Here are the parameter values ​​of the model at a certain point, where Vp, Vs, and ρ represent the transverse and longitudinal wave velocity models, the wave velocity model, and the density model, respectively. For the normalized result of the parameters, This represents the spatial average value of the corresponding parameter. , and These are the spatial average values ​​of P-wave velocity, S-wave velocity, and density, respectively. To smooth out the regularization loss.

Images