VSP sound wave full waveform inversion method of double-branch physical drive recurrent neural network

Abstract

The invention discloses a VSP sound wave full-waveform inversion method of a double-branch physical driving recurrent neural network, which is applied to the field of seismic data processing and aims to solve the problems of low interpretability, low calculation cost and the like due to the adoption of a complete data driving mode in the existing seismic velocity inversion technology. As the wave equation describes the change of the wave field in time and space, a time stepping method is adopted during numerical simulation, namely, the state of the wave field at the next moment is calculated according to the state at the previous moment. Taking a wave field state at a previous moment as a hidden layer of the RNN, performing finite difference solution on an acoustic wave equation, taking a speed parameter in the equation as a trainable parameter in a physically-driven RNN forward modeling network, and taking a shot gather record obtained by outputting through the physically-driven RNN forward modeling network each time as a predicted value; taking shot gather records obtained through physical driving RNN forward modeling output by setting real speed parameters as observation values; the process of solving the loss between the predicted value and the observed value to carry out back propagation to correct the speed parameter is called as a speed inversion process.

Description

Technical Field

[0001] The invention belongs to the field of seismic data processing, and particularly relates to a full waveform inversion technology. Background Art

[0002] Full Waveform Inversion (FWI) uses the maximum wavelength information contained in seismic data to reconstruct underground parameters by minimizing the error between observed data and predicted data. Full waveform inversion was first proposed by Lailly and Tarantola. The full waveform inversion problem is a strongly nonlinear problem, so FWI is highly dependent on the initial model. In addition, the data acquisition during seismic exploration lacks low-frequency components, which makes the inversion easy to fall into local minima and produce cycle jump problems. Warner and Guasch proposed adaptive waveform inversion (AWI) to reduce dependence on data and cycle jump problems, but it may reduce the resolution of the inversion results and increase the computational complexity. In general, the FWI method faces challenges such as how to balance computational efficiency, result accuracy, and dependence on the initial model.

[0003] Full waveform inversion is a strongly nonlinear method, and its results are multi-knotted. Therefore, it is highly dependent on the accuracy of the initial data. The key to traditional full waveform inversion is how to optimize the objective function, that is, the error between the observed data and the predicted data. The adjoint state method is generally used to calculate the gradient of the objective function relative to the velocity parameter, and then a suitable optimization algorithm is selected to update the velocity parameter. Common optimization algorithms include the steepest descent method, conjugate gradient method, Newton method, trust region method, etc.

[0004] With the improvement of large-scale computing power and the rapid development of deep learning technology, deep learning technology is increasingly widely used in the field of geophysical inversion. Deep learning has shown great potential in dealing with nonlinear problems and provides a new research direction for solving seismic velocity inversion problems. Traditional deep learning inversion methods generally adopt a completely data-driven approach, allowing neural networks to learn the mapping relationship between seismic observation data and geophysical models from a pre-constructed training data set. This approach makes deep learning easy to implement in the field of seismic exploration, but it faces problems such as poor interpretability and computational cost. In recent years, physics-guided neural networks (PINNs) have begun to be widely used. The forward module of wave equation physical modeling is integrated into the network loop, and the physically driven data residual is used as the loss function of full waveform inversion. Compared with the completely data-driven deep learning method, this method is more flexible and can merge multiple data types and structural information of any structure. It has stronger anti-noise ability, thereby improving the accuracy and reliability of the inversion results. However, since the PIRNN full waveform inversion method often consumes a lot of GPU computing resources, especially when dealing with large-scale three-dimensional problems. Summary of the invention

[0005] In order to solve the above technical problems, the present invention proposes a VSP acoustic wave full waveform inversion method based on a dual-branch physical-driven recurrent neural network, combines the physical-driven RNN full waveform inversion network with scale decomposition, and establishes a Dual-PIRNN (Dual Physical-Informed Recurrent Neural Network, based on a dual-branch physical-driven recurrent neural network), realizing a VSP acoustic wave full waveform inversion method with low GPU resource consumption based on machine learning and physics guidance.

[0006] The technical solution adopted by the present invention is: a VSP acoustic full waveform inversion method of a dual-branch physical-driven recurrent neural network, constructing a dual-branch physical-driven recurrent neural network, which specifically includes: an original-scale true velocity model forward modeling branch, an original-scale initial velocity model forward modeling branch, and a downscaled velocity model inversion branch; the original-scale true velocity model forward modeling branch, the original-scale initial velocity model forward modeling branch, and the downscaled velocity model inversion branch are implemented by a physically-driven RNN forward modeling network;

[0007] The acoustic wave equation is solved by finite difference, and the velocity parameters in the equation are used as trainable parameters in the physically driven RNN forward modeling network. The output of the forward modeling branch of the original-scale true velocity model is used as the observation value, and the output of the forward modeling branch of the original-scale initial velocity model is used as the prediction value. The loss between the prediction value and the observation value is input into the downscaled velocity model inversion branch to update the velocity parameters. After each iterative inversion, the velocity parameters on the downscaled velocity model inversion branch are updated to the forward modeling branch of the original-scale initial velocity model through upsampling for forward modeling again, and the update is cyclic.

[0008] Beneficial effects of the present invention: The Dual-PIRNN full waveform inversion network structure proposed in the present invention combines physical wave equation constraints to ensure the consistency of synthetic seismic data and observed seismic data, thereby improving the accuracy of the inversion results. The present invention introduces scale decomposition and constructs a dual-branch physical-driven recurrent neural network full waveform inversion architecture, which effectively reduces the GPU resource consumption in the inversion while ensuring the accuracy of the inversion results. The scale decomposition in the present invention can reduce the local minimum in the inversion and reduce the situation of falling into cycle jumps, thereby converging faster to obtain a velocity model. BRIEF DESCRIPTION OF THE DRAWINGS

[0009] Figure 1 It is the expanded form of the PIRNN architecture;

[0010] Figure 2 The internal structure of the staggered grid finite difference RNN operator;

[0011] Figure 3 This is the Dual-PRINN full waveform inversion flow chart;

[0012] Figure 4 are the correct model and the initial model;

[0013] Among them, (a) is the true velocity model, (b) is the initial model under Gaussian smoothing with variance sigma = 5, and (c) is the initial model under Gaussian smoothing with variance sigma = 10;

[0014] Figure 5 is the inversion result;

[0015] Among them, (a) is the velocity model obtained under the initial model of sigma = 5 and the time scale inversion, (b) is the velocity model obtained under the initial model of sigma = 5 and the time and space scale inversion, (c) is the velocity model obtained under the initial model of sigma = 10 and the time scale inversion, (d) is the velocity model obtained under the initial model of sigma = 10 and the time and space scale inversion;

[0016] Figure 6 is the velocity profile;

[0017] Among them, (a) is the velocity curve of the initial model with different smoothing degrees after time scale reduction and inversion at x=300m, (b) is the velocity curve of the initial model with different smoothing degrees after time scale reduction and inversion at x=700m, (c) is the velocity curve of the initial model with different smoothing degrees after time and space scale reduction and inversion at x=300m, (d) is the velocity curve of the initial model with different smoothing degrees after time and space scale reduction and inversion at x=700m. DETAILED DESCRIPTION

[0018] To facilitate those skilled in the art to understand the technical content of the present invention, the present invention is further explained below with reference to the accompanying drawings.

[0019] 1. First, explain the recurrent neural network:

[0020] Recurrent Neural Networks (RNN) is a type of artificial neural network. RNN can be divided into two categories: finite impulse response networks with directed acyclic graphs and infinite impulse response networks with directed cyclic graphs. Both types of networks can represent dynamic time behavior. Unlike deep neural networks, RNN can use its internal state (memory) to process input sequences so that the current output results are affected by past results. This makes them suitable for time-dependent signal processing tasks, such as text and audio processing. Since the wave equation describes the changes in the wave field in time and space, the time-stepping method is generally used in numerical simulation solutions, that is, the wave field state at the next moment is calculated based on the wave field state at the previous moment. The core feature of RNN is that it can process time series data, and its hidden state can carry the information of the previous moment to the next moment, which is highly consistent with the time-stepping solution method of the wave equation. Therefore, the seismic forward modeling calculation (time-stepping solution of the wave equation) can be integrated into the RNN unit through finite differences. The RNN network forward modeling architecture is as follows: Figure 1 As shown, it takes the wave field information of the previous moment as input, generates a wave field record at the current moment, and saves the wave field record generated at this moment in the GPU and passes it to the next moment to continue generating wave field records. Under the VSP Vertical Seismic Profiling (vertical seismic profiling) observation system, the space-time field information about the velocity model is obtained by collecting the wave field records at each time point.

[0021] 2. Then the staggered grid finite difference RNN operator network structure adopted by the present invention is described:

[0022] A high-order finite difference format with a staggered grid form is used to solve the first-order velocity-pressure acoustic wave equation. Existing studies have shown that the process of integrating seismic forward modeling into the forward process of the RNN network to achieve full waveform inversion is equivalent to the traditional full waveform inversion method. The velocity parameters are used as trainable parameters, and seismic forward modeling is performed based on the velocity parameters. The shot records at each time point are collected as predicted values. The error between the observed value and the predicted value is used for back propagation and updating the velocity parameters. This process is called the velocity inversion process. The automatic differentiation mechanism of the neural network is used, and other information can be easily added in the PIRNN network to improve the inversion accuracy and robustness. The internal structure of the staggered grid finite difference RNN operator is shown in Figure 2. Figure 2 shown.

[0023] Figure 2The operator demonstrated here is based on PIRNN (Physically-Informed Recurrent Neural Network) to achieve high-precision forward modeling. The input is velocity parameters, source location, and wave field information at the previous moment, and the output is the VSP shot gather record at the current moment. is the input at the current moment. In this structure, As the internal parameter of each SGFD RNN Operator (Staggered Grid Finite Difference Recurrent Neural Network Operator) unit, it is used to propagate between different RNN units, which contains the chain relationship related to the speed parameter v at each moment. The relationship in the entire time field will calculate the gradient related to v according to the automatic differentiation property in the back propagation, and then update the speed parameter through the Adam optimizer.

[0024] This embodiment takes the first-order pressure-velocity fluctuation equation as an example to construct the RNN recursive formula with second-order time accuracy in the time domain and second-order staggered grid finite difference accuracy:

[0025]

[0026] where p x and p z are the pressure components in the x and z directions, v x ,v z is the particle velocity in the x and z directions, and v is the sound velocity that needs to be inverted. Accurate forward modeling is the basis of full waveform inversion. The present invention will use high-order finite differences in a staggered grid format to solve the above-mentioned first-order pressure-velocity acoustic wave equation, and add SPML boundaries to ensure a high-precision forward modeling process.

[0027] Discretize the above equation and replace the pressure p x ,p z Translated to the integer grid points, the particle velocity v x , v z Defined on half grid points, Expanding at time and [ixΔt,izΔz], we get the following formula:

[0028]

[0029]

[0030] Among them, △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, which represents the interval between two adjacent time points in the grid;

[0031] Further processing of the above discrete equations gives x , p z ,v x , v z The recursive formula is:

[0032]

[0033] Combination Figure 2 It can be seen that in each time cycle, the present invention first uses the pressure p at time kΔt k and The particle velocity at time and Obtain The particle velocity at time and Then according to the obtained The particle velocity at time and and the pressure of kΔt and To calculate the pressure at time (k+1)Δt and In this process, each time step to solve the wave equation uses the previous moment kΔt and The wave field data at a certain moment, this process makes it possible to perform seismic forward modeling through RNN, where each layer of the RNN can be constructed to calculate and store the wave field information at a certain moment. This structure allows us to map the time evolution process of the wave field to the hierarchical structure of the RNN. Specifically, each layer of the RNN can represent the state of the wave field at a specific time step, including pressure and particle velocity. Under this framework, each layer of the RNN not only stores the wave field information at the current moment, but also receives the output of the previous layer as input. This structure simulates the time dependence of the wave equation well, where the solution at each moment depends on the state at the previous moment. The present invention sets the above-mentioned forward solution process as the SGFD RNNOperator operator for forward solution at each moment.

[0034] 3. Explain the dual-branch physical-driven recurrent neural network based on the physical embedded recurrent neural network:

[0035] The VSP acoustic full waveform inversion structure of the dual-branch physics-driven recurrent neural network established based on the physics-embedded recurrent neural network is shown in the figure. Figure 3First, we will have a real original scale velocity model, and use the SGFD RNN Operator to forward model the VSP observation shot record d obs In each iterative inversion process, the present invention first performs forward modeling on the original scale initial model through the SGFD RNN Operator to obtain the simulated shot gathering data d cal , through the L1 norm of d obs and d cal Find the loss and get the loss on the original scale Loss original The loss is used for inversion on the downscaled branch model through scale conversion to update the velocity parameters. After each iterative inversion, the velocity parameters on the downscaled model are updated to the original scale model through upsampling for forward modeling again, and the update is cyclical.

[0036] Compared with the completely data-driven deep learning method, the PIRNN-based full waveform inversion method is more flexible, can incorporate multiple data types and structural information of any structure, and has stronger noise resistance, thereby improving the accuracy and reliability of the inversion results. Since the PIRNN full waveform inversion method often consumes a lot of GPU computing resources, especially when dealing with large-scale three-dimensional problems.

[0037] The Dual-PIRNN full waveform inversion network structure proposed in the present invention combines physical wave equation constraints to ensure the consistency of synthetic seismic data and observed seismic data, thereby improving the accuracy of the inversion results. The present invention introduces scale decomposition and constructs a dual-branch physical-driven recurrent neural network full waveform inversion architecture, which effectively reduces the GPU resource consumption in the inversion while ensuring the accuracy of the inversion results. The scale decomposition in the present invention can reduce the local minimum in the inversion and reduce the situation of falling into cycle jumps, thereby converging faster to obtain a velocity model.

[0038] The implementation process of the present invention includes the following contents:

[0039] 1. Staggered grid finite difference acoustic wave equation embedded in recurrent neural network

[0040] In the time domain, the first-order velocity-pressure form of the acoustic wave equation is:

[0041]

[0042] Where p is the pressure, p x and p z represent the pressure in the x and z directions respectively, v x ,v z are the x and z components of the particle velocity, v is the speed of sound, and (r, t) represents the value of the physical quantity at position r at time t.

[0043] In order to make the network suitable for velocity model inversion, the RNN structure contains trainable parameters, such as Figure 2 As shown, that is, the velocity parameter v. Then, the seismic shot collection record (observation record) at the wave field detector simulated by the network at time t obtained according to the set correct velocity model is used as a label; the initial velocity model is set with reference to the correct velocity model, and can be smoothed or overall offset, and the seismic shot collection obtained by inputting the initial velocity model into the network is used as the predicted seismic record. Information is passed Figure 1 The trainable RNN network structure shown completes forward propagation and simulates the forward modeling of the discrete wave field of the acoustic wave model.

[0044] The standard network training problem is to minimize the error between the network output (predicted seismic records) and the labeled data (observed seismic records). The trainable network weights are discrete, they are the sampled values ​​of the velocity field. The inversion process of Dual-PIRNN can be expressed as an optimization problem that minimizes the difference between the observed data and the predicted data. The L2 norm is chosen to measure the data difference:

[0045]

[0046] In the formula, ‖·‖ 2 represents the L2 norm, Ns is the number of earthquake sources, and T is the total time. cal and d obs Represents predicted data and observed data respectively. In the deep learning framework, the market relationship between variables is automatically recorded, and the gradient is calculated through automatic differentiation during the back propagation process, and then the parameters are updated through the optimizer, such as the Adam optimizer.

[0047] 2. Dual-PIRNN network structure

[0048] By converting the loss obtained by forward modeling under the original scale branch to the downscale branch for inversion according to the scale, the high resource consumption in PIRNN inversion can be effectively reduced. Specifically, the difference between the original scale model and the downscale model lies in the different model parameter settings, such as the time sampling rate dt of the source function, the number of time points NT of the numerical simulation, and the multiple relationship in the number of grids. In each iterative inversion process, such as Figure 2 As shown in the figure, firstly, the SGFDRNN Operator is used to perform forward modeling on the original scale initial model to obtain the simulated shot gathering data d cal , through the L1 norm of d obs and d cal Find the loss and get the loss on the original scale Loss original :

[0049] Loss original =dobs -d cal

[0050] Here, the present invention obtains d by the VSP observation system cal and d obs is a time series matrix of shape NTxnz, so the calculated Loss original It is also a time series matrix with a shape of NT×nz. Then, the present invention uses a downsampling method based on average pooling to reduce the loss original Extract the loss Loss that can match the downscaling model dowm Finally, the VSP shot collection record d on the downscaled model is obtained by forward modeling on the other downscaled model through the SGFD RNN Operator. Dcal . down and d Dcal Add together to get the observation record d on the downscaled size model Dobs :

[0051] Loss down = ExtractionOperator(Loss original )

[0052] d Dobs =d Dcal +Loss down

[0053] Get the downscaled observation record d Dobs Afterwards, the inversion process based on PIRNN can be performed on the downscaled model to update the parameters of the downscaled velocity model. The loss function used in the present invention is the mean square error function MSE, and the Adam optimizer is used for parameter update. The loss Loss on the downscaled model is:

[0054] Loss = MSE(d Dobs ,d Dcal )

[0055] The calculated loss is back-propagated on the downscaled model to update its velocity parameter V, and finally the updated velocity parameter V is converted to the original scale initial model through the upscaling operation. The above operation is the process of the next iterative inversion experience of the Dual-PIRNN network architecture. In this process, since the inversion operation of the present invention is performed on the downscaled model, the gradient that needs to be saved in the back propagation will be less than the original scale, thereby reducing the computational cost.

[0056] 3. Scale decomposition

[0057] The purpose of scale conversion is to enable the information obtained by forward modeling on the original scale model to be used on the downscaled model. This embodiment uses the Marmousi velocity model as an example and divides it into two different scales, a time scale down model and a space-time scale down model. The parameter setting example is shown in Table 1. Figure 5 As shown in the figure, the two scale models are called OriginalScale model and DownScale model. In this example, the grid points of the original scale model are 100x100, the grid spacing dx=dz=10, and the area size of 1000m x 1000m is simulated. The source function used is the sampling rate f s =4000Hz or f s = 8000 Hz Ricker wavelet, time point NT = 1600, the grid spacing dx, dz of the downscaled model is half of that of the original scale model, so the time point NT required for forward modeling is also half of the original scale, and the sampling rate of the source function used is f s =2000Hz, and the other parameter settings are shown in Table 1.

[0058] Table 1. Marmousi model parameter settings at different scales

[0059]

[0060] It should be noted that in the above embodiment, the original scale source sampling rate corresponding to the time scale reduction is 4000Hz, and the original scale source sampling rate corresponding to the time and space scale reduction is 8000Hz. In actual applications, the parameter settings of time scale reduction and time and space scale reduction are determined according to how much resource consumption needs to be reduced.

[0061] The technical effects of the present invention are described below in conjunction with specific data:

[0062] The Marmousi true velocity model and initial velocity model are established. The true velocity model is obtained by resampling the Marmousi public data set, such as Figure 4 As shown in (a), the velocity range is [1500,5500]. The initial velocity model is obtained by smoothing the real Marmousi model with different variances. Figure 4 (b) and Figure 4 (c) The initial model settings when the variance σ = 5 and σ = 10 are respectively, and the synthetic data set is calculated by forward propagation using the VSP observation system. Two downscaling strategies are used for inversion, namely downscaling in time and downscaling in space and time. The parameter settings are shown in Table 1, and Adam is used as the optimization algorithm. The inversion results are shown in Figure 5 As shown by Figure 5It can be seen that the inversion results under the reduced time scale can accurately identify the shallow velocity interface and also have a high resolution in the deep area. Figure 5 (a) and (c) show that good inversion results can be obtained under different smoothing initial models, which reflects the robustness of the method. The inversion results of reducing the spatiotemporal scale can accurately predict the velocity boundary of the shallow area and the velocity structure and change trend of the deep area. Overall, due to the loss of more information in spatiotemporal decomposition, the effect is not as good as reducing the time scale, but it saves more GPU resources.

[0063] The distances of 300m and 700m were selected as slices for speed curve comparison. The comparison results are shown in Figure 6 As shown:

[0064] Figure 6 The three curves are given in Figure 1: True value, result value under the initial model of σ = 10, and result value under the initial model of σ = 5. Under the time-reducing inversion, Figure 6 (a) and (b) show that in the high wrinkle area, the velocity values ​​in the shallow area at x = 300m and x = 700m are very close to the true value, and at x = 700m, the two are almost consistent. Under the inversion of down-scale time and space, Figure 6 (c) and (d) show that the velocity values ​​and changing trends can be accurately obtained in the shallow areas of x = 300m and x = 700m. There are periodic jumps in the deep areas, but the overall velocity trend is correct. The conclusion is consistent with the previous discussion. Figure 6 Here, DownTime corresponds to the down time, and DownTimeSpace corresponds to the down time and space.

[0065] Table 2 shows the GPU resource usage under different inversion strategies. Under the PIRNN full waveform inversion structure, most of the memory resource usage comes from saving the intermediate wave field data related to the velocity parameter v. We use the Dual-PIRNN architecture to reduce the numerical simulation time points and grid points according to different downscaling strategies. As shown in Table 2, under the time-scale downscaling inversion strategy, the intermediate activation value usage is reduced from 32Gb to 16Gb, reducing the resource usage by about 50%. Under the time-space downscaling, the intermediate activation value usage is reduced from 32Gb to about 2Gb.

[0066] Table 2 GPU resource usage under different inversion strategies

[0067]

[0068] Model tests and actual data tests prove that the Dual-PIRNN full waveform inversion network structure proposed in the present invention significantly reduces the GPU resources required in the inversion process while ensuring a certain inversion accuracy. This proves that the method effectively combines physical constraints with neural networks and scale decomposition, reduces computational costs and ensures full waveform inversion accuracy.

[0069] 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 VSP acoustic full waveform inversion method using a dual-branch physical driven recurrent neural network, characterized in that: Construct a dual-branch physics-driven recurrent neural network, which specifically includes: original scale true velocity model forward model branch, original scale initial velocity model forward model branch, downscale velocity model inversion branch; original scale true velocity model forward model branch, original scale initial velocity model forward model branch, downscale velocity model inversion branch are each implemented by a physics-driven RNN forward modeling network; The acoustic wave equation is solved by finite difference, and the velocity parameters in the equation are used as trainable parameters in the physics-driven RNN forward modeling network. The output of the original-scale true velocity model forward modeling branch is used as the observation value, and the output of the original-scale initial velocity model forward modeling branch is used as the predicted value. The loss between the predicted value and the observed value is input into the downscaled velocity model inversion branch to update the velocity parameters. After each iterative inversion, the velocity parameters on the downscaled velocity model inversion branch are updated to the original-scale initial velocity model forward model branch through upsampling for forward modeling again, and the update is cyclic.

2. The VSP acoustic full waveform inversion method of a dual-branch physical driven recurrent neural network according to claim 1 is characterized in that: The input of the physics-driven RNN forward modeling network is the velocity parameters, source location, and wave field information at the previous moment, and the output is the VSP shot gather record at the current moment.

3. The VSP acoustic full waveform inversion method of a dual-branch physical driven recurrent neural network according to claim 2 is characterized in that: The loss between the predicted value and the observed value is recorded as Loss original ; The Loss is reduced by downsampling method based on average pooling. original Extract and obtain the loss Loss that matches the downscaled velocity model down ; In the downscaled velocity model inversion branch, the VSP shot gather records d on the downscaled velocity model are obtained by forward modeling based on the physically driven RNN forward modeling network. Dcal ; Loss down and d Dcal Add together to get the observation record d on the downscaled velocity model Dobs .

4. The VSP acoustic full waveform inversion method of a dual-branch physical driven recurrent neural network according to claim 3 is characterized in that: Loss original The loss is obtained by calculating the L1 norm between the predicted value and the observed value.

5. The VSP acoustic full waveform inversion method of a dual-branch physical driven recurrent neural network according to claim 4 is characterized in that: The loss function used in the downscaling velocity model is the mean square error function.

6. The VSP acoustic full waveform inversion method of a dual-branch physical driven recurrent neural network according to claim 5 is characterized in that: The downscaling velocity model is specifically a time-downscaling model or a space-time-downscaling model.

Images