Seismic waveform inversion method and system based on deep learning image registration

Abstract

The invention relates to the technical field of seismic data processing, in particular to a seismic waveform inversion method and system based on deep learning image registration, and the method comprises the steps: carrying out the least square migration based on small offset observation data and a background velocity field, and constructing a reflection coefficient model; performing Born forward modeling according to the reflection coefficient and the background velocity field to obtain a simulated seismic record image; acquiring a seismic record deformation field by adopting deep learning image registration based on the simulated seismic record image and the seismic observation image; calculating a normalized deformation field, and judging whether the normalized deformation field meets a preset condition or not; if yes, terminating the iteration process; otherwise, executing the next step; and according to the seismic record deformation field and the accompanying seismic source, updating the background velocity field through a conjugate gradient algorithm, and repeating the steps. According to the method, high efficiency and high precision of deep learning are utilized, and the problem that a cycle jump problem which is prone to occurring in waveform matching falls into a local minimum value is avoided.

Description

Technical Field

[0001] The present invention relates to the technical field of seismic data processing, and specifically to a seismic waveform inversion method and system based on deep learning image registration. Background Art

[0002] Seismic exploration refers to a geophysical exploration method that uses elastic differences in underground media caused by artificially excited elastic waves to observe and analyze the propagation laws of seismic waves generated by artificial earthquakes underground to infer the properties and morphology of underground rock formations. In the data processing of seismic exploration, velocity modeling is a key step, and its accuracy directly affects the quality of seismic wave imaging. Full waveform inversion (FWI) is a nonlinear parameter inversion method that uses all the information of seismic waves (including travel time, amplitude, multiple waves, scattered waves, etc.) to synthesize data and fit observed data to estimate underground velocity density parameters. This method can obtain low, medium, and high full wave number velocity information, and the inversion results can accurately characterize and describe the complex underground structural morphology.

[0003] However, the traditional full waveform inversion method uses the least squares norm (L2) to measure the difference between the observed data and the calculated data. This measurement method will cause cycle jumping when there is a large difference between the simulated data and the observed data, especially when the simulated data moves more than half a cycle relative to the observed data. This requires FWI to have a good enough initial model. If the model can predict the data within half a cycle, FWI converges to the correct solution. Otherwise, FWI will be affected by cycle jumping and be guided to a local minimum.

[0004]

[0005] Among them, (t) represents the time direction; (h) represents the spatial direction; (t, h) represents that the data is a two-dimensional matrix; d cal (t, h) represents simulated data, which is seismic data obtained by computer numerical simulation; d obs (t,h) represents the observed data, which are the real seismic data collected in the field by the geophysical exploration ship; Indicates the data difference, d cal (t, h) and d obs The squared difference of (t,h).

[0006] In order to avoid the inversion falling into the local minimum, some scholars have proposed methods such as cross-correlation, dynamic time warping, and optimal transport. To calculate the one-dimensional travel time difference of the wave field, the cross-correlation function is often used to calculate the time shift as the optimization target, that is:

[0007]

[0008] in, Δτ represents the arrival time difference between the observed data and the simulated data as the waveform time of the image sequence; s represents the data index of an artificial source (shot point); r represents the data index of a seismic detector; J t It represents the sum of the squares of the arrival time differences of the waveforms of all seismic traces in the work area.

[0009] At a higher level, the process of finding differences and optimizing seismic records can be solved using image registration methods. The registration problem is simply to align one image to another so that the aligned images are as similar as possible. That is, given a floating image I m and a fixed image I f , the predicted deformation field Φ(t, h), I w represents the floating image after the floating image is restored by the deformation field Φ(t, h).

[0010] I w =I m *Φ(t,h) (3)

[0011] The deformation field Φ(t, h) represents the position movement of the sampling point in the t direction and the h direction.

[0012] Registration methods can be classified into traditional methods and machine learning-based methods. Traditional image registration methods use local image similarity optimization methods or feature registration methods. The former uses continuous iteration to find the optimal solution for spatial transformation, and has better results when dealing with image registration problems with small deformations; the latter mostly uses feature points for matching, which performs better when the image is severely deformed, but the computational cost is high. Traditional methods require manual repeated adjustment of parameters and are not suitable for automated large-scale data operations. Machine learning-based registration methods refer to registration methods that use neural network training. In this method, parameters are shared. First, a model is trained using a large amount of data, and then the trained model is used to register new images. Although the training process of this method is computationally expensive, the registration process is an order of magnitude faster than traditional methods.

[0013] Currently, full waveform inversion FWI technology faces several key technical problems in practice:

[0014] (1) When low-frequency data is missing and the L2 objective function is used to calculate the data residual, cycle jumps (waveform misalignment) are likely to occur. Velocity sections or diffraction points appear underground, and reflected waves may have complex cross-overlapping, which further aggravates the cycle jump problem.

[0015] (2) The amplitude of seismic records is affected by velocity, Lamé constant, and density coupling, and has a high degree of uncertainty. However, due to the limitations of limited observation aperture and acoustic wave equation simulation, only the phase accuracy of the wave field can be guaranteed. On this basis, it is difficult to obtain an accurate velocity structure based on the L2 objective function due to the influence of amplitude.

[0016] (3) Methods for calculating phase difference include cross-correlation (CC), dynamic time warping (DTW), and optimal transport method (0T). If cross-correlation is used to calculate phase difference, it is necessary to divide the time window and pick the maximum time movement difference. If the analysis time window is too large, the resolution of the cross-correlation phase difference result will be very low. If the analysis time window is too small, the stability will be poor. The dynamic time warping (DTW) algorithm can only calculate one-dimensional differences in the time direction. When the shallow velocity changes laterally, the apex of the diffraction wave hyperbola moves in time and space between the simulated data and the observed data. The one-dimensional method cannot describe this lateral movement. The optimal transport objective function has a high computational cost (0T). The calculation time of the transport scheme for a two-dimensional image may exceed the forward shot gathering time. Summary of the invention

[0017] The purpose of the present invention is to provide a seismic waveform inversion method and system based on deep learning image registration in order to solve at least one of the above technical problems.

[0018] The present invention provides a seismic waveform inversion method based on deep learning image registration, comprising the following steps:

[0019] Step S102, performing least square migration based on small offset observation data and background velocity field to construct a reflection coefficient model, and performing Born forward modeling based on the reflection coefficient and the background velocity field to obtain a simulated seismic recording image;

[0020] Step S104, based on the simulated earthquake record image and the earthquake observation image, using deep learning image registration to obtain the earthquake record deformation field;

[0021] Step S106, calculate the normalized deformation field using the following formula:

[0022] Φ_nor(t,h)=Φ(t,h) / wavelet_length(t,h) (4)

[0023] Among them, Φ_nor(t, h) represents the normalized deformation field; Φ(t, h) represents the earthquake record deformation field; wavelet_length(t, h) represents the wavelength of the main frequency wave of the current data, (t) represents the time direction, and (h) represents the space direction;

[0024] Determine whether the normalized deformation field meets the preset condition; if so, the iteration process is terminated; otherwise, execute step S108;

[0025] Step S108, calculating the accompanying seismic source according to the earthquake record deformation field, loading the accompanying seismic source into the wave equation for wave field propagation, using cross-correlation to calculate the low wave number gradient component to update the background velocity field, and repeating steps S102 to S106.

[0026] Furthermore, in step S102, the least squares shift is implemented by using a shift deconvolution method.

[0027] Furthermore, in step S104, the simulated earthquake record image is used as a floating image, and the earthquake observation image is used as a fixed image, and deep learning image registration is performed to obtain the earthquake record deformation field.

[0028] Further, step S104 includes the following steps:

[0029] Input an image to be registered, and output a deformation vector field after learning with a multi-scale parallel fully convolutional neural network; wherein the image to be registered includes the floating image and the fixed image;

[0030] The floating image generates a deformable image based on the deformation vector field using a spatial transformation network;

[0031] The difference between the fixed image and the deformed image is calculated according to the loss function, and the parameters of the seismic record deformation field are updated, and multiple iterations are performed to obtain the final seismic record deformation field.

[0032] Further, a three-layer multi-scale parallel downsampling module is used to reduce the image size of the image to be registered and obtain a multi-scale downsampling feature map;

[0033] The multi-scale down-sampled feature map is sequentially passed through a pyramid dilated convolution module and an adaptive channel attention module for feature extraction and feature processing;

[0034] The feature map processed by the adaptive channel attention module is input into the direction perception layer and the upsampling layer to obtain the deformation vector field of the input image pair.

[0035] Furthermore, the loss function adopts the following formula:

[0036] L(I f , I m ,Φ)=L(I f , I m ,Φ)+λL smoth (I f , I m , Φ) (5)

[0037] Among them, I f Represents a fixed image; I mrepresents floating image; Φ represents earthquake record deformation field Φ(t,h); L(I f , I m , Φ) represents I m The grayscale similarity of the transformed image; λ represents the spatial smoothness coefficient; L smoth (If,I m , Φ) indicates smoothing of image grayscale similarity.

[0038] Further, in step S106:

[0039] The preset condition is: the normalized deformation field Φ _ nor(t,h)<1%.

[0040] Further, step S108 includes the following steps:

[0041] The accompanying earthquake source is calculated according to the following formula:

[0042]

[0043] Where, (t) represents the time direction; (h) represents the spatial direction; Δτ represents the arrival time difference between the observed data and the simulated data as the waveform time of the image sequence; Δs represents the arrival time difference between the observed data and the simulated data as a waveform space of the image sequence; gt(t+Δτ,h) represents the gradient used for velocity update; represents the first-order time derivative of the observed data; represents the second-order time derivative of the observed data; d cal Represents simulation data; W is the proportionality coefficient, which is used to balance the proportional relationship between time and space; represents the time-delayed accompanying source; It indicates that the accompanying source of the earthquake record causes the vertices of the diffraction waves to move laterally due to the lateral change of velocity;

[0044] The companion source is loaded into the wave equation for wave field propagation, and the low wave number gradient component is calculated by cross-correlation to update the background velocity field, specifically through the following formula:

[0045] The gradient is:

[0046]

[0047] Where, f(t, h) represents the source wave field; g t (t+Δτ, h) represents the accompanying source wave field; represents the gradient, represents the cross-correlation of the two at each time in the underground; J represents the objective function, and m represents the particle velocity.

[0048] The present invention also provides a seismic waveform inversion system based on deep learning image registration, comprising:

[0049] A simulated seismic record image acquisition module is used to perform least square migration based on small offset observation data and background velocity field to build a reflection coefficient model, and to perform Born forward modeling based on the reflection coefficient and the background velocity field to acquire a simulated seismic record image;

[0050] A seismic record deformation field acquisition module, used to acquire the seismic record deformation field by deep learning image registration based on the simulated seismic record image and the seismic observation image;

[0051] The seismic waveform inversion module is used to calculate the normalized deformation field using the following formula,

[0052] Φ_nor(t,h)=Φ(t,h) / wavelet_length(t,h) (4)

[0053] Among them, Φ_nor(t, h) represents the normalized deformation field; Φ(t, h) represents the earthquake record deformation field; wavelet_length(t, h) represents the wavelength of the main frequency wave of the current data, (t) represents the time direction, and (h) represents the space direction;

[0054] The seismic waveform inversion module is also used to determine whether the normalized deformation field meets the preset conditions; if so, the iteration process is terminated; otherwise, the background velocity field is updated through the velocity update module;

[0055] The velocity update module is used to calculate the accompanying earthquake source according to the earthquake record deformation field, load the accompanying earthquake source into the wave equation for wave field propagation, and calculate the low wave number gradient component by cross-correlation to update the background velocity field.

[0056] The present invention also provides an electronic device, comprising: a processor, a memory, and a computer program stored in the memory and executable on the processor, wherein when the computer program is executed by the processor, any seismic waveform inversion method based on deep learning image registration is implemented.

[0057] The beneficial effects of the present invention are:

[0058] (1) For the local extrema problem caused by frequency jumps in classical full waveform inversion, the present invention performs deep learning global matching to improve the dependence on the initial velocity model and adopts a deformation field to avoid the uncertainty problem caused by matching amplitude.

[0059] (2) For the calculation of deformation vector fields, the present invention adopts a multi-scale parallel downsampling model and an adaptive channel attention module, which is not only suitable for the calculation of large deformation fields but also has high registration accuracy for small deformations. After training, the efficiency of the present invention is significantly improved compared with traditional image registration methods.

[0060] (3) The deformation vector field of the present invention is a two-dimensional matrix, which is more reasonable than the one-dimensional time delay obtained by the current methods such as cross-correlation, dynamic time warping, and transport, and is easy to obtain a more accurate velocity model.

[0061] (4) The present invention uses Born forward modeling to obtain reflected wave data for matching, which is not only suitable for transmitted waves but also for reflected waves, thus avoiding the problem of insufficient speed update depth of classical full waveform inversion.

[0062] The present invention utilizes the high efficiency and high precision of deep learning to avoid the problem of frequency jump that waveform matching is prone to fall into the local minimum. The present invention can construct a more stable velocity model than conventional methods, thereby significantly improving the quality of seismic migration imaging. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 The process of the seismic waveform inversion method based on deep learning image registration according to one embodiment of the present invention is as follows Figure 1 ;

[0064] Figure 2 The process of the seismic waveform inversion method based on deep learning image registration according to one embodiment of the present invention is as follows Figure 2 ;

[0065] Figure 3 A comparison diagram of non-rigid image registration shot gather test according to an embodiment of the present invention;

[0066] Figure 4 A network structure diagram of a multi-scale downsampling module according to an embodiment of the present invention;

[0067] Figure 5 A schematic diagram of a multi-scale parallel fully convolutional neural network according to an embodiment of the present invention;

[0068] Figure 6 An initial model observation data shot gather according to one embodiment of the present invention;

[0069] Figure 7 A de-migration simulation shot gather according to an embodiment of the present invention;

[0070] Figure 8 The initial model data shot gather deformation field of one embodiment of the present invention;

[0071] Fig. 9An initial model data shot gather accompanied by a seismic source in one embodiment of the present invention;

[0072] Fig.10 An initial model velocity update gradient according to an embodiment of the present invention;

[0073] Fig.11 : Fig.11 (a) is a true velocity profile of a two-dimensional depression model; Fig.11 (b) is the initial velocity model; Fig.11 (c) is the velocity model updated by model inversion;

[0074] Fig.12 : Fig.12 (a) is the migration result of the true velocity model; Fig.12 (b) is the migration result of the initial velocity model; Fig.12 (c) is the migration result of the inverted velocity model;

[0075] Fig.13 This is a structural diagram of a seismic waveform inversion system based on deep learning image registration according to an embodiment of the present invention. DETAILED DESCRIPTION

[0076] The content of the present invention will now be discussed with reference to exemplary embodiments. It should be understood that the embodiments discussed are only to enable those skilled in the art to better understand and thus implement the content of the present invention, rather than implying any limitation on the scope of the present invention.

[0077] As used herein, the term “including” and variations thereof are to be interpreted as open-ended terms meaning “including, but not limited to.” The term “based on” is to be interpreted as “based, at least in part, on.” The terms “one embodiment” and “an embodiment” are to be interpreted as “at least one embodiment.”

[0078] Embodiment 1

[0079] Figure 1 The process of the seismic waveform inversion method based on deep learning image registration according to one embodiment of the present invention is as follows Figure 1 .like Figure 1 As shown, according to one embodiment of the present invention, a seismic waveform inversion method based on deep learning image registration includes the following steps:

[0080] Step S102, performing least square migration based on small offset observation data and background velocity field to construct a reflection coefficient model, and performing Born forward modeling based on the reflection coefficient and background velocity field to obtain a simulated seismic recording image;

[0081] Step S104, based on the simulated earthquake record image and the earthquake observation image, deep learning image registration is used to obtain the earthquake record deformation field;

[0082] Step S106, calculate the normalized deformation field using the following formula:

[0083] Φ -nor(t, h )= Φ (t, h ) / wave l et_ l engt h (t, h )(4)

[0084] Among them, Φ _ nor(t, h) represents the normalized deformation field; Φ(t, h) represents the earthquake record deformation field; wavelet_length(t, h) represents the wavelength of the main frequency wave of the current data, (t) represents the time direction, and (h) represents the space direction;

[0085] Determine whether the normalized deformation field meets the preset conditions; if so, the iteration process is terminated; otherwise, execute step S108;

[0086] Step S108, calculate the accompanying earthquake source according to the earthquake record deformation field, load the accompanying earthquake source into the wave equation for wave field propagation, use cross-correlation to calculate the low wave number gradient component to update the background velocity field, and repeat steps S102 to S106.

[0087] In this embodiment, the seismic waveform inversion method is combined with deep learning image registration technology and traditional seismic data processing technology.

[0088] First, using small offset observation data and background velocity field, least squares migration and Born forward modeling are used to simulate seismic recording images, providing basic data for subsequent image registration. Preferably, in this embodiment, least squares migration is implemented by using a migration deconvolution method. The migration deconvolution method constructs a frequency domain wave field propagation operator in the least squares migration imaging method, and constructs mutually conjugate forward modeling operators and migration operators according to the deconvolution imaging conditions, so that the imaging results are more accurate. This method can effectively eliminate the crosstalk noise in the multiple wave imaging results and improve the resolution and signal-to-noise ratio of the imaging results.

[0089] Then, through deep learning image registration technology, the method can accurately obtain the deformation field of the earthquake record, which is a key step in realizing seismic waveform inversion. Preferably, in this embodiment, the simulated earthquake record image is used as a floating image, and the observation image (seismic observation image) collected by the earthquake is used as a fixed image, and deep learning image registration is performed to obtain the earthquake record deformation field.

[0090] After obtaining the deformation field of the seismic record, the normalized deformation field is further obtained through normalization processing, and the preset conditions are used to determine whether the iterative process is terminated, which reflects the flexibility and adaptability of the algorithm. If the preset conditions are not met, the iteration continues, and the wave field propagation is carried out by calculating the accompanying earthquake source and loading it into the wave equation, and then the low wave number gradient component is calculated by cross-correlation to update the background velocity field, thereby realizing the gradual optimization of waveform inversion. Preferably, in this embodiment, the preset conditions are: the normalized deformation field Φ - nor(t, h) < 1%; when the observed data and the simulated data are consistent, the normalized deformation field is 0, and when the observed data and the simulated data are significantly different, the normalized deformation field exceeds 1. Here we set Φ _ nor(t, h)<1% indicates that the deformation range is less than 1% and the iteration is terminated. Otherwise, the process goes to step S108 to update the speed.

[0091] The entire process fully utilizes the advantages of deep learning in image registration and combines it with expertise in seismic data processing to form an efficient and accurate seismic waveform inversion method.

[0092] The present invention significantly improves the accuracy and efficiency of seismic data processing: through deep learning image registration technology, the deformation field of seismic records can be accurately obtained, providing a more accurate basis for waveform inversion; by iteratively updating the background velocity field, the waveform inversion results are continuously optimized, so that the final seismic image is closer to the actual situation; the present invention has good flexibility and adaptability, and can adjust parameters according to specific data characteristics and requirements to meet the seismic data processing needs in different scenarios.

[0093] According to one embodiment of the present invention, step S104 includes the following steps:

[0094] Step S1042, inputting the image to be registered, and outputting the deformation vector field after learning by the multi-scale parallel full convolutional neural network; wherein the image to be registered includes a floating image and a fixed image;

[0095] The images to be registered include floating images (simulated seismic recording images) and fixed images (seismic observation images). Fully convolutional neural network (FCN) is a deep learning model for image-to-image conversion tasks, suitable for tasks such as image segmentation and image registration. The multi-scale parallel structure means that the network processes images in parallel at different scales, which helps to capture deformation information at different scales. The deformation vector field output by the network after learning describes the displacement vector of each pixel from the floating image to the fixed image, that is, the deformation information between the images.

[0096] Step S1044, the floating image generates a deformable image using a spatial transformation network based on the deformation vector field;

[0097] The spatial transformer network is a learnable module that can spatially transform an image based on an input deformation vector field. The deformed image obtained by transforming the floating image using the deformation vector field is morphologically closer to the fixed image.

[0098] Step S1046, calculating the difference between the fixed image and the deformed image according to the loss function and updating the parameters of the earthquake record deformation field, performing multiple iterations to obtain the final earthquake record deformation field.

[0099] Through optimization methods such as back propagation algorithm and gradient descent, the parameters of the multi-scale parallel full convolutional neural network are updated according to the gradient of the loss function, so that the deformation vector field output by the network is more accurate. The difference between the fixed image and the deformed image is gradually reduced through iterative optimization, and finally an accurate earthquake record deformation field is obtained. Preferably, in this embodiment, the loss function adopts the following formula,

[0100] L(I f , I m ,Φ)=L(I f , I m ,Φ)+λL smoth (I f , I m , Φ) (5)

[0101] Among them, I f Represents a fixed image; I m represents floating image; Φ represents earthquake record deformation field Φ(t,h); L(I f , I m , Φ) represents I m The grayscale similarity of the transformed image; λ represents the spatial smoothness coefficient; L smoth (I f , I m , Φ) indicates smoothing of image grayscale similarity.

[0102] The present invention realizes the precise registration of simulated seismic recording images and seismic acquisition observation images by introducing a multi-scale parallel fully convolutional neural network and a spatial transformation network. It can capture deformation information of different scales, gradually reduce the differences between images through iterative optimization, and finally obtain an accurate seismic recording deformation field, which provides strong technical support for seismic waveform inversion and data analysis, and improves the accuracy and efficiency of seismic data analysis.

[0103] According to one embodiment of the present invention, in step S1042:

[0104] A three-layer multi-scale parallel downsampling module is used to reduce the image size and obtain a multi-scale downsampling feature map for the image to be registered;

[0105] The multi-scale downsampled feature maps are sequentially passed through the pyramid dilated convolution module and the adaptive channel attention module for feature extraction and feature processing;

[0106] The feature map processed by the adaptive channel attention module is input into the direction perception layer and the upsampling layer to obtain the deformation vector field of the input image pair.

[0107] In this implementation, based on the structure of a fully convolutional neural network (FCN), a three-layer multi-scale parallel downsampling module is used to reduce the image size and obtain a multi-scale downsampling feature map. Multi-scale feature extraction is performed through a pyramid dilated convolution module, and then three layers of adaptive channel attention are used for feature processing. Finally, a dense vector deformation field of the input image pair is obtained through a direction perception layer and an upsampling layer.

[0108] Multi-scale parallel downsampling model: In order to consider the need for extracting larger deformation features in the seismic image registration task, a pyramid pooling module (PPM) is used to form a multi-scale parallel downsampling module to extract multi-scale features. This module first performs average pooling of the input image with a step size of 1, 2, and 4 to obtain feature maps of three different scales. Then, through the BN layer and ReLu activation function, the three feature maps and the input image are simultaneously bilinearly interpolated to obtain feature maps with a resolution of 1 / 2 of the original image at different scales. Finally, the three feature maps of different scales and the original feature map are feature-connected to output a cascade feature map with a channel number of 4xc. This technology mainly adopts a multi-scale strategy and can sample large-scale and small-scale deformation features.

[0109] Pyramid convolution module: After multi-scale parallel downsampling, ordinary convolution is difficult to process feature maps containing many different scales of information at the same time, and the downsampling uses the pooling layer to concentrate the features, which easily causes information loss and affects the registration effect. Therefore, the pyramid hole convolution module is used for feature extraction to solve the above problems. It needs to set a dilation factor to control the holes in the convolution kernel. When the dilation factor is 1, the hole convolution is the same as the ordinary convolution operation; when the dilation factor is greater than 1, each operator in the convolution kernel will skip the corresponding number of points, forming a larger and sparse convolution kernel. This technology uses a multi-layer convolution structure in conjunction with a multi-scale strategy to capture image deformation features of different scales.

[0110] Adaptive attention module: The attention mechanism can guide the model to focus on feature maps with large amounts of information, improving model performance while avoiding increasing network depth. Since the model extracts features through a multi-scale parallel downsampling module and a pyramid hole convolution layer, the number of feature maps expressing large deformations (such as strong reflection axes, etc.) is significantly higher than that of small features (such as high-frequency earthquake axes, etc.). In order to express the features of small deformations to improve the registration accuracy, an adaptive channel attention module is used to pass the adaptive channel attention to the feature images after three downsamplings, appropriately increase the weight of small deformations, and improve the registration accuracy.

[0111] According to one embodiment of the present invention, step S108 includes the following steps:

[0112] The accompanying earthquake source is calculated according to the following formula:

[0113]

[0114] Where, (t) represents the time direction; (h) represents the spatial direction; Δτ represents the arrival time difference between the observed data and the simulated data as the waveform time of the image sequence; Δs represents the arrival time difference between the observed data and the simulated data as a waveform space of the image sequence; g t (t+Δτ, h) represents the gradient used for velocity update; represents the first-order time derivative of the observed data; represents the second-order time derivative of the observed data; d cal represents simulation data; w is the proportional coefficient, which is used to balance the proportional relationship between time and space;

[0115] represents the time-delayed accompanying source; It indicates that the accompanying source of the earthquake record causes the vertices of the diffraction waves to move laterally due to the lateral change of velocity;

[0116] The companion source is loaded into the wave equation for wave field propagation, and the low wave number gradient component is calculated by cross-correlation to update the background velocity field. Specifically, the following formula is used:

[0117]

[0118] Where, f(t, h) represents the source wave field; g t (t+Δτ, h) represents the accompanying source wave field; represents the gradient, represents the cross-correlation of the two at each time in the underground; J represents the objective function, and m represents the particle velocity.

[0119] The present invention utilizes a deep learning network to learn the deformation vector field between seismic images, and realizes high-precision inversion of underground velocity structure through accompanying earthquake source calculation and wave equation wave field propagation. It can significantly improve the accuracy and efficiency of seismic waveform inversion, and provide strong technical support for seismic exploration.

[0120] Embodiment 2

[0121] In the process of seismic full waveform inversion, it is necessary to minimize the difference between observed records and simulated records. The key is to solve the geometric transformation from simulated data to observed data. Figure 2 The process of the seismic waveform inversion method based on deep learning image registration according to one embodiment of the present invention is as follows Figure 2 ; Figure 3 This is a comparison diagram of a non-rigid image registration shot collection test according to an embodiment of the present invention, wherein the left side of the image is the original seismic shot collection, the middle of the image is the distorted deformed shot collection, and the right side of the image is the distorted restored shot collection; Figure 4 A network structure diagram of a multi-scale downsampling module according to an embodiment of the present invention; Figure 5 Schematic diagram of a multi-scale parallel fully convolutional neural network according to an embodiment of the present invention. Figure 2-5 As shown, according to one embodiment of the present invention, the present invention provides a seismic waveform inversion method based on deep learning image registration, in which the observed seismic record is set as a reference image, and the simulated seismic record is set as a floating image, wherein the objective function is defined as the matching residual of the floating and reference images, and a deep learning fully convolutional neural network is used to calculate the geometric transformation field for full waveform inversion.

[0122] The program includes four key steps:

[0123] Step 1: Use simulated data to calculate velocity updates, use small offset observation data and background velocity field to perform least square migration to establish a reflection coefficient model, and use reflection coefficient and background velocity field to perform Born forward modeling to obtain simulated seismic record images;

[0124] Step 2: Registration, using the simulated seismic record image as a floating image and the observation image collected by seismic acquisition as a fixed image for deep learning image registration to obtain the seismic record deformation field Φ(t, h), which is the key point of the present invention.

[0125] Step 3: Calculate the normalized deformation field Φ_nor(t,h) = Φ(t,h) / wavelet_length(t,h), where wavelet_length(t,h) represents the wavelength of the main frequency wave of the current data. When the observed data and the simulated data are consistent, the normalized deformation field is 0. When the observed data and the simulated data are significantly different, the normalized deformation field exceeds 1. Here, Φ_nor(t,h) is set to <1%, which means that the iteration is terminated if the deformation range is less than 1%. Otherwise, go to step 4 to update the velocity.

[0126] Step 4: Velocity update: calculate the conjugate source of the deformation field of the seismic record, load the accompanying source into the wave equation for wave field propagation, use cross-correlation to calculate the low-wavenumber gradient component to update the background velocity field, and separate and calculate the high-frequency seismic image to update the high-frequency reflection coefficient.

[0127] According to one embodiment of the present invention, step 1 comprises the following steps:

[0128] Step 1.1: Use least squares migration to improve the consistency of amplitude and waveform between observed seismic images and simulated seismic images; the least squares migration is implemented by using migration deconvolution method.

[0129] Step 1.2: Use the Born forward modeling technique to input the reflection coefficient image and background velocity model to simulate the seismic record image.

[0130] According to one embodiment of the present invention, step 2 comprises the following steps:

[0131] Step 2.1: First, input the images to be registered (floating image and fixed image), and after learning with a multi-scale parallel fully convolutional neural network, output a deformation vector field, which represents the displacement of the corresponding positions in the floating image and the fixed image;

[0132] In step 2.1, based on the structure of the Fully Convolutional Networks (FCN), a three-layer multi-scale parallel downsampling module is used to reduce the image size and obtain a multi-scale down-sampled feature map. Multi-scale feature extraction is performed through a pyramid hole convolution module, and then three layers of adaptive channel attention are used for feature processing. Finally, a dense vector deformation field of the input image pair is obtained through a direction perception layer and an upsampling layer. Specifically,

[0133] Multi-scale parallel downsampling model: In order to consider the need for extracting larger deformation features in the seismic image registration task, a pyramid pooling module (PPM) is used to form a multi-scale parallel downsampling module to extract multi-scale features. This module first performs average pooling of the input image with a step size of 1, 2, and 4 to obtain feature maps of three different scales. Then, through the BN layer and ReLu activation function, the three feature maps and the input image are simultaneously bilinearly interpolated to obtain feature maps with a resolution of 1 / 2 of the original image at different scales. Finally, the three feature maps of different scales and the original feature map are feature-connected to output a cascade feature map with a channel number of 4xc. This technology mainly adopts a multi-scale strategy and can sample large-scale and small-scale deformation features.

[0134] Pyramid convolution module: After multi-scale parallel downsampling, ordinary convolution is difficult to process feature maps containing many different scales of information at the same time, and the downsampling uses the pooling layer to concentrate the features, which easily causes information loss and affects the registration effect. Therefore, the pyramid hole convolution module is used for feature extraction to solve the above problems. It needs to set a dilation factor to control the holes in the convolution kernel. When the dilation factor is 1, the hole convolution is the same as the ordinary convolution operation; when the dilation factor is greater than 1, each operator in the convolution kernel will skip the corresponding number of points to form a larger and sparse convolution kernel. This technology uses a multi-layer convolution structure in conjunction with a multi-scale strategy to capture image deformation characteristics of different scales.

[0135] Adaptive attention module: The attention mechanism can guide the model to focus on feature maps with large amounts of information, improving model performance while avoiding increasing network depth. Since the model extracts features through a multi-scale parallel downsampling module and a pyramid hole convolution layer, the number of feature maps expressing large deformations (such as strong reflection axes, etc.) is significantly higher than that of small features (such as high-frequency earthquake axes, etc.). In order to express the features of small deformations to improve the registration accuracy, an adaptive channel attention module is used to pass the adaptive channel attention to the feature images after three downsamplings, appropriately increase the weight of small deformations, and improve the registration accuracy.

[0136] Step 2.2: Then the floating image generates a new deformed image using the spatial transformation network according to the deformation vector field;

[0137] Step 2.3: Finally, the difference between the fixed image and the deformed image is calculated based on the loss function (including similarity measurement and smoothness constraint) and the parameters of the model are updated. After multiple iterations, the final deformable model is obtained.

[0138] The main goal of image registration is to obtain the deformation field Φ(t, h) and find the floating image I m and a fixed image I f, the pixel-by-pixel correspondence between them. The present invention designs the grayscale value of the image to maximize an objective function about the matching degree between images, L(I f , I m , Φ) represents I m The grayscale similarity of the transformed image, L smoth (I f , I m , Φ) represents the smoothing of image grayscale similarity, and λ is a coefficient representing the degree of spatial smoothness.

[0139] L(I f , I m ,Φ)=L(I f , I m ,Φ)+λL smoth (I f , I m , Φ)

[0140] According to an embodiment of the present invention, step 4 includes: calculating the adjoint source according to the deformation field Φ(t, h) obtained in step 2; specifically,

[0141]

[0142] Where, (t) represents the time direction; (h) represents the spatial direction; Δτ represents the arrival time difference between the observed data and the simulated data as the waveform time of the image sequence; Δs represents the arrival time difference between the observed data and the simulated data as a waveform space of the image sequence; gt(t+Δτ,h) represents the gradient used for velocity update; represents the first-order time derivative of the observed data; represents the second-order time derivative of the observed data; d cal Represents simulation data; W is the proportionality coefficient, which is used to balance the proportional relationship between time and space;

[0143] After derivation, the earthquake record with the source can be obtained, where the first term in the above formula is represents the time-delayed companion source, the second term It indicates that the companion source of the seismic record causes the vertices of the diffraction waves to move laterally due to the lateral change of velocity;

[0144] The gradient is:

[0145]

[0146] Where, f(t, h) represents the source wave field; g t (t+Δτ, h) represents the companion source wave field; represents the gradient, and represents the cross-correlation of the two at each time in the underground.

[0147] After determining the objective function in step 4, the nonlinear optimization problem is linearized using Taylor expansion, and the result is obtained step by step through iteration. This step uses the conjugate gradient method to complete the speed update, and selects the reverse direction of the gradient and the direction of the last update to construct the conjugate gradient direction. The update step size is obtained by parabolic interpolation.

[0148] The beneficial effects of the present invention are:

[0149] (1) For the local extrema problem caused by frequency jumps in classical full waveform inversion, the present invention performs deep learning global matching to improve the dependence on the initial velocity model and adopts a deformation field to avoid the uncertainty problem caused by matching amplitude.

[0150] (2) For the calculation of deformation vector fields, the present invention adopts a multi-scale parallel downsampling model and an adaptive channel attention module, which is not only suitable for the calculation of large deformation fields but also has high registration accuracy for small deformations. After training, the efficiency of this method is significantly improved compared with traditional image registration methods.

[0151] (3) The deformation vector field of the present invention is a two-dimensional matrix, which is more reasonable than the one-dimensional time delay obtained by the current methods such as cross-correlation, dynamic time warping, and transport, and is easy to obtain a more accurate velocity model.

[0152] (4) The reflected wave data is obtained by Born forward modeling for matching. The present invention is not only applicable to transmitted waves but also to reflected waves, thus avoiding the problem of insufficient speed update depth of the classical full waveform inversion.

[0153] Taking all factors into consideration, the present invention utilizes the high efficiency and high precision of deep learning to avoid the frequency skipping problem that is prone to occur in waveform matching and falls into the local minimum; the present invention can construct a more stable velocity model than conventional methods, thereby significantly improving the quality of seismic migration imaging.

[0154] Embodiment 3

[0155] like Figure 6-Figure 12 As shown, this embodiment is based on the practical application of embodiment 2.

[0156] Fig.11 (a) is a true velocity profile of a 2D depression model, 10000 m long and 8000 m deep, with a 2D seismic observation system, 200 m shot spacing, and 10000 m offset length; Fig. 9 An initial model data shot gather associated with a seismic source in one embodiment of the present invention, Fig. 9 right Fig.11(a) The model simulates the point source acoustic wave equation. The source function is a Ricker wavelet with a main frequency of 8 Hz. The observation data is 40 shot points evenly distributed on the surface. The shot point spacing is 250 m. Each shot point contains 801 detectors with a spacing of 12.5 m. The earthquake record formed by surface observation.

[0157] Test use Figure 2 The technical flow chart shown in the figure performs velocity inversion test on the data. Fig.11 (b) is the initial velocity model, which is the initial velocity model used for testing. The fully convolutional neural network structure used in this experiment is shown in the figure. Figure 5 shown.

[0158] Figure 6 The initial model observation data shot gather of one embodiment of the present invention is an original shot gather from which the transmitted wave has been cut off, wherein the curvature of the reflection event axis is complex and there are a large number of intersections, and the complexity is relatively high.

[0159] Figure 7 The de-migration simulated shot gather in one embodiment of the present invention is a de-migration simulated shot gather obtained by performing least square migration in step S202 and forward modeling using the initial velocity Born, and has a shape similar to the original shot gather;

[0160] Figure 8 The initial model data shot gather deformation field of one embodiment of the present invention is a deformation field calculated using a fully convolutional neural network, which describes the transformation relationship between the reverse migration simulation data and the original shot gather fixed data. It can be seen from the figure that the deformation field changes very little at near offset distances, while the deformation field values ​​at bilateral far offset distances are relatively large.

[0161] Fig. 9 The deformation field is used to calculate the accompanying source. In this experiment, the accompanying source will be propagated in the detector wave field in the reverse direction and the velocity gradient field will be calculated.

[0162] Fig.10 An initial model velocity update gradient according to an embodiment of the present invention; Fig.11 (c) is the velocity model after model inversion update. After completing 5 velocity iterations using this method, the reflection wave waveform inversion results are as follows: Fig.11 As shown in (c), it can be seen that this method can better restore the low-frequency structural information of the velocity model. In this experiment, the three velocity models were migrated and imaged. Fig.12 (a) is the migration result of the true velocity model. Fig.12 (b) is the migration result of the initial velocity model, and the initial velocity model has poor imaging; Fig.12 (c) is the migration result of the inverted velocity model. It is the migration result of the inverted velocity of this experiment. It can be seen that there is a significant improvement compared with the initial model imaging 11(b), and most of the structures are well imaged.

[0163] Embodiment 4

[0164] Fig.13 FIG. 1 is a structural diagram of a seismic waveform inversion system based on deep learning image registration according to an embodiment of the present invention. Fig.13 As shown, a seismic waveform inversion system based on deep learning image registration includes:

[0165] The simulated seismic record image acquisition module 10 is used to perform least square migration based on small offset observation data and background velocity field to construct a reflection coefficient model, and to perform Born forward modeling based on the reflection coefficient and background velocity field to obtain simulated seismic record images;

[0166] A seismic record deformation field acquisition module 20 is used to acquire the seismic record deformation field by using deep learning image registration based on simulated seismic record images and seismic observation images;

[0167] The seismic waveform inversion module 30 is used to calculate the normalized deformation field using the following formula:

[0168] Φ_nor(t,h)=Φ(t,h) / wavelet_length(t,h)

[0169] Among them, Φ_nor(t, h) represents the normalized deformation field; Φ(t, h) represents the earthquake record deformation field; wavelet_length(t, h) represents the wavelength of the main frequency wave of the current data, (t) represents the time direction, and (h) represents the space direction;

[0170] The seismic waveform inversion module 30 is also used to determine whether the normalized deformation field meets the preset conditions; if so, the iteration process is terminated; otherwise, the background velocity field is updated through the velocity update module;

[0171] The velocity updating module 40 is used to calculate the accompanying earthquake source according to the deformation field of the earthquake record, load the accompanying earthquake source into the wave equation for wave field propagation, and calculate the low wave number gradient component by cross-correlation to update the background velocity field.

[0172] An electronic device comprises: a processor, a memory, and a computer program stored in the memory and executable on the processor, wherein when the computer program is executed by the processor, any seismic waveform inversion method based on deep learning image registration is implemented.

[0173] Those skilled in the art can clearly understand that, for the convenience and brevity of description, the specific working process of the system and device described above can refer to the corresponding process in the aforementioned method embodiment, and will not be repeated here.

[0174] The above description is only a preferred embodiment of the present application and an explanation of the technical principles used. Those skilled in the art should understand that the scope of the invention involved in the present application is not limited to the technical solution formed by a specific combination of the above technical features, but should also cover other technical solutions formed by any combination of the above technical features or their equivalent features without departing from the inventive concept. For example, the above features are replaced with (but not limited to) technical features with similar functions disclosed in the present application.

[0175] It should be understood that the size of the serial numbers of each step in the content of the invention and the embodiments of the present invention does not absolutely mean the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of the present invention.

Claims

1. A seismic waveform inversion method based on deep learning image registration, characterized in that: The following steps are involved: Step S102, performing least square migration based on small offset observation data and background velocity field to construct a reflection coefficient model, and performing Born forward modeling based on the reflection coefficient and the background velocity field to obtain a simulated seismic recording image; Step S104, based on the simulated earthquake record image and the earthquake observation image, using deep learning image registration to obtain the earthquake record deformation field; Step S106, calculate the normalized deformation field using the following formula: Φ_nor(t,h)=Φ(t,h) / wavelet_length(t,h) Among them, Φ_nor(t, h) represents the normalized deformation field; Φ(t, h) represents the earthquake record deformation field; wave1et_1ength(t, h) represents the wavelength of the main frequency wave of the current data, (t) represents the time direction, and (h) represents the space direction; Determine whether the normalized deformation field meets the preset condition; if so, the iteration process is terminated; otherwise, execute step S108; Step S108, calculating the accompanying seismic source according to the earthquake record deformation field, loading the accompanying seismic source into the wave equation for wave field propagation, using cross-correlation to calculate the low wave number gradient component to update the background velocity field, and repeating steps S102 to S106.

2. The seismic waveform inversion method based on deep learning image registration according to claim 1, characterized in that: In step S102, the least squares shift is implemented by using a shift deconvolution method.

3. The seismic waveform inversion method based on deep learning image registration according to claim 1, characterized in that: In step S104, the simulated earthquake record image is used as a floating image, and the earthquake observation image is used as a fixed image, and deep learning image registration is performed to obtain the earthquake record deformation field.

4. The seismic waveform inversion method based on deep learning image registration according to claim 3, characterized in that: Step S104 includes the following steps: Input an image to be registered, and output a deformation vector field after learning with a multi-scale parallel fully convolutional neural network; wherein the image to be registered includes the floating image and the fixed image; The floating image generates a deformable image based on the deformation vector field using a spatial transformation network; The difference between the fixed image and the deformed image is calculated according to the loss function, and the parameters of the seismic record deformation field are updated, and multiple iterations are performed to obtain the final seismic record deformation field.

5. The seismic waveform inversion method based on deep learning image registration according to claim 4, characterized in that: Using a three-layer multi-scale parallel downsampling module to reduce the image size and obtain a multi-scale downsampling feature map for the image to be registered; The multi-scale down-sampled feature map is sequentially passed through a pyramid dilated convolution module and an adaptive channel attention module for feature extraction and feature processing; The feature map processed by the adaptive channel attention module is input into the direction perception layer and the upsampling layer to obtain the deformation vector field of the input image pair.

6. The seismic waveform inversion method based on deep learning image registration according to claim 4, characterized in that: The loss function adopts the following formula: THERE f ,THE m ,Φ)=I(I f ,THE m ,Φ)+λL smotμ (THE f ,THE m (Φ) Among them, I f Represents a fixed image; I m represents floating image; Φ represents earthquake record deformation field Φ(t,h); L(I f , I m , Φ) represents I m The grayscale similarity of the transformed image; λ represents the spatial smoothness coefficient; L smoth (I f , I m , Φ) indicates smoothing of image grayscale similarity.

7. The seismic waveform inversion method based on deep learning image registration according to claim 1, characterized in that: In step S106: The preset condition is: the normalized deformation field Φ_nor(t, h)<1%.

8. The seismic waveform inversion method based on deep learning image registration according to claim 1, characterized in that: Step S108 includes the following steps: The accompanying earthquake source is calculated according to the following formula: Where, (t) represents the time direction; (h) represents the spatial direction; Δτ represents the arrival time difference between the observed data and the simulated data as the waveform time of the image sequence; Δs represents the arrival time difference between the observed data and the simulated data as the waveform space of the image sequence; gt(t+Δτ, h) represents the gradient used for velocity update; represents the first-order time derivative of the observed data; represents the second-order time derivative of the observed data; d cal Represents simulation data; W is the proportionality coefficient, which is used to balance the proportional relationship between time and space; represents the time-delayed accompanying source; It indicates that the accompanying source of the earthquake record causes the vertices of the diffraction waves to move laterally due to the lateral change of velocity; The companion source is loaded into the wave equation for wave field propagation, and the low wave number gradient component is calculated by cross-correlation to update the background velocity field, specifically through the following formula: The gradient is: Among them, f (t，h) represents the source wave field; g t (t+Δτ, h) represents the accompanying source wave field; represents the gradient, represents the cross-correlation of the two at each time in the underground; J represents the objective function, and m represents the particle velocity.

9. A seismic waveform inversion system based on deep learning image registration, characterized in that: include: A simulated seismic record image acquisition module is used to perform least square migration based on small offset observation data and background velocity field to construct a reflection coefficient model, and to perform Born forward modeling based on the reflection coefficient and the background velocity field to obtain a simulated seismic record image; A seismic record deformation field acquisition module, used to acquire the seismic record deformation field by using deep learning image registration based on the simulated seismic record image and the seismic observation image; The seismic waveform inversion module is used to calculate the normalized deformation field using the following formula, Φ_nor(t,h)=Φ(t,h) / wavelet_length(t,h) Among them, Φ_nor(t, h) represents the normalized deformation field; Φ(t, h) represents the earthquake record deformation field; wavelet_length(t, h) represents the wavelength of the main frequency wave of the current data, (t) represents the time direction, and (h) represents the space direction; The seismic waveform inversion module is also used to determine whether the normalized deformation field meets the preset conditions; if so, the iteration process is terminated; otherwise, the background velocity field is updated through the velocity update module; The velocity update module is used to calculate the accompanying earthquake source according to the earthquake record deformation field, load the accompanying earthquake source into the wave equation for wave field propagation, and calculate the low wave number gradient component by cross-correlation to update the background velocity field.

10. An electronic device, characterized in that: include: A processor, a memory, and a computer program stored in the memory and executable on the processor, wherein the computer program, when executed by the processor, implements the seismic waveform inversion method based on deep learning image registration as described in any one of claims 1 to 8.

Images