Robust full waveform inversion method with first arrival constraint physical embedded recurrent neural network

Abstract

The application discloses a kind of first arrival constraint physical embedded recurrent neural network robust full waveform inversion methods, first by the forward propagation of physical embedded recurrent neural network to solve wave equation, generate forward shot record, then use first arrival information to filter out the seismic trace with higher matching degree with observation data, participate in loss calculation, and gradually expand the seismic trace range participating in inversion along with the progress of inversion process, ensure that inversion result gradually tends to global optimal solution, and utilize the recorded space-time field gradient information of recurrent neural network forward propagation to carry out reverse propagation, further optimize velocity parameter, improve inversion precision.The method of the application embeds the physical process of seismic wave propagation into the recurrent neural network, introduces first arrival constraint and progressive inversion strategy, significantly reduces the dependence of FWI on the initial model, even in the case of significant deviation between the initial model and the true model, the underground velocity field structure can still be accurately reconstructed.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of seismic data processing, and particularly relates to a robust full waveform inversion method of first arrival constraint physical embedded recurrent neural network. BACKGROUND

[0002] Full-Waveform Inversion (FWI) is an important geophysical method for obtaining information of underground medium parameters. The core goal is to invert the medium parameters such as velocity and density of the underground by minimizing the difference between the observed data and the synthetic data. FWI relies on the numerical solution of the wave equation and can utilize the complete waveform information of seismic waves to provide high-resolution images of the underground structure. However, FWI is essentially a highly nonlinear optimization problem, and the accuracy of the inversion result is extremely sensitive to the accuracy of the initial velocity model. When the initial model is significantly different from the true model, the inversion process is easily trapped in local minimum, leading to the deviation of the inversion result from the true solution.

[0003] Existing FWI methods usually rely on low-frequency information to construct the initial model, or gradually invert from low frequency to high frequency through multi-scale strategy to reduce the influence of nonlinear problems. However, these methods face many challenges in practical applications. First, the acquisition of low-frequency information is often difficult in field exploration, especially in complex geological conditions, and the lack of low-frequency data will seriously affect the construction of the initial model. Second, although the multi-scale method can alleviate the nonlinear problem to some extent, its effect depends on the reasonable division of frequency bands, and the inversion process may still lead to local optimal solution due to the inaccuracy of the initial model.

[0004] In recent years, the rise of deep learning provides a new way to solve the nonlinear problem in FWI. Although the existing data-driven methods can train the model through a large amount of data, the results lack physical interpretability, and the requirements for data quality and scale are high. In order to overcome this problem, Physics-Informed Neural Networks (PINN) is introduced into FWI, which enhances the physical interpretability of the inversion result by embedding the wave equation and other physical constraints into the neural network. However, the existing PINN method still has limitations when dealing with inaccurate initial models, especially in the absence of prior information and low-frequency data, the inversion result may still be unsatisfactory. SUMMARY

[0005] To solve the above technical problems, the application provides a robust full waveform inversion method of first arrival constraint physical embedded recurrent neural network (FAFWI), which embeds the physical process of seismic wave propagation into a recurrent neural network (RNN), generates synthetic seismic data through physical-driven forward simulation, and combines the first arrival time difference constraint to dynamically screen seismic traces and gradually optimize the inversion result.

[0006] The technical scheme adopted by the application is as follows:

[0007] S1, a physical embedded recurrent neural network (PIRNN) framework based on first arrival constraint is constructed.

[0008] The framework comprises a first arrival information extraction module, a physical embedded recurrent neural network module, a first arrival constraint objective function construction module, and an inversion optimization module.

[0009] The physical embedded recurrent neural network module comprises an initial model, an SGFD RNN operator module, and the initial model is a velocity model.

[0010] S2, based on the framework constructed in step S1, the physical embedded recurrent neural network module maps the finite difference solution process of the wave equation to the forward propagation of the RNN, records the gradient information of the entire space-time field, and takes the velocity model parameters as trainable parameters in the network, to complete the physical constraint-based seismic forward modeling;

[0011] S3, the first arrival information extraction module extracts and utilizes the first arrival wave information from the observed data and the simulated data, constructs the first arrival constraint, and then the first arrival constraint objective function construction module constructs a dynamic objective function using the first arrival constraint, and only selects seismic traces with a first arrival time difference less than a time threshold to participate in loss calculation;

[0012] The first arrival constraint is constructed by calculating the cross-correlation function between the observed data and the simulated data, extracting the first arrival time difference, and dynamically selecting seismic traces participating in the inversion based on the first arrival time difference.

[0013] S4, based on the dynamic objective function constructed in step S3, the inversion optimization module updates the velocity model parameters through back propagation using the space-time field gradient information recorded by the forward propagation of the physical embedded recurrent neural network module, gradually optimizes the inversion result, and realizes robust full waveform inversion.

[0014] Further, step S2 is specifically as follows:

[0015] In the time domain, the initial model uses the first-order pressure-velocity acoustic wave equation for forward modeling, and its mathematical expression is as follows:

[0016]

[0017] where p represents the pressure, p x and p z represent the pressure components in the x and z directions respectively, v represents the acoustic velocity, v x and v z represent the particle velocities in the x and z directions respectively, v(r) represents the acoustic velocity to be inverted, and (r, t) represents the value of the physical quantity at position r at time t.

[0018] Then the velocity model is embedded as a trainable parameter into the RNN network, the wave equation is solved by using the staggered grid high-order finite difference method, and the perfect matching layer boundary condition is introduced, that is, the forward calculation is gradually performed through the SGFD RNN operator module.

[0019] Discretize equation (1), that is, define the pressure p x , p z on the integer grid points, and define the particle velocities v x , v z on the half grid points. For the second-order time accuracy and the second-order finite difference, expand at time and [ixΔt, izΔz] to obtain the following expression:

[0020]

[0021] where k represents the index of the time step, ix and iz represent the discrete number indexes of the spatial grid in the x and z directions respectively, p k represents the discrete value of the pressure p at the kth time step, Δx represents the discrete step length in the x direction on the finite difference space, Δz represents the discrete step length in the z direction on the finite difference space, Δx and Δz are also called grid spacing; and Δt represents the discrete step length in time, that is, the interval between two adjacent time points in the grid.

[0022] Further process the discrete equations of equations (2)-(5) to obtain the recursive formulas about p x , p z , v x , v z , and the expression is as follows:

[0023]

[0024] Then solve the wave field state at the next time step by time step method, that is, the forward modeling of the seismic wave field can be performed through the RNN.

[0025] Each layer of the RNN is used to calculate and store the wave field information at a certain time, and each RNN unit contains a chain relationship about the velocity parameter v. During the forward process, the neural network records the chain relationship of the entire space-time field, and calculates the gradient of the loss function with respect to the velocity parameter v according to the automatic differentiation mechanism during the back propagation, and updates the velocity.

[0026] Further, the step S3 is specifically as follows:

[0027] During the inversion process, first, for each seismic trace, the cross-correlation function between the observed data and the simulated data is calculated The expression is as follows:

[0028]

[0029] Where τ represents the time delay, s represents the seismic trace, and t represents the time. The peak position of the cross-correlation function corresponds to the best time alignment position between the observed data and the simulated data, and the expression is as follows:

[0030]

[0031] Then, the first arrival time difference is extracted, and the first arrival time difference is extracted through the peak position of the cross-correlation function The expression is as follows:

[0032]

[0033] Where T represents the total number of time points, Δt s represents the time phase difference between the observed data and the simulated data corresponding to the s seismic trace.

[0034] Based on the first arrival time difference Δt s , the seismic trace participating in the loss calculation is dynamically selected, and the first arrival constraint is constructed, that is, a time threshold Δt threshold is set, and when |Δt s |<Δt threshold , it is considered that the underground velocity information carried by the seismic trace is well matched with the current initial model; otherwise, the seismic trace will be temporarily excluded from the loss calculation. That is, only the seismic trace that satisfies |Δt s |<Δt threshold is selected to participate in the inversion, and the expression is as follows:

[0035]

[0036] Where, represents the set of all seismic traces, denotes the selected seismic trace set.

[0037] The first arrival constraint target function construction module utilizes the first arrival constraint to construct a dynamic target function, only selects seismic traces with a first arrival time difference less than a time threshold to participate in loss calculation, and then the dynamic target function J(v) based on the first arrival constraint is expressed as follows:

[0038]

[0039] wherein, denotes the selected seismic trace set, denotes the simulated data of the s-th seismic trace at time t, denotes the corresponding observed data, and T denotes the time step number.

[0040] Further, the step S4 is specifically as follows:

[0041] In the inversion process, a progressive strategy is adopted, and the number of seismic traces satisfying the first arrival time difference threshold is The inversion process is dynamically adjusted.

[0042] wherein, I denotes an indicator function, which is used to select appropriate seismic traces, and Shots denotes the total number of seismic traces.

[0043] Then, iterative optimization is performed, in each iteration process, steps S2-S3 are repeated, the initial model is forward calculated through the SGFDRNN operator module to obtain simulated data and record the gradient information of the entire space-time field, then the simulated data and the observed data are screened according to the first arrival constraint, the loss of the screened seismic traces is calculated, and the initial model is updated by back propagation, and the process is repeated until most seismic traces are selected to participate in the operation, the total first arrival difference is close to 0, and the loss function converges to a preset threshold, thereby realizing robust full waveform inversion.

[0044] wherein, the inversion optimization module updates the velocity model parameter v through the formula γ denotes a learning rate, k denotes an iteration number, denotes the gradient of the loss function with respect to the velocity parameter.

[0045] The beneficial effects of this invention are as follows: The method of this invention first solves the wave equation through forward propagation of a physically embedded recurrent neural network to generate forward shot gather records. Then, it uses first arrival information to select seismic traces with high matching degree with the observation data to participate in loss calculation, avoiding waveform mismatch problems caused by inaccurate initial models. As the inversion process progresses, the range of seismic traces participating in the inversion is gradually expanded to ensure that the inversion results gradually tend to the global optimal solution. Furthermore, the spatiotemporal field gradient information recorded by the forward propagation of the recurrent neural network is used for back propagation to further optimize velocity parameters and improve inversion accuracy. The method of this invention embeds the physical process of seismic wave propagation into a recurrent neural network, introduces first-arrival constraints and a progressive inversion strategy, effectively avoids the periodic jump problem caused by inaccurate initial models in existing full-waveform inversion, significantly reduces the dependence of FWI on the initial model, and can still accurately reconstruct the underground velocity field structure even when there are significant deviations between the initial model and the real model. Moreover, the first-arrival constraint strategy can screen out seismic traces with high matching degree with the initial model in the early stage of inversion, avoid error accumulation, and ensure stable convergence of the inversion process. At the same time, the progressive strategy gradually introduces more seismic traces to participate in the inversion, ensuring that the inversion results gradually tend to the global optimum, and significantly improves the robustness of the inversion results. Attached Figure Description

[0046] Figure 1 This is a flowchart of a robust full waveform inversion method for first-arrival constraint physical embedded recurrent neural network according to the present invention.

[0047] Figure 2 This is a flowchart of the PIRNN forward modeling process in an embodiment of the present invention.

[0048] Figure 3 This is a schematic diagram of the real velocity model and the constant velocity model in the embodiments of the present invention.

[0049] Figure 4 This is a schematic diagram of the inversion results in an embodiment of the present invention.

[0050] Figure 5 This is a schematic diagram of the velocity profile in an embodiment of the present invention.

[0051] Figure 6 This is a schematic diagram showing the loss function, total phase difference Δφ, and number of selected seismic traces on the Marmousi model in this embodiment of the invention, as a function of the method of the present invention, and the curves of these parameters as a function of the number of iterations. Detailed Implementation

[0052] The method of the present invention will be further described below with reference to the accompanying drawings and embodiments.

[0053] like Figure 1 The flowchart of a robust full-waveform inversion method for first-arrival constraint physical embedded recurrent neural network according to the present invention is shown below. The specific steps are as follows:

[0054] S1, constructing a physical embedding recurrent neural network PIRNN framework based on first arrival constraints;

[0055] The framework comprises a first arrival information extraction module, a physical embedding recurrent neural network module, a first arrival constraint objective function construction module, and an inversion optimization module.

[0056] The physical embedding recurrent neural network module comprises an initial model and an SGFD RNN operator module. The initial model is a velocity model.

[0057] S2, based on the framework constructed in step S1, the physical embedding recurrent neural network module maps the finite difference solving process of the wave equation to the forward propagation of the RNN, records the gradient information of the entire space-time field, and takes the velocity model parameters as trainable parameters in the network, to complete the physical constraint-based seismic forward modeling;

[0058] S3, based on the first arrival information extraction module, first arrival wave information is extracted and utilized from the observed data and the simulated data, a first arrival constraint is constructed, and then the first arrival constraint objective function construction module constructs a dynamic objective function using the first arrival constraint, and only selects seismic traces with a first arrival time difference less than a time threshold to participate in loss calculation;

[0059] Wherein, the first arrival constraint is constructed by calculating the cross-correlation function between the observed data and the simulated data, extracting the first arrival time difference, and dynamically selecting seismic traces participating in inversion based on the first arrival time difference.

[0060] S4, based on the dynamic objective function constructed in step S3, the space-time field gradient information recorded by the forward propagation of the physical embedding recurrent neural network module, the inversion optimization module updates the velocity model parameters through back propagation, gradually optimizes the inversion result, and realizes robust full waveform inversion.

[0061] In this embodiment, step S2 is specifically as follows:

[0062] First, the seismic forward embedding recurrent neural network is described, and accurate seismic forward modeling is the basis for full waveform inversion. Compared with existing full waveform inversion, the physical embedding recurrent neural network (Physical-Informed Recurrent Neural Network, PIRNN) can seamlessly integrate physical constraints such as wave equations into the neural network structure. In this embodiment, in the time domain, the initial model uses the first-order pressure-velocity acoustic wave equation for forward modeling, and its mathematical expression is as follows:

[0063]

[0064] Wherein, p represents pressure, px and p z denote the pressure components in x and z directions, v denotes the sound speed, v x and v z denote the particle velocities in x and z directions, v(r) denotes the sound speed to be inverted, (r,t) denotes the value of a physical quantity at r position at t time.

[0065] As Figure 2 shown, in order to ensure high accuracy of the forward process, the velocity model is embedded into the RNN network as a trainable parameter (mapping the time dependence of the wave equation into the hierarchical structure of recurrent neural network (RNN)), the wave equation is solved by using the staggered grid high-order finite difference method, and the perfect matched layer (PML) boundary condition is introduced to effectively absorb the outgoing wave and reduce the boundary reflection, that is, the forward calculation is carried out step by step through the SGFD RNN operator module (Staggered Grid Finite Difference Recurrent Neural Network Operator).

[0066] Figure 2 The workflow of PIRNN forward modeling is shown. Among them, Input(i) represents the wave field information at the i-th time step, and Output(i) represents the simulated shot record generated at the i-1 time step. Through this structure, the present application can efficiently realize seismic forward modeling and provide reliable physical constraints and gradient information for full waveform inversion.

[0067] In this process, the underground medium parameters are embedded into each SGFD RNN operator as trainable parameters, and the gradient information of these parameters is included in the loss calculation. This enables the neural network to calculate the gradient of the objective function on the underground parameters through back propagation and iteratively update these parameters. In this way, the present embodiment realizes the seismic forward modeling based on physical constraints and lays a foundation for subsequent full waveform inversion.

[0068] First, formula (1) is discretized, that is, the pressure p x , p z is defined on the integer grid point, the particle velocity v x , v z is defined on the half grid point, for the second-order time accuracy and the second-order finite difference, the expansion is carried out at time and [ixΔt,izΔz], and the expression is as follows:

[0069]

[0070] where k denotes the index of time step, ix, iz denote the discrete index of spatial grid in x and z direction respectively, p k denotes the discrete value of pressure p at the kth time step, Δx denotes the discrete step in x direction in finite difference space, Δz denotes the discrete step in z direction in finite difference space, Δx and Δz are also called grid spacing; Δt denotes the discrete step in time, i.e. the interval between two adjacent time points in the grid.

[0071] Further processing the discrete equations of formula (2)-(5) obtains the recursive formula about p x , p z , v x , v z , the expression is as follows:

[0072]

[0073]

[0074] Then the wave field state at the next time is solved in turn according to the time stepping method, i.e. the forward modeling of seismic wave field is performed through RNN.

[0075] The recurrent neural network (RNN) is a kind of artificial neural network capable of processing time series data, and its core feature is to use internal state (memory) to process sequence input, so that the current output can be affected by previous results. This feature makes RNN particularly suitable for processing time-dependent signal processing tasks. In the numerical simulation of wave equation, the time stepping method is widely used to calculate the change of wave field state. The hidden state of RNN can carry the information of the previous time, which is highly consistent with the time stepping method in the solution of wave equation. Specifically, each layer of RNN is used to calculate and store the wave field information at a certain time, including pressure and particle velocity, and each RNN unit contains a chain relationship about the velocity parameter v. In the process of forward modeling, the neural network records the chain relationship of the whole space-time field, and calculates the gradient of the loss function with respect to the velocity parameter v according to the automatic differentiation mechanism during back propagation, and updates the velocity. This structure allows the time evolution process of the wave field to be mapped into the hierarchical structure of RNN, and each layer not only stores the wave field information at the current time, but also receives the output of the previous layer as input, thereby well simulating the time dependence of the wave equation.

[0076] In the embodiment, the step S3 is specifically as follows:

[0077] In view of the problem of waveform matching period jump in full waveform inversion (FWI), the embodiment introduces a first arrival constraint strategy to alleviate the nonlinearity problem. The first arrival wave is the wave that arrives at the receiver first in the process of seismic wave propagation, and is usually located near the time point with the strongest energy of the seismic trace. Under a vertical seismic profiling (VSP) observation system, the time arrival characteristics of the first arrival wave can provide key information about the velocity structure. By extracting and utilizing the first arrival wave information, the present application can effectively reduce the influence of the initial velocity model on the inversion result. In the VSP observation system, the time difference Δt of the first arrival wave is an important indicator for judging the matching degree of the observation data and the synthetic data. When the initial velocity model is greatly different from the true model, the time difference of the first arrival wave may exceed half a period, resulting in the period jump phenomenon of the FWI in waveform matching, and further causing the model update to fall into a local optimal solution. To avoid this problem, the embodiment proposes a first arrival constraint strategy, the core of which is to only select the seismic trace with a first arrival time difference Δt less than half a period to participate in the construction and optimization of the objective function. Through this strategy, the period jump phenomenon can be effectively avoided, and the model parameters can be guided to iteratively update in the correct direction.

[0078] Specifically, the embodiment extracts the first arrival time difference by calculating the cross-correlation function between the observation data and the synthetic data.

[0079] In the inversion process, first, for each seismic trace, the cross-correlation function between the observation data and the simulation data is calculated. The expression is as follows:

[0080]

[0081] Where τ represents the time delay, s represents the seismic trace, and t represents the time. The peak position of the cross-correlation function corresponds to the best time alignment position between the observation data and the simulation data, and the expression is as follows:

[0082]

[0083] Then, the first arrival time difference is extracted by the peak position of the cross-correlation function , and the expression is as follows:

[0084]

[0085] Where T represents the total number of time points, and Δt s represents the time phase difference between the observation data and the simulation data corresponding to the s seismic trace.

[0086] Based on the first arrival time difference Δts , dynamically select seismic traces involved in loss calculation, construct first arrival constraint, i.e. set a time threshold Δt threshold , when |Δt s | < Δt threshold , it is considered that the underground velocity information carried by the seismic trace is well matched with the current initial model; otherwise, the seismic trace will be temporarily excluded from the loss calculation. That is, only seismic traces that meet |Δt s | < Δt threshold participate in inversion, and the expression is as follows:

[0087]

[0088] wherein, denotes the set of all seismic traces, denotes the set of selected seismic traces. This dynamic selection mechanism can ensure that, in the early stage of inversion, when the initial model is greatly different from the true model, only seismic traces with small first arrival time difference participate in loss calculation, thereby effectively reducing error accumulation and improving the stability of inversion. With iterative optimization of model parameters, the number of seismic traces that meet the first arrival time difference threshold condition will gradually increase, thereby providing more abundant underground medium structure information, and finally converging to the global optimal solution.

[0089] The first arrival constraint objective function construction module constructs a dynamic objective function using the first arrival constraint, and only selects seismic traces with a first arrival time difference less than a time threshold to participate in loss calculation. Then, the expression of the dynamic objective function J(v) based on the first arrival constraint is as follows:

[0090]

[0091] wherein, denotes the set of selected seismic traces, denotes the simulation data of the s-th seismic trace at time t, denotes the corresponding observation data, and T denotes the time step.

[0092] Through this dynamic objective function, the embodiment can gradually introduce more seismic traces in the inversion process, ensuring that the inversion result starts from reliable data and gradually improves the accuracy and comprehensiveness.

[0093] In the embodiment, the step S4 is specifically as follows:

[0094] In the inversion process, a progressive strategy is adopted, and only seismic traces with a first arrival time difference less than a time threshold are selected to participate in inversion in the early stage. With the progress of inversion, the number of seismic traces participating in inversion is gradually increased, ensuring that the inversion result gradually tends to be globally optimal. That is, the number of seismic traces that meet the first arrival time difference threshold is gradually increased in the inversion process. The inversion process is dynamically adjusted.

[0095] where I denotes the indicator function for selecting appropriate seismic traces, and Shots denotes the total number of seismic traces.

[0096] Then iterative optimization is performed, in each iteration, steps S2-S3 are repeated, the initial model is forward modeled by the SGFDRNN operator module to obtain the simulated data and record the gradient information of the entire space-time field, then the simulated data and the observed data are screened according to the first arrival constraint, the selected seismic traces are calculated to obtain the loss, and the initial model is updated by back propagation, and the process is repeated until most of the seismic traces are selected for operation, the total first arrival difference is close to 0, and the loss function converges to the preset threshold, and robust full waveform inversion is realized.

[0097] where the inversion optimization module updates the velocity model parameter v by the formula , γ denotes the learning rate, k denotes the iteration number, represents the gradient of the loss function with respect to the velocity parameter.

[0098] The embodiment is further verified by experiments, as shown in Figure 3 Fig. 2, a comparison is made between the original FWI and the method of the present application in terms of initial model dependence in the Marmousi model. The constant velocity model is used as the initial model for full waveform inversion, and the initial model has a significant difference from the true velocity model, without low frequency information and well logging information, and therefore is a challenging initial model for FWI. In addition, in order to ensure the consistency of the experiment, the Adam optimization algorithm is used for different FWI methods, the time threshold Δt threshold is set to 2.5 ms, and the learning rate is uniformly set to 40 and remains unchanged.

[0099] Figure 3 (a) is the setting of the true velocity model and the observation system, the red pentagram indicates the placement position of the seismic source, and the black dashed line indicates the placement position of the geophone, Figure 3 (b) is a constant velocity model with a velocity value of 1800 m / s, Figure 3 (b) shows the constant velocity initial model used in the Marmousi inversion, and the velocity value of each point is 1800 m / s. If no prior structural information is added as a constraint, such a velocity model is easy to fall into a local minimum value and obtain an incorrect physical explanation under the original FWI.

[0100] Figure 4 is a schematic diagram of the inversion result, as shown in Figure 4 (b), the inversion result under the original FWI presents a high degree of uncertainty and completely deviates from the true velocity structure. Compared with the original FWI, the inversion result obtained by the method proposed in the present application is shown in Figure 4(a) as shown, it can be seen that under the same initial model, the FAFWI method of the present application can still obtain more accurate inversion results within the effective range of the VSP observation system, not only reconstructing the velocity structure of the shallow and middle layers, but also better inverting the velocity structure and change trend in the deep layer region.

[0101] As Figure 5 shown, in order to quantitatively evaluate the velocity field characteristics of the wellbore adjacent area, the vertical velocity profile with the wellbore as the center and the distance of ±75m is extracted. Among them Figure 5 (a) is the velocity profile 75m away from the left side of the wellbore under the method (FAFWI) of the present application, Figure 5 (b) is the velocity profile 75m away from the right side of the wellbore under the method (FAFWI) of the present application, and the profile comparison analysis shows that the inversion result has high consistency with the true model in most areas from the shallow layer to the deep layer.

[0102] From the spatial distribution characteristics of the velocity field inversion result, the effective inversion area presents a clear inverted triangular shape, which is due to the fact that the geophones in the VSP observation system are vertically arranged along the wellbore, and the seismic reflection waves in the lateral area far away from the shot-hole pair are difficult to be effectively received, so the inversion result conforms to the physical law. On this basis, the R2 score, the structure similarity index SSIM and the normalized correlation coefficient of the inversion result in this range and the true model are calculated to comprehensively evaluate the inversion result obtained by the method of the present application in terms of linear correlation, structural correlation and the like. The results obtained are shown in Table 1 (quantitative comparison of the inversion results of the FAFWI method on the Marmousi model). The values in Table 1 show that the model inverted by the method of the present application presents strong linear relationship and structural similarity with the true model, and has an advantage in the normalized correlation.

[0103] Table 1

[0104]

[0105] As Figure 6 shown, with the continuous progress of the inversion process, the velocity parameters will be updated to the correct direction, so in each iteration, the number of seismic traces that satisfy |Δt s |<Δt threshold will be more and more, and the total phase difference will gradually tend to 0, and the number of selected seismic traces will show a gradually rising trend. This dynamic adjustment mechanism realizes a dynamic objective function. With the progress of the inversion, more seismic traces are allowed to participate in the inversion, so that the inversion result is more comprehensive and accurate. Among them, Figure 6 (a) is a dynamic loss curve diagram, Figure 6(b) a plot of the trend of the total phase difference Δφ for all seismic traces, Figure 6 (c) a plot of the trend of the number of seismic traces that satisfy the time threshold Δt threshold .

[0106] Model tests and actual data tests prove that the FAFWI full waveform inversion network structure using the method of the present application significantly reduces the dependence on the initial model under the premise of ensuring a certain inversion accuracy. This proves that the method of the present application effectively combines physical constraints with neural networks and first arrival constraints, reduces the sensitivity to the initial model, and ensures the accuracy of full waveform inversion.

[0107] In summary, through numerical experiments, the method of the present application can accurately reconstruct the subsurface velocity field structure even if there is a significant deviation between the initial model and the true model, significantly reducing the dependence of FWI on the initial model. By introducing first arrival wave information, the method dynamically selects seismic traces with high matching degree to the initial model to participate in inversion, effectively reducing the error accumulation caused by inaccurate initial models. At the same time, by combining the spatiotemporal field gradient information of recurrent neural network (RNN), the finite difference solution of the wave equation is realized, further improving the stability and accuracy of the inversion.

[0108] Those skilled in the art will realize that the embodiments described herein are for the purpose of aiding the reader in understanding the principles of the present application and should be construed as not limiting the scope of the present application to such specific details and embodiments. The present application can have various modifications and changes for those skilled in the art. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principles of the present application should be included in the scope of the claims of the present application.

Claims

1. A robust full-waveform inversion method for first-arrival constraint-based physical embedding of a recurrent neural network, the specific steps of which are as follows: S1. Construct the Physical Embedded Recurrent Neural Network (PIRNN) framework based on initial arrival constraints; The framework includes: The module includes: initial arrival information extraction module, physical embedding recurrent neural network module, initial arrival constraint objective function construction module, and inversion optimization module. The physical embedded recurrent neural network module includes: an initial model and an SGFD RNN operator module; the initial model is the velocity model. S2. Based on the framework constructed in step S1, the physical embedded recurrent neural network module maps the finite difference solution process of the wave equation to the forward propagation of the RNN, records the gradient information of the entire spatiotemporal field, and uses the velocity model parameters as trainable parameters in the network to complete the seismic forward modeling based on physical constraints. S3. The first arrival information extraction module extracts and utilizes first arrival information from observation data and simulation data to construct first arrival constraints. Then, the first arrival constraint objective function construction module constructs a dynamic objective function using the first arrival constraints, and only selects seismic traces with first arrival time differences less than the time threshold to participate in loss calculation. Among them, constructing first arrival constraints involves calculating the cross-correlation function between observed data and simulated data, extracting the first arrival time difference, and dynamically selecting seismic traces to participate in the inversion based on the first arrival time difference; S4. Based on the dynamic objective function constructed in step S3, the spatiotemporal field gradient information recorded by the forward propagation of the physically embedded recurrent neural network module is used. The inversion optimization module updates the velocity model parameters through back propagation, gradually optimizes the inversion results, and achieves robust full waveform inversion. Step S4 is as follows: In the inversion process, a progressive strategy is adopted, which is determined by the number of seismic traces that meet the first arrival time difference threshold. The inversion process is dynamically adjusted. in, This indicates the indicator function used to select appropriate seismic traces; Shots represents the total number of seismic traces. Indicates the first The time phase difference between the observed data and the simulated data corresponding to each seismic trace; This indicates the set time threshold; Then, iterative optimization is performed. In each iteration, steps S2-S3 are repeated. The initial model is forward modeled through the SGFD RNN operator module to obtain the simulation data and record the gradient information of the entire spatiotemporal field. Then, the simulation data and observation data are screened for seismic traces according to the first arrival constraint. The loss of the screened seismic traces is calculated and backpropagated to update the initial model. This process is repeated until most seismic traces are selected to participate in the calculation, the total first arrival difference is close to 0, and the loss function converges to the preset threshold, thus achieving robust full waveform inversion. The inversion optimization module uses the formula Update speed model parameters , Indicates the learning rate. Indicates the number of iterations. This represents the gradient of the loss function with respect to the velocity parameter.

2. The robust full-waveform inversion method for first-arrival constraint-embedded recurrent neural networks according to claim 1, characterized in that, Step S2 is as follows: In the time domain, the initial model is modeled using the first-order pressure-velocity acoustic wave equation, and its mathematical expression is as follows: （1）； in, Indicates pressure. and They represent and Pressure component in the direction, Indicates the speed of sound. and They represent and The velocity of the particle in the direction of [direction] This indicates the velocity of the sound wave that needs to be inverted. Indicates the physical quantity in time The value at the location; Then, the velocity model is embedded into the RNN network as a trainable parameter, and the wave equation is solved by the staggered grid high-order finite difference method. The boundary conditions of the perfectly matched layer are introduced, that is, the forward calculation is performed step by step through the SGFD RNN operator module. Discretize equation (1), that is, pressure , Defined on integer grid points, the velocity of a particle , Defined at half-grid points, for second-order time precision and second-order finite difference, in Time Expanding at this point, we get the following expression: （2）； （3）； （4）； （5）； in, Indicates the index of the time step. They represent and Discrete indexing of spatial grids in orientation. Indicates the first At each time step, the pressure discrete values, Denotes the dispersion measure in the x-direction on a finite difference space. Denotes the distance step size in the z-direction on the finite difference space. and Also known as grid spacing; It represents the distance between two adjacent time points in the grid. Further processing of the discrete equations (2)-(5) yields information about... , , , The recursive formula is expressed as follows: （6）； （7）； （8）； （9）； Then, by solving the wave field state at the next time step in sequence using the time step method, the forward modeling of the seismic wave field can be performed using RNN; Each layer of the RNN is used to calculate and store wavefield information at a certain moment, and each RNN unit contains information about velocity parameters. The neural network records the chain-like relationships of the entire spatiotemporal field during forward propagation, and calculates the loss function with respect to the velocity parameters during backpropagation based on the automatic differentiation mechanism. The gradient is calculated, and the velocity is updated.

3. The robust full-waveform inversion method for first-arrival constraint-embedded recurrent neural networks according to claim 1, characterized in that, Step S3 is as follows: In the inversion process, the observation data are first calculated for each seismic trace. With simulation data Cross-correlation function between The expression is as follows: （10）； in, Indicates a time delay. Indicates earthquake path, Indicates time; peak position of the cross-correlation function The expression corresponding to the optimal time alignment between observed and simulated data is as follows: （11）； Then, the first arrival time difference is extracted, and the peak position of the cross-correlation function is used. The expression for extracting the difference between arrival and departure times is as follows: （12）； in, Indicates the total number of time points. Indicates the first The time phase difference between the observed data and the simulated data corresponding to each seismic trace; Based on the difference in arrival time Dynamically select seismic traces to participate in loss calculation and construct first arrival constraints, i.e., set a time threshold. ,when If the seismic trace is considered to have a good match with the current initial model, then the seismic trace will be temporarily excluded from the loss calculation; otherwise, only those seismic traces that meet the criteria will be selected. The seismic traces are used in the inversion, and the expression is as follows: （13）； in, Represents the set of all seismic traces. This represents the selected set of seismic traces; The first-arrival constraint objective function construction module uses the first-arrival constraints to construct a dynamic objective function. Only seismic traces with first-arrival time differences less than a time threshold are selected for loss calculation. Therefore, the dynamic objective function based on the first-arrival constraints... The expression is as follows: （14）； in, This represents the selected set of seismic traces. Indicates the first The earthquake track in time Simulation data, This represents the corresponding observation data. Indicates the number of time steps.

Images