Robust full waveform inversion method of first arrival constraint physical embedding recurrent neural network

Abstract

The invention discloses a robust full-waveform inversion method of a first-arrival constraint physical embedded recurrent neural network, which comprises the following steps of: solving a wave equation through forward propagation of the physical embedded recurrent neural network, generating a forward shot gather record, screening out a seismic trace with high matching degree with observation data by using first-arrival information, and participating in loss calculation to obtain a seismic trace with high matching degree with observation data; and along with the proceeding of the inversion process, the range of seismic traces participating in the inversion is gradually expanded, it is ensured that the inversion result gradually tends to a globally optimal solution, back propagation is carried out by using space-time field gradient information recorded by forward propagation of the recurrent neural network, the speed parameter is further optimized, and the inversion precision is improved. According to the method, the physical process of seismic wave propagation is embedded into the recurrent neural network, the first arrival constraint and the progressive inversion strategy are introduced, the dependence of the FWI on the initial model is remarkably reduced, and the underground velocity field structure can still be accurately reconstructed even if the initial model and a real model have remarkable deviation.

Description

Technical Field

[0001] The invention belongs to the technical field of seismic data processing, and in particular relates to a robust full waveform inversion method of first-arrival constraint physical embedded recurrent neural network. Background Art

[0002] Full-Waveform Inversion (FWI) is an important geophysical method for obtaining underground medium parameter information. Its core goal is to invert underground medium parameters such as velocity and density by minimizing the difference between observed data and synthetic data. FWI relies on the numerical solution of the wave equation and can use the complete waveform information of seismic waves to provide high-resolution underground structure images. However, FWI is essentially a highly nonlinear optimization problem, and its inversion results are 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 prone to fall into a local minimum, causing the inversion results to deviate 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 a multi-scale strategy to reduce the impact 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 under complex geological conditions. The lack of low-frequency data will seriously affect the construction of the initial model. Secondly, although the multi-scale method can alleviate the nonlinear problem to a certain extent, its effect depends on the reasonable division of frequency bands, and the inversion process may still lead to a local optimal solution due to inaccurate initial models.

[0004] In recent years, the rise of deep learning has provided new ideas for solving nonlinear problems in FWI. Although existing data-driven methods can train models through a large amount of data, their results lack physical interpretability and have high requirements on data quality and scale. In order to overcome this problem, Physics-Informed Neural Networks (PINN) were introduced into FWI. By embedding physical constraints such as wave equations into neural networks, the physical interpretability of inversion results is enhanced. 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 results may still be unsatisfactory. Summary of the invention

[0005] To solve the above technical problems, the present invention provides a robust full waveform inversion method (FAFWI) with first arrival constraint physically embedded in a recurrent neural network. The physical process of seismic wave propagation is embedded in a recurrent neural network (RNN), synthetic seismic data is generated through physically driven forward simulation, and seismic traces are dynamically screened in combination with first arrival time difference constraints to gradually optimize the inversion results.

[0006] The technical solution adopted by the present invention is: a robust full waveform inversion method with first arrival constraint physics embedded in recurrent neural network, and the specific steps are as follows:

[0007] S1. Construct a physical embedding recurrent neural network (PIRNN) framework based on first arrival constraints;

[0008] The framework includes: 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.

[0009] The physical embedding recurrent neural network module includes: an initial model and an SGFD RNN operator module. The initial model is a speed model.

[0010] S2. Based on the framework constructed in step S1, the physically 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 uses the velocity model parameters as trainable parameters in the network to complete the seismic forward modeling based on physical constraints;

[0011] S3, extracting and utilizing the first arrival wave information from the observed data and the simulated data based on the first arrival information extraction module, constructing the first arrival constraint, and then constructing the dynamic objective function using the first arrival constraint by the first arrival constraint objective function construction module, and only selecting the seismic traces whose first arrival time difference is less than the time threshold to participate in the loss calculation;

[0012] Among them, constructing the first arrival constraint is to extract the first arrival time difference by calculating the cross-correlation function between the observed data and the simulated data, and dynamically select the 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 space-time field gradient information recorded by the forward propagation of the physically embedded recurrent neural network module is utilized, and the inversion optimization module updates the velocity model parameters through back propagation, gradually optimizes the inversion results, and realizes robust full waveform inversion.

[0014] Furthermore, the 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 pressure, p x and p z represent the pressure components in the x and z directions respectively, v represents the speed of sound, and v x and v z They represent the particle velocity in the x and z directions respectively, v(r) represents the sound wave velocity that needs to be inverted, and (r,t) represents the value of the physical quantity at position r at time t.

[0018] The velocity model is then embedded into the RNN network as a trainable parameter, the staggered grid high-order finite difference method is used to solve the wave equation, and the perfectly matched 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, the pressure p x , p z Defined at integer grid points, the particle velocity v x , v z Defined at half-grid points, for second-order time accuracy and second-order finite differences, Expanding at time and [ixΔt,izΔz], the expression is as follows:

[0020]

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

[0022] Further processing of the discrete equations of equations (2)-(5) yields the following equations about p: x , p z , v x , v z The recursive formula is as follows:

[0023]

[0024] Then, the wave field state at the next moment is solved in sequence according to the time stepping method, and the forward modeling of the seismic wave field can be performed through RNN.

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

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

[0027] In the inversion process, first for each seismic trace, the observed data is calculated With simulated data The cross-correlation function 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 is Corresponding to the best time alignment between the observed and simulated data, the expression is as follows:

[0030]

[0031] Then the first arrival time difference is extracted by the peak position of the cross-correlation function Extract the first arrival time difference, 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 sth seismic trace.

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

[0035]

[0036] in, represents the set of all seismic traces, Represents the selected seismic trace set.

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

[0038]

[0039] in, represents the set of selected seismic traces, represents the simulated data of the sth seismic trace at time t, Represents the corresponding observation data, and T represents the number of time steps.

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

[0041] In the inversion process, a progressive strategy is adopted to calculate the number of seismic traces that meet the threshold of the first arrival time difference. Dynamically adjust the inversion process.

[0042] Where I represents the indicator function, which is used to select the appropriate seismic trace, and Shots represents the total number of seismic traces.

[0043] Then iterative optimization is performed. In each iteration, steps S2-S3 are repeated, and the initial model is forward modeled through the SGFDRNN operator module to obtain the simulated data while recording the gradient information of the entire space-time field. The simulated data and the observed data are then screened for seismic traces according to the first arrival constraints. The losses of the screened seismic traces are calculated and back-propagated to update the initial model. This process is repeated until most of the seismic traces are selected for 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.

[0044] Among them, the inversion optimization module uses the formula Update speed model parameter v, γ represents the learning rate, k represents the number of iterations, Represents the gradient of the loss function with respect to the speed parameter.

[0045] The beneficial effects of the present invention are as follows: the method of the present invention first solves the wave equation by forward propagation of a physically embedded recurrent neural network to generate a forward shot gather record, and then uses the first arrival information to screen out seismic traces with a high degree of match with the observed data to participate in the loss calculation, thereby avoiding the waveform mismatch problem caused by an inaccurate initial model, and gradually expands the range of seismic traces participating in the inversion as the inversion process proceeds to ensure that the inversion result gradually approaches the global optimal solution, and uses the space-time field gradient information recorded by the forward propagation of the recurrent neural network for reverse propagation to further optimize the velocity parameters and improve the inversion accuracy. The method of the present 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 cycle jump problem caused by an inaccurate initial model in the existing full waveform inversion, and significantly reduces the dependence of FWI on the initial model. Even when there is a significant deviation between the initial model and the true model, the underground velocity field structure can still be accurately reconstructed. The first-arrival constraint strategy can screen out seismic traces with a high degree of matching with the initial model at the beginning of the 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 result gradually approaches the global optimum, and significantly improving the robustness of the inversion result. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] Figure 1 The present invention is a flowchart of a robust full waveform inversion method of first arrival constrained physics embedded in a recurrent neural network.

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

[0048] Figure 3 Schematic diagram of the real speed model and the constant speed model in an embodiment of the present invention.

[0049] Figure 4 Schematic diagram of inversion results in an embodiment of the present invention.

[0050] Figure 5 Schematic diagram of velocity profile in an embodiment of the present invention.

[0051] Figure 6 Schematic diagram of the loss function of the method of the present invention on the Marmousi model in an embodiment of the present invention, the total phase difference Δφ and the variation curve of the number of selected seismic traces with the number of iterations. DETAILED DESCRIPTION

[0052] The method of the present invention is further described below in conjunction with the accompanying drawings and embodiments.

[0053] like Figure 1 As shown, a flow chart of a robust full waveform inversion method of first arrival constraint physical embedded recurrent neural network of the present invention, the specific steps are as follows:

[0054] S1. Construct a physical embedding recurrent neural network (PIRNN) framework based on first arrival constraints;

[0055] The framework includes: 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 includes: an initial model and an SGFD RNN operator module. The initial model is a speed model.

[0057] S2. Based on the framework constructed in step S1, the physically 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 uses the velocity model parameters as trainable parameters in the network to complete the seismic forward modeling based on physical constraints;

[0058] S3, extracting and utilizing the first arrival wave information from the observed data and the simulated data based on the first arrival information extraction module, constructing the first arrival constraint, and then constructing the dynamic objective function using the first arrival constraint by the first arrival constraint objective function construction module, and only selecting the seismic traces whose first arrival time difference is less than the time threshold to participate in the loss calculation;

[0059] Among them, constructing the first arrival constraint is to extract the first arrival time difference by calculating the cross-correlation function between the observed data and the simulated data, and dynamically select the seismic traces participating in the 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 physically embedded recurrent neural network module is utilized, and the inversion optimization module updates the velocity model parameters through back propagation, gradually optimizes the inversion results, and realizes robust full waveform inversion.

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

[0062] First, the seismic forward modeling embedded recurrent neural network is explained. Establishing an accurate seismic forward modeling process is the basis of full waveform inversion. Compared with the existing full waveform inversion, the physical embedded 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] Where p represents pressure, px and p z represent the pressure components in the x and z directions respectively, v represents the speed of sound, and v x and v z They represent the particle velocity in the x and z directions respectively, v(r) represents the sound wave velocity that needs to be inverted, and (r,t) represents the value of the physical quantity at position r at time t.

[0065] like Figure 2 As shown in the figure, in order to ensure the high accuracy of the forward modeling process, the velocity model is embedded in the RNN network as a trainable parameter (the time dependence of the wave equation is mapped into the hierarchical structure of the recurrent neural network (RNN)), the staggered grid high-order finite difference method is used to solve the wave equation, and the perfectly matched layer (PML) boundary condition is introduced to effectively absorb the outgoing wave and reduce boundary reflections, that is, the forward modeling is performed step by step through the SGFD RNN operator module (Staggered Grid Finite DifferenceRecurrent Neural Network Operator).

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

[0067] In this process, underground medium parameters are embedded as trainable parameters into each SGFD RNN operator, 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 to the underground parameters through back propagation and iteratively update these parameters. In this way, this embodiment realizes seismic forward modeling based on physical constraints, laying the foundation for subsequent full waveform inversion.

[0068] First, discretize equation (1), that is, the pressure p x , p z Defined at integer grid points, the particle velocity v x , v z Defined at half-grid points, for second-order time accuracy and second-order finite differences, Expanding at time and [ixΔt,izΔz], the expression is as follows:

[0069]

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

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

[0072]

[0073]

[0074] Then, the wave field state at the next moment is solved in sequence according to the time stepping method, and the forward modeling of the seismic wave field can be performed through RNN.

[0075] Recurrent neural network (RNN) is a type of artificial neural network that can process time series data. 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 equations, 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 moment, which is highly consistent with the time stepping method in solving wave equations. Specifically, each layer of RNN is used to calculate and store the wave field information at a certain moment, including pressure and particle velocity, and each RNN unit contains a chain relationship about the velocity parameter v. During the forward modeling process, the neural network will record the chain relationship of the entire space-time field, and will calculate the gradient of the loss function about the velocity parameter v according to the automatic differentiation mechanism during back propagation, and update the velocity. This structure allows the time evolution process of the wave field to be mapped to the hierarchical structure of RNN. Each layer not only stores the wave field information at the current moment, but also receives the output of the previous layer as input, thereby well simulating the time dependence of the wave equation.

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

[0077] In view of the problem of waveform matching cycle jump in full waveform inversion (FWI), this embodiment introduces a first arrival constraint strategy to alleviate the nonlinear problem. The first arrival wave is the first wave to arrive at the detector during the propagation of seismic waves, and is usually located near the time point when the seismic trace has the strongest energy. Under the 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 invention 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 degree of matching between the observed 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 cycle, resulting in a cycle jump phenomenon in FWI during waveform matching, thereby causing the model update to fall into a local optimal solution. To avoid this problem, this embodiment proposes a first arrival constraint strategy, the core of which is to select only seismic traces with a first arrival time difference Δt less than half a cycle to participate in the construction and optimization of the objective function. Through this strategy, the cycle jump phenomenon can be effectively avoided and the model parameters can be guided to iteratively update in the correct direction.

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

[0079] In the inversion process, first for each seismic trace, the observed data is calculated With simulated data The cross-correlation function 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 is Corresponding to the best time alignment between the observed and simulated data, the expression is as follows:

[0082]

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

[0084]

[0085] 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 sth seismic trace.

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

[0087]

[0088] in, represents the set of all seismic traces, Represents the set of selected seismic traces. This dynamic selection mechanism ensures that in the early stage of inversion, when the initial model differs greatly from the true model, only seismic traces with smaller first arrival time differences participate in the loss calculation, thereby effectively reducing error accumulation and improving the stability of the inversion. With the iterative optimization of model parameters, the number of seismic traces that meet the first arrival time difference threshold condition will gradually increase, thereby providing richer underground medium structure information and eventually converging to the global optimal solution.

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

[0090]

[0091] in, represents the set of selected seismic traces, represents the simulated data of the sth seismic trace at time t, Represents the corresponding observation data, and T represents the number of time steps.

[0092] Through this dynamic objective function, this embodiment can gradually introduce more seismic traces during the inversion process, ensuring that the inversion results start from reliable data and gradually improve accuracy and comprehensiveness.

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

[0094] In the inversion process, a progressive strategy is adopted. In the initial stage, only seismic traces with first arrival time difference less than the time threshold are selected for inversion. As the inversion progresses, the number of seismic traces participating in the inversion is gradually increased to ensure that the inversion results gradually tend to the global optimum. That is, by Dynamically adjust the inversion process.

[0095] Where I represents the indicator function, which is used to select the appropriate seismic trace, and Shots represents the total number of seismic traces.

[0096] Then iterative optimization is performed. In each iteration, steps S2-S3 are repeated, and the initial model is forward modeled through the SGFDRNN operator module to obtain the simulated data while recording the gradient information of the entire space-time field. The simulated data and the observed data are then screened for seismic traces according to the first arrival constraints. The losses of the screened seismic traces are calculated and back-propagated to update the initial model. This process is repeated until most of the seismic traces are selected for 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.

[0097] Among them, the inversion optimization module uses the formula Update speed model parameter v, γ represents the learning rate, k represents the number of iterations, Represents the gradient of the loss function with respect to the speed parameter.

[0098] This embodiment is further verified by experiments. Figure 3 As shown in the figure, the difference between the original FWI and the method of the present invention in terms of initial model dependence is compared and analyzed in the Marmousi model. The constant velocity model is used as the initial model for full waveform inversion. This initial model is significantly different from the true velocity model and has no low-frequency information and logging information. Therefore, it 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, and the time threshold Δt threshold The learning rate is set to 2.5ms and is uniformly set to 40 and remains unchanged.

[0099] Figure 3 (a) is the real velocity model and observation system settings. The red five-pointed star indicates the location of the source, and the black dotted line indicates the location of the detector. Figure 3 (b) is the constant speed model with a speed value of 1800m / s. Figure 3 (b) shows the constant velocity initial model used in Marmousi inversion. The velocity value of each point is 1800m / s. If such a velocity model does not add any prior structural information as a constraint, it is easy to fall into the local minimum and obtain an erroneous physical solution under the original FWI.

[0100] Figure 4 is a schematic diagram of the inversion results, such as Figure 4 As shown in (b), the inversion results under the original FWI show high uncertainty and completely deviate from the true velocity structure. Compared with the original FWI, the inversion results obtained by the method proposed in this invention are as follows: Figure 4As shown in (a), it can be seen that under the same initial model, the FAFWI method of the present invention can still obtain relatively accurate inversion results within the effective range of the VSP observation system. It not only reconstructs the velocity structure of the shallow and middle layers, but also can better invert the velocity structure and change trend in the deep area.

[0101] like Figure 5 As shown in the figure, in order to quantitatively evaluate the velocity field characteristics in the vicinity of the wellbore, vertical velocity profiles at distances of ±75 m with the wellbore as the center were extracted. Figure 5 (a) is the velocity profile at 75 m from the left side of the wellbore under the method of the present invention (FAFWI), Figure 5 (b) is the velocity profile 75 m from the right side of the wellbore using the method of the present invention (FAFWI). The profile comparison analysis shows that the inversion results are highly consistent with the true model in most areas from shallow to deep layers.

[0102] From the spatial distribution characteristics of the velocity field inversion results, the effective inversion area presents an obvious inverted triangle shape. This is because the detectors in the VSP observation system are arranged vertically along the wellbore, and the seismic reflection waves in the lateral area far away from the gun-well pair are difficult to be effectively received, so the inversion results conform to the physical laws. On this basis, the R2 score, structural similarity index SSIM and normalized correlation coefficient of the inversion results within this range and the true model are calculated to comprehensively evaluate the inversion results obtained by the method of the present invention in terms of linear correlation, structural correlation, etc. The results obtained are shown in Table 1 (quantitative comparison of inversion results on the Marmousi model using the FAFWI method). The values ​​in Table 1 show that the model inverted by the method of the present invention presents a strong linear relationship and structural similarity with the true model, and also has advantages in normalized correlation.

[0103] Table 1

[0104]

[0105] like Figure 6 As shown in the changes of Marmousi's training indicators, as the inversion process continues, the velocity parameters will be updated in the correct reverse direction, so that |Δt is satisfied in each iteration s |<Δt threshold The number of seismic channels is increasing, and the total phase difference will gradually approach 0, and the number of selected seismic traces will gradually increase. This dynamic adjustment mechanism realizes the dynamic objective function. As the inversion proceeds, more seismic traces are allowed to participate in the inversion, making the inversion results more comprehensive and accurate. Figure 6 (a) is the dynamic loss curve. Figure 6(b) is a schematic diagram showing the changing trend of the total phase difference Δφ of all seismic traces. Figure 6 (c) To meet the time threshold Δt threshold Schematic diagram of the changing trend of the number of seismic traces.

[0106] Model tests and actual data tests prove that the FAFWI full waveform inversion network structure using the method of the present invention significantly reduces the dependence on the initial model while ensuring a certain inversion accuracy. This proves that the method of the present invention 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, numerical experiments have verified that the method of the present invention can still accurately reconstruct the underground velocity field structure when 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 the first arrival wave information, the seismic traces with a high degree of matching with the initial model are dynamically selected to participate in the inversion, which effectively reduces the error accumulation caused by the inaccurate initial model. At the same time, combined with the spatiotemporal field gradient information of the 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 appreciate that the embodiments described herein are intended to help readers understand the principles of the present invention, and should be understood that the scope of protection of the present invention is not limited to such specific statements and embodiments. For those skilled in the art, the present invention may have various changes and variations. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of the claims of the present invention.

Claims

1. A robust full waveform inversion method with first-arrival constraint physics embedded in a recurrent neural network, the specific steps are as follows: S1. Construct a physical embedding recurrent neural network (PIRNN) framework based on first arrival constraints; The framework includes: First arrival information extraction module, physical embedding recurrent neural network module, first arrival constraint objective function construction module, inversion optimization module; The physical embedding recurrent neural network module includes: an initial model and an SGFD RNN operator module; the initial model is a speed model; S2. Based on the framework constructed in step S1, the physically 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 uses the velocity model parameters as trainable parameters in the network to complete the seismic forward modeling based on physical constraints; S3, extracting and utilizing the first arrival wave information from the observed data and the simulated data based on the first arrival information extraction module, constructing the first arrival constraint, and then constructing the dynamic objective function using the first arrival constraint by the first arrival constraint objective function construction module, and only selecting the seismic traces whose first arrival time difference is less than the time threshold to participate in the loss calculation; Among them, constructing the first arrival constraint is to extract the first arrival time difference by calculating the cross-correlation function between the observed data and the simulated data, and dynamically select the seismic traces participating in the inversion based on the first arrival time difference; S4. Based on the dynamic objective function constructed in step S3, the space-time field gradient information recorded by the forward propagation of the physically embedded recurrent neural network module is utilized, and the inversion optimization module updates the velocity model parameters through back propagation, gradually optimizes the inversion results, and realizes robust full waveform inversion.

2. The robust full waveform inversion method of first arrival constraint physical embedded recurrent neural network according to claim 1 is characterized in that: The step S2 is specifically as follows: 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: Where p represents pressure, p x and p z represent the pressure components in the x and z directions respectively, v represents the speed of sound, and v x and v z They represent the particle velocity in the x and z directions respectively, v(r) represents the sound wave velocity to be inverted, and (r,t) represents the value of the physical quantity at position r at time t; Then, the velocity model is embedded into the RNN network as a trainable parameter, and the wave equation is solved by using the staggered grid high-order finite difference method. The perfectly matched layer boundary condition is introduced, that is, the forward calculation is performed step by step through the SGFD RNN operator module. Discretize equation (1), that is, the pressure p x , p z Defined at integer grid points, the particle velocity v x , v z Defined at half-grid points, for second-order time accuracy and second-order finite differences, Expanding at time and [ixΔx,izΔz], the expression is as follows: Where k is the index of the time step, ix and iz are the discrete number indexes of the spatial grid in the x and z directions respectively, and p k represents the discrete value of 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; Δt represents the discrete step length in time, that is, the interval between two adjacent time points in the grid; Further processing of the discrete equations of equations (2)-(5) yields the following equations about p: x , p z , v x , v z The recursive formula is as follows: Then, the wave field state at the next moment is solved in sequence according to the time stepping method, and the forward modeling of the seismic wave field can be performed through RNN; Among them, each layer of the RNN is used to calculate and store the wave field information at a certain moment, and each RNN unit contains a chain relationship about the speed parameter v. During the forward modeling process, the neural network will record the chain relationship of the entire space-time field, and will calculate the gradient of the loss function with respect to the speed parameter v according to the automatic differentiation mechanism during back propagation, and perform speed updates.

3. The robust full waveform inversion method of first arrival constraint physical embedded recurrent neural network according to claim 1 is characterized in that: The step S3 is specifically as follows: In the inversion process, first for each seismic trace, the observed data is calculated With simulated data The cross-correlation function The expression is as follows: Where τ represents the time delay, s represents the seismic trace, and t represents the time; the peak position of the cross-correlation function Corresponding to the best time alignment position between the observed data and the simulated data, the expression is as follows: Then the first arrival time difference is extracted by the peak position of the cross-correlation function Extract the first arrival time difference, the expression is as follows: Where T represents the total number of time points, Δt s Represents the time phase difference between the observed data and simulated data corresponding to the sth seismic trace; Based on the first arrival time difference Δt s , dynamically select the seismic traces involved in the loss calculation, and construct the first arrival constraint, that is, set a time threshold Δt threshold , when |Δt s |<Δt threshold , it is considered that the underground velocity information carried by the seismic trace matches the current initial model well; otherwise, the seismic trace will be temporarily excluded from the loss calculation; that is, only the seismic traces that satisfy |Δt s |<Δt threshold The seismic traces of are involved in the inversion, and the expression is as follows: in, represents the set of all seismic traces, Represents the selected seismic trace set; The first-arrival constraint objective function construction module uses the first-arrival constraint to construct a dynamic objective function, and only selects seismic traces with first-arrival time differences less than the time threshold to participate in the loss calculation. The expression of the dynamic objective function J(v) based on the first-arrival constraint is as follows: in, represents the set of selected seismic traces, represents the simulated data of the sth seismic trace at time t, Represents the corresponding observation data, and T represents the number of time steps.

4. The robust full waveform inversion method of first arrival constraint physical embedded recurrent neural network according to claim 1 is characterized in that: The step S4 is specifically as follows: In the inversion process, a progressive strategy is adopted to calculate the number of seismic traces that meet the threshold of the first arrival time difference. Dynamically adjust the inversion process; Where I represents the indicator function, which is used to select the appropriate seismic trace, and Shots represents the total number of seismic traces; Then, iterative optimization is performed. In each iteration, steps S2-S3 are repeated, and the initial model is forward-modeled through the SGFD RNN operator module to obtain the simulated data while recording the gradient information of the entire space-time field. Then, the simulated data and the observed data are screened for seismic traces according to the first arrival constraint, and the loss of the screened seismic traces is calculated and back-propagated to update the initial model. This process is repeated until most of the seismic traces are selected for 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. Among them, the inversion optimization module uses the formula Update speed model parameter v, γ represents the learning rate, k represents the number of iterations, Represents the gradient of the loss function with respect to the speed parameter.

Images