Elastic wave full waveform inversion method based on vertical seismic profile

Abstract

The application provides an elastic wave full waveform inversion method based on a vertical seismic profile, and the method comprises the following steps: obtaining original shot gather data; performing horizontal component rotation on a first original component and a second original component to obtain a first radial component and a tangential component; performing wave field separation on a third original component and the first radial component based on multi-channel vector median filtering to obtain a first downgoing wave, a second downgoing wave, a first upgoing wave and a second upgoing wave; performing amplitude compensation on the first downgoing wave, the second downgoing wave, the first upgoing wave and the second upgoing wave, and then stacking to obtain a fourth original component and a second radial component; and performing consistency correction on the fourth original component and the second radial component based on a root mean square value to obtain a target original component and a target radial component. The application improves the accuracy and reliability of subsequent elastic wave full waveform inversion through a fine pretreatment step.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of seismic exploration, and particularly relates to an elastic wave full waveform inversion method based on a vertical seismic profile. BACKGROUND

[0002] The elastic wave full waveform inversion (EFWI) technology can provide high-fidelity full wave field information, including uplink and downlink transverse waves and longitudinal waves, and this characteristic enables the technology to realize accurate velocity modeling, thereby improving the accuracy of depth migration.

[0003] However, the EFWI has very strict requirements for data preprocessing, and the preprocessing technology must be able to maintain the vector characteristics of the seismic wave amplitude to ensure the accuracy of the inversion, in addition, the EFWI is also very sensitive to the initial velocity model, and an accurate initial velocity model is needed to avoid deviation of the inversion result. SUMMARY

[0004] Therefore, the present application aims to provide an elastic wave full waveform inversion method based on a vertical seismic profile.

[0005] As an aspect of the present application, an elastic wave full waveform inversion method based on a vertical seismic profile is provided, comprising:

[0006] obtaining original shot gather data; the original shot gather data includes a first original component, a second original component and a third original component obtained based on three-component orientation;

[0007] rotating the first original component and the second original component to obtain a first radial component and a tangential component;

[0008] based on multi-channel vector median filtering, performing wave field separation on the third original component and the first radial component to obtain a first downlink wave, a second downlink wave, a first uplink wave and a second uplink wave;

[0009] amplitude compensation is performed on the first downlink wave, the second downlink wave, the first uplink wave and the second uplink wave to obtain a fourth original component and a second radial component;

[0010] based on a root mean square value, the fourth original component and the second radial component are corrected for consistency to obtain a target original component and a target radial component.

[0011] Optionally, an initial velocity model is determined according to the target original component and the target radial component; the initial velocity includes a longitudinal wave velocity model and a transverse wave velocity model.

[0012] determine an initial density model according to the P-wave velocity model;

[0013] determine a target initial model according to the initial velocity model and the initial density model;

[0014] determine first wave field data according to the target initial model;

[0015] obtain second wave field data; the second wave field data is actual wave field data;

[0016] determine a residual value according to the first wave field data and the second wave field data;

[0017] the residual value satisfies the following expression:

[0018]

[0019] wherein E(c) is the residual value, x s is a source position, x r is a receiver position, t is time, u cal (x s , x r , t) is the first wave field data, u obs (x s , x r , t) is the second wave field data;

[0020] iteratively update model parameters of the target initial model based on the residual value;

[0021] stop iteration in response to the target initial model reaching a preset condition; the preset condition includes at least one of the following: the number of iterations reaches a first threshold value and the residual value is less than or equal to a second threshold value.

[0022] Optionally, the determining an initial velocity model according to the target original component and the target radial component comprises:

[0023] obtain a first arrival time of the target original component and the target radial component;

[0024] determine the initial velocity model based on the first arrival time.

[0025] Optionally, the iteratively updating model parameters of the target initial model based on the residual value comprises:

[0026] obtain model vectors of the nth iteration and the n+1th iteration;

[0027] determine a step length of the nth iteration based on the model vectors of the nth iteration and the n+1th iteration and the residual value;

[0028] The step length of the nth iteration satisfies the following expression:

[0029] p n = m n+1 - m n ;

[0030]

[0031]

[0032] wherein m n and m n+1 are model vectors of the nth and (n+1)th iterations, p n is a model update vector of the (n+1)th iteration, q n is a gradient difference vector of the (n+1)th iteration, a n is a step length of the nth iteration, is a transpose of p n-1 , ‖q n-1 ‖2 is a two-norm of q n-1 ;

[0033] Based on the step length of the nth iteration and the residual value, the model parameters of the target initial model are iteratively updated based on a steepest gradient descent method.

[0034] The model vector of the (n+1)th iteration satisfies the following expression:

[0035]

[0036] Optionally, the initial velocity model is verified based on a corridor stacking method and a pre-stack depth migration method.

[0037] Optionally, the third original component and the first radial component are denoised based on a same noise attenuation parameter.

[0038] Optionally, the amplitude compensation includes spherical divergence compensation and absorption attenuation compensation.

[0039] As a second aspect of the present application, an elastic wave full waveform inversion device based on a vertical seismic profile is provided, comprising an acquisition module, a processing module and a correction module.

[0040] The acquisition module is configured to acquire original shot gather data; the original shot gather data includes a first original component, a second original component and a third original component obtained based on three-component orientation;

[0041] The processing module is configured to perform horizontal component rotation on the first original component and the second original component to obtain a first radial component and a tangential component;

[0042] The processing module is further configured to perform wave field separation on the third original component and the first radial component based on a multi-channel vector median filter to obtain a first downgoing wave, a second downgoing wave, a first upgoing wave and a second upgoing wave;

[0043] The processing module is further configured to perform amplitude compensation on the first downgoing wave, the second downgoing wave, the first upgoing wave and the second upgoing wave and superimpose to obtain a fourth original component and a second radial component;

[0044] The correction module is configured to perform consistency correction on the fourth original component and the second radial component based on a root mean square value to obtain a target original component and a target radial component.

[0045] Optionally, the elastic wave full waveform inversion device based on a vertical seismic profile further comprises a determination module and an updating module.

[0046] The determination module is configured to determine an initial velocity model according to the target original component and the target radial component; the initial velocity comprises a P-wave velocity model and a S-wave velocity model.

[0047] The determination module is further configured to determine an initial density model according to the P-wave velocity model.

[0048] The determination module is further configured to determine a target initial model according to the initial velocity model and the initial density model.

[0049] The determination module is further configured to determine first wave field data according to the target initial model.

[0050] The acquisition module is configured to acquire second wave field data; the second wave field data is actual wave field data.

[0051] The determination module is further configured to determine a residual value according to the first wave field data and the second wave field data.

[0052] The residual value satisfies the following expression:

[0053]

[0054] wherein E(c) is the residual value, x s is a source position, x r is a receiver position, t is time, u cal (x s ,x r ,t) is the first wave field data, u obs (x s ,x r ,t) is the second wave field data.

[0055] The updating module is configured to iteratively update model parameters of the target initial model based on the residual value.

[0056] The processing module is further configured to stop iteration in response to the target initial model reaching a preset condition, the preset condition including at least one of a first threshold of iteration times and a second threshold of the residual value being less than or equal to.

[0057] Optionally, the determining module is specifically configured to acquire a first arrival time of the target original component and the target radial component, and determine the initial velocity model based on the first arrival time.

[0058] Optionally, the updating module is specifically configured to acquire model vectors of the n th iteration and the n+1 th iteration.

[0059] The updating module is configured to iteratively update model parameters of the target initial model based on the residual value.

[0060] The n th iteration step length satisfies the following expression:

[0061] p n =m n+1 -m n ;

[0062]

[0063]

[0064] wherein, m n and m n+1 are the model vectors of the n th iteration and the n+1 th iteration, p n is a model update vector of the n+1 th iteration, q n is a gradient difference vector of the n+1 th iteration, a n is the n th iteration step length, is a transpose of p n-1 , ||q n-1 ||2 is a two norm of q n-1 ;

[0065] The updating module is configured to iteratively update model parameters of the target initial model based on the residual value.

[0066] The n+1 th iteration model vector satisfies the following expression:

[0067]

[0068] Optionally, the vertical seismic profile based elastic wave full waveform inversion device further comprises a verification module.

[0069] The verification module is configured to verify the initial velocity model based on a corridor stack method and a pre-stack depth migration method.

[0070] Optionally, the third original component and the first radial component are denoised based on a same noise attenuation parameter.

[0071] Optionally, the amplitude compensation comprises spherical divergence compensation and absorption attenuation compensation.

[0072] As a third aspect of the present application, an electronic device comprises a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor implements the above-mentioned vertical seismic profile based elastic wave full waveform inversion method when executing the program.

[0073] As a fourth aspect of the present application, a non-transitory computer readable storage medium is provided, which stores computer instructions for causing the computer to execute the above-mentioned vertical seismic profile based elastic wave full waveform inversion method provided by the present application.

[0074] As can be seen from the above, the vertical seismic profile based elastic wave full waveform inversion method provided by the present application replaces three-component directional processing of original shot gather data with horizontal component rotation, converts seismic wave components from an original coordinate system to a radial and tangential coordinate system which is more convenient for analysis, provides convenience for subsequent wave field separation, and reduces errors that may be caused by mismatch of coordinate systems. Furthermore, the multi-channel vector median filtering technology is used to perform wave field separation on the third original component and the first radial component, which can effectively distinguish different types of wave vectors, including downgoing waves and upgoing waves, and helps to more clearly identify and analyze the propagation path and characteristics of seismic waves. Further, the downgoing waves and the upgoing waves obtained by separation are subjected to amplitude compensation and stacking to obtain a fourth original component and a second radial component. The amplitude compensation can correct amplitude attenuation caused by factors such as propagation distance and medium attenuation, and the stacking helps to retain the original vector characteristics and improve the signal-to-noise ratio of the data. Finally, the fourth original component and the second radial component are subjected to consistency correction based on a root mean square value to obtain a target original component and a target radial component. This step helps to eliminate differences and errors that may be caused in the data acquisition and processing process, and ensures the consistency between different components. Based on the above-mentioned fine preprocessing steps, high-quality input data are provided for subsequent elastic wave full waveform inversion, which helps to reduce the influence of noise, retains the vector characteristics of the target original component and the target radial component, improves the accuracy and reliability of inversion, and thus avoids deviation of the inversion result. BRIEF DESCRIPTION OF DRAWINGS

[0075] In order to more clearly illustrate the technical solutions in the application or the related art, the accompanying drawings needed to be used in the embodiments or the related art description will be briefly introduced. Obviously, the accompanying drawings in the following description only constitute the embodiments of the application, and for those skilled in the art, other drawings can also be obtained without creative labor based on these drawings.

[0076] Figure 1 A flowchart of an elastic wave full waveform inversion method based on a vertical seismic profile is provided for the embodiments of the application.

[0077] Figure 2 A schematic diagram of a raw shot record is provided for the embodiments of the application.

[0078] Figure 3 A schematic diagram of a shot record after horizontal rotation is provided for the embodiments of the application.

[0079] Figure 4 A schematic diagram of a filtered shot record and filtered noise is provided for the embodiments of the application.

[0080] Figure 5 A schematic diagram of up-and-down wave field records of R components is provided for the embodiments of the application.

[0081] Figure 6 A schematic diagram of up-and-down wave field records of Z components is provided for the embodiments of the application.

[0082] Figure 7 A schematic diagram of source energy consistency correction is provided for the embodiments of the application.

[0083] Figure 8 Another flowchart of an elastic wave full waveform inversion method based on a vertical seismic profile is provided for the embodiments of the application.

[0084] Figure 9 Still another flowchart of an elastic wave full waveform inversion method based on a vertical seismic profile is provided for the embodiments of the application.

[0085] Figure 10 A schematic diagram of a target initial model is provided for the embodiments of the application.

[0086] Figure 11 A schematic diagram of an elastic wave full waveform inversion model is provided for the embodiments of the application.

[0087] Figure 12 A schematic diagram of forward recording of an elastic wave full waveform inversion model is provided for the embodiments of the application.

[0088] Figure 13 A flowchart of still another elastic wave full waveform inversion method based on a vertical seismic profile is provided for an embodiment of the present application;

[0089] Figure 14 A pre-stack depth migration profile and corridor stack calibration diagram is provided for an embodiment of the present application;

[0090] Figure 15 A flowchart of still another elastic wave full waveform inversion method based on a vertical seismic profile is provided for an embodiment of the present application;

[0091] Figure 16 A composition diagram of an elastic wave full waveform inversion device based on a vertical seismic profile is provided for an embodiment of the present application;

[0092] Figure 17 A composition diagram of an electronic device is provided for an embodiment of the present application. DETAILED DESCRIPTION

[0093] In order to make the objectives, technical solutions and advantages of the present application clearer, the present application is further described in detail below with reference to the drawings and in conjunction with specific embodiments. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative work fall within the scope of protection of the present application.

[0094] It should be noted that in the embodiments of the present application, the words such as "exemplarily" or "for example" are used to represent as an example, illustration or explanation. Any embodiment or design scheme described as "exemplarily" or "for example" in the embodiments of the present application should not be interpreted as more preferred or more advantageous than other embodiments or design schemes. Rather, the words such as "exemplarily" or "for example" are intended to present the relevant concept in a specific manner. Unless otherwise defined, the technical terms or scientific terms used in the embodiments of the present application should be the usual meanings understood by those of ordinary skill in the art to which the present application belongs. SUMMARY

[0096] In the related art, the EFWI technology can provide more detailed underground velocity and density distribution, which is very valuable for understanding and describing complex underground structures such as faults, fractures and lithology changes. Compared with the traditional seismic exploration method, EFWI can better reveal the geological details of the underground, thereby helping to improve the success rate of exploration and reduce the development risk. However, the EFWI has very strict requirements for data preprocessing, and the preprocessing technology must be able to maintain the vector characteristics of the seismic wave amplitude, which means that not only the amplitude information of the wave but also the polarization direction of the wave should be preserved during the processing. If these vector characteristics are damaged during the preprocessing, it will directly affect the accuracy of the inversion result.

[0097] The inventors of the present application found that although the current first arrival time inversion of P-wave velocity based on vertical seismic profiling (VSP) technology is quite mature, it faces challenges in processing S-wave first arrival. Because the S-wave first arrival is often covered by other wave field signals, it is difficult to extract the S-wave first arrival with high precision. In addition, noise interference is another important factor affecting the stability of inversion. Many existing VSP prestack processing technologies cannot effectively maintain the polarization relationship of elastic wave field in the vertical component and the radial component, which seriously weakens the reliability and accuracy of the inversion result.

[0098] To solve the above problems, the present application provides an elastic wave full waveform inversion method based on vertical seismic profiling. By replacing the three-component directional processing of the original shot gather data with horizontal component rotation, the seismic wave components are converted from the original coordinate system to the radial and tangential coordinate systems which are more convenient for analysis, providing convenience for subsequent wave field separation and reducing the errors that may be caused by coordinate system mismatch. Further, the multi-channel vector median filtering technology is used to separate the wave field of the third original component and the first radial component, which can effectively distinguish different types of wave vectors, including downgoing waves and upgoing waves, helping to more clearly identify and analyze the propagation path and characteristics of seismic waves. Furthermore, the downgoing waves and upgoing waves obtained by separation are amplitude compensated and superimposed to obtain the fourth original component and the second radial component. Amplitude compensation can correct the amplitude attenuation caused by factors such as propagation distance and medium attenuation, while superposition helps to preserve the original vector characteristics and improve the signal-to-noise ratio of the data. Finally, the fourth original component and the second radial component are uniformly corrected based on the root mean square value to obtain the target original component and the target radial component. This step helps to eliminate the differences and errors that may be generated during data acquisition and processing, ensuring the consistency between different components. Based on the above fine preprocessing steps, high-quality input data are provided for subsequent elastic wave full waveform inversion, which helps to reduce the influence of noise while preserving the vector characteristics of the target original component and the target radial component, improving the accuracy and reliability of the inversion, thereby avoiding deviation of the inversion result.

[0099] After introducing the basic principles of the present application, the various non-limiting embodiments of the present application will be specifically introduced.

[0100] Figure 1 A flowchart of an elastic wave full waveform inversion method based on a vertical seismic profile is provided for an embodiment of the present application. As shown in the figure, the elastic wave full waveform inversion method based on the vertical seismic profile provided by the present application specifically includes the following steps: Figure 1

[0101] S101, obtaining original shot gather data.

[0102] The original shot gather data includes a first original component, a second original component and a third original component obtained based on three-component orientation.

[0103] In some embodiments, in seismic exploration, the original shot gather data received by the receiver after being excited by the seismic source and propagating through the underground medium records the amplitude, frequency and phase of the seismic wave at different times and different positions, wherein the three-component geophone can measure the vibration of the seismic wave in three-dimensional space, thereby obtaining three directional components.

[0104] It should be noted that the three directional components refer to the vibration records of the seismic wave in three mutually perpendicular directions, and the three components are X component (first original component), Y component (second original component) and Z component (third original component).

[0105] S102, rotating the first original component and the second original component to obtain a first radial component and a tangential component.

[0106] In some embodiments, as shown in the figure, the X component (a) and the Y component (b) of the original shot gather data are converted by horizontal component, as shown in the figure, to obtain the R component (first radial component) and the T component (tangential component), wherein the longitudinal wave (P wave) and the vertical shear transverse wave (SV wave) are mainly projected on the R component, and the horizontal shear transverse wave (SH wave) is mainly polarized on the T component. Figure 2 Figure 3

[0107] ​​​It should be noted that the horizontal component rotation is a process of rotating two horizontal components (X component and Y component) perpendicular to each other in seismic data from the observation coordinate system to the source-receiver coordinate system. In this process, the data itself does not change, only the observation angle changes. In three-component seismic exploration, due to the change of surface shot point position and the fixed underground receiver, the X component and Y component of the receiver receive SH wave and SV wave. Through horizontal rotation, these wave fields can be separated into R and T two components, so as to more easily identify and analyze the characteristics of seismic waves.

[0108] S103, based on the multi-channel vector median filtering, wave field separation is performed on the third original component and the first radial component to obtain a first downgoing wave, a second downgoing wave, a first upgoing wave and a second upgoing wave.

[0109] In some embodiments, the same noise attenuation parameter is applied to the third original component and the first radial component in the wave field separation process to maintain the polarization relationship of P wave and SV wave in Z and R component amplitude, as shown in Figure 4 As shown in the R component noise (a), Z component noise (b), R component (c) and Z component (d), it can be seen that there is no strong energy of upgoing and downgoing P and S waves on the residual record, which shows that the filtering is amplitude-preserving.

[0110] In some embodiments, the wave field separation decomposes seismic wave data into wave field components in different directions, which is helpful for subsequent data analysis and processing. The multi-channel vector median filtering is used for step-by-step wave field separation. According to the characteristics that the downgoing wave first arrival can be better distinguished, while the upgoing wave phase axis is often difficult to identify, the energy of different types of wave field is sorted from strong to weak for step-by-step separation. The separation order is: downgoing PP wave (first downgoing wave), downgoing PS wave (second downgoing wave), upgoing PP wave (first upgoing wave) and upgoing PS wave (second upgoing wave). The separated R component wave field record: (a) downgoing PP wave, (b) downgoing PS wave, (c) upgoing PP wave and (d) upgoing PS wave as shown in Figure 5 The separated Z component wave field record: (a) downgoing PP wave, (b) downgoing PS wave, (c) upgoing PP wave and (d) upgoing PS wave as shown in Figure 6 .

[0111] It should be noted that the multi-channel vector median filtering is a nonlinear filtering method for processing multi-channel data, which removes noise and interference by sorting and taking median value of the data of each channel.

[0112] It should be understood that the upgoing wave is a wave propagating from the underground medium to the ground surface, and the downgoing wave is a wave propagating from the ground surface to the underground medium.

[0113] S104, amplitude compensation is performed on the first downgoing wave, the second downgoing wave, the first upgoing wave and the second upgoing wave, and superposition is performed to obtain a fourth original component and a second radial component.

[0114] In some embodiments, amplitude compensation is performed on the P wave and the SV wave of the upgoing and downgoing waves separated from the Z and R components, including spherical divergence compensation and absorption attenuation compensation. In addition, for the Z and R components after the wave field separation in the same mode, the same parameters are used for amplitude compensation. After the amplitude compensation, the P wave and the SV wave of the upgoing and downgoing waves have similar waveforms in different time windows, so as to ensure that the wave field with polarization characteristics is obtained.

[0115] It should be noted that the spherical divergence compensation is a compensation method considering the attenuation of the seismic wave due to spherical divergence during propagation, and the absorption attenuation compensation is a compensation method considering the attenuation of the seismic wave due to medium absorption during propagation.

[0116] In some embodiments, according to the propagation distance and medium characteristics of the seismic wave, an amplitude compensation factor is calculated to compensate the amplitudes of the downgoing wave and the upgoing wave, so as to correct the amplitude attenuation caused by the propagation distance and medium attenuation. The wave field components after the compensation are superimposed, that is, the P wave and the SV wave of the upgoing and downgoing waves separated from the Z component are superimposed to obtain a fourth original component with preserved vector characteristics, and the P wave and the SV wave of the upgoing and downgoing waves separated from the R component are superimposed to obtain a second radial component with preserved vector characteristics.

[0117] S105, based on the root mean square value, the fourth original component and the second radial component are uniformly corrected to obtain a target original component and a target radial component.

[0118] In some embodiments, the root mean square (RMS) is a statistical quantity used to measure the fluctuation or dispersion degree of data, which represents the square root of the average value of the square sum of data. By calculating the root mean square value of the fourth original component and the second radial component, the average is calculated, and as shown in Figure 7 the root mean square value is used to uniformly correct the fourth original component and the second radial component, so as to ensure that the energy of the fourth original component and the second radial component is consistent. Further, the corrected data is verified and corrected to ensure the accuracy and reliability of the data.

[0119] The elastic wave full waveform inversion method based on the vertical seismic profile provided in the application is characterized in that: the horizontal component rotation is used to replace the three-component directional processing of the original shot gather data, the seismic wave component is converted from the original coordinate system to the radial and tangential coordinate system which is more convenient for analysis, and the subsequent wave field separation is facilitated, and the error caused by the mismatch of the coordinate system is reduced. Further, the multi-channel vector median filtering technology is used for wave field separation of the third original component and the first radial component, which can effectively distinguish different types of wave vectors, including downgoing waves and upgoing waves, and is helpful for more clearly identifying and analyzing the propagation path and characteristics of the seismic wave. Furthermore, the amplitude compensation is performed on the separated downgoing wave and upgoing wave, and the fourth original component and the second radial component are obtained by superposition. The amplitude compensation can correct the amplitude attenuation caused by the propagation distance, medium attenuation and other factors, and the superposition is helpful for retaining the original vector characteristics and improving the signal-to-noise ratio of the data. Finally, the consistency correction is performed on the fourth original component and the second radial component based on the root mean square value to obtain the target original component and the target radial component. This step is helpful for eliminating the differences and errors that may be generated in the data acquisition and processing process, ensuring the consistency between different components, providing high-quality input data for the subsequent elastic wave full waveform inversion based on the above-mentioned fine preprocessing steps, reducing the influence of noise, retaining the vector characteristics of the target original component and the target radial component, improving the accuracy and reliability of the inversion, and thus avoiding the deviation of the inversion result.

[0120] In some embodiments, as shown in Figure 8 The elastic wave full waveform inversion method based on the vertical seismic profile provided in the embodiments of the application further includes the following S201-S208.

[0121] S201, determining an initial velocity model according to the target original component and the target radial component.

[0122] The initial velocity includes a P-wave velocity model and a S-wave velocity model.

[0123] In some embodiments, the initial velocity model is the basis of seismic wave inversion, which describes the propagation velocity of P-waves and S-waves in the underground medium. By obtaining the first arrival time (i.e., the time when the wave reaches the detector) of the target original component and the target radial component, the initial P-wave velocity model and the initial S-wave velocity model can be determined.

[0124] In some embodiments, as shown in Figure 9 S201 can be specifically implemented as the following S2011-S2012.

[0125] S2011, obtaining the first arrival time of the target original component and the target radial component.

[0126] In some embodiments, the first arrival time of the down-going P-wave is extracted from the Z-component and the first arrival time of the down-going S-wave is extracted from the R-component. The first arrival time refers to the first time that the seismic wave reaches the receiver from the source. For P-wave (P-wave) and S-wave (S-wave), the first arrival time can provide basic information about the velocity of the underground medium

[0127] S2012, determining an initial velocity model based on the first arrival time.

[0128] In some embodiments, the P-wave and S-wave velocity models of the underground medium are determined by waveform inversion techniques, such as ray tracing or wave equation method, using the first arrival time of the down-going P-wave of the Z-component and the first arrival time of the down-going S-wave of the R-component.

[0129] S202, determining an initial density model according to the P-wave velocity model.

[0130] In some embodiments, the density distribution of the medium is estimated according to the P-wave velocity using an empirical formula or a rock physics model (Gardner), and the initial density model is obtained,

[0131] S203, determining a target initial model according to the initial velocity model and the initial density model.

[0132] In some embodiments, the target initial model is a complete underground medium model combining the initial velocity model and the initial density model, such as Figure 10 as shown, which provides a basis for subsequent waveform inversion.

[0133] S204, determining first wave field data according to the target initial model.

[0134] In some embodiments, the first wave field data is the seismic wave field data simulated based on the target initial model. The first wave field data is generated by seismic forward simulation, such as finite difference method or spectral element method.

[0135] S205, obtaining second wave field data.

[0136] In some embodiments, the second wave field data is the actual received seismic wave field data, i.e. the wave field data obtained based on the target original component and the target radial component.

[0137] S206, determining a residual value according to the first wave field data and the second wave field data.

[0138] In some embodiments, the residual value is the difference between the simulated wave field data and the actual wave field data, which reflects the difference between the target initial model and the actual underground medium. The residual value is determined by calculating the mean square error (MSE) between the two wave field data, and the residual value satisfies the following expression:

[0139]

[0140] wherein E(c) is a residual value, x s is a source position, x r is a receiver position, t is time, u cal (x s , x r , t) is first wavefield data, u obs (x s , x r , t) is second wavefield data.

[0141] S207, based on the residual value, iteratively update the model parameters of the target initial model.

[0142] In some embodiments, based on the residual value, an iterative algorithm such as the steepest descent method, the conjugate gradient method, the Newton method, etc. is used to continuously update the model parameters of the target initial model to reduce the residual value and improve the global search ability and convergence speed.

[0143] In some embodiments, the model vector of the nth iteration and the model vector of the n+1th iteration are obtained, and based on the model vector of the nth iteration and the model vector of the n+1th iteration, and the residual value, the step size of the nth iteration is determined.

[0144] In some embodiments, the step size of the nth iteration satisfies the following expression:

[0145] p n = m n+1 - m n ;

[0146]

[0147] wherein m n and m n+1 are the model vector of the nth iteration and the model vector of the n+1th iteration, p n is the model update vector of the n+1th iteration, q n is the gradient difference vector of the n+1th iteration, a n is the step size of the nth iteration, is the transpose of p n-1 , ||q n-1 ||2 is the two-norm of q n-1 ;

[0148] Further, based on the steepest descent method, the model parameters of the target initial model are iteratively updated according to the step size of the nth iteration and the residual value.

[0149] In some embodiments, the model vector of the n+1th iteration satisfies the following expression:

[0150]

[0151] S208, in response to the target initial model reaching the preset condition, stopping iteration.

[0152] In some embodiments, the preset condition includes at least one of the following: the number of iterations reaches a first threshold value and the residual value is less than or equal to a second threshold value, setting the upper limit of the number of iterations and the lower limit of the residual value as the condition for stopping iteration, and checking in real time during the iteration process whether the stopping condition is met, and if so, stopping iteration.

[0153] In some embodiments, as Figure 11 indicated, including the P-wave velocity model (a), the S-wave velocity model (b), and the density model (c), from Figure 11 The inversion result shows that the update amount of the shallow layer and the well model is larger, and the update amount of the deep layer is relatively small, which is related to the maximum well source distance of the shot record. At the same time, as Figure 12 indicated, the forward R (a) and Z (b) component records of the initial model and the forward R (c) and Z (d) component records of the model obtained by elastic wave full waveform inversion are included, by comparing the forward records of the initial model and the forward records of the model updated by elastic wave full waveform inversion, it can be directly seen that the reflection wave in the record of the initial model is less and not obvious, and after the elastic wave full waveform inversion, the seismic record can better approximate the seismic record of the real model. By Figure 12 comparing c, d and Figure 4 c, d, the final elastic wave forward record has a high similarity with the actual record above 1000ms; the wave group relationship on the forward record corresponds to the actual record, but the actual record still has a relatively obvious viscoelastic effect. Therefore, in order to accelerate convergence, the stopping condition of EFWI is set to reach the preset minimum residual above 1000ms.

[0154] In some embodiments, as Figure 13 indicated, the elastic wave full waveform inversion method based on the vertical seismic profile provided by the embodiments of the present application further includes the following S209:

[0155] S209, verifying the initial velocity model based on the corridor stacking method and the prestack depth migration method.

[0156] In some embodiments, based on the zero-offset VSP data and the corridor stacking method, the corridor stacking seismic profile is obtained, which is helpful for directly marking the horizon on the seismic profile and plays a bridging role between the seismic profile and the geological horizon; the prestack depth migration method is used to convert the seismic data from the time domain to the depth domain and improve the resolution and signal-to-noise ratio of the seismic image. After converting the corridor stacking profile to the depth domain, the accuracy of the velocity model is verified by comparing with the prestack depth migration profile.

[0157] In some embodiments, the inversion obtained P-SV velocity model is subjected to Walkaway VSP PP wave and PS wave Kirchhoff prestack depth migration to verify the reliability of the elastic wave full waveform inversion velocity modeling. As shown in Figure 14 (a) PP wave depth domain profile, (b) PP wave corridor, (c) PS wave depth domain profile, the PP wave prestack depth migration profile is calibrated by the corridor stacking profile, it can be seen that the wave group relationship of the two is almost consistent, and each seismic reflection event has good consistency, which shows that the P wave velocity modeling is accurate. In addition, the stratum occurrence reflected by the PP wave and PS wave prestack depth migration is close, which further shows the accuracy of the elastic wave full waveform inversion result and the reliability of the inversion process.

[0158] As shown in Figure 15 , the embodiment of the application provides a flowchart of an elastic wave full waveform inversion method based on a vertical seismic profile, first, the three-component orientation is replaced by horizontal rotation to obtain shot record of R component and T component, and the Z component and R component obtained by the three-component orientation are used for subsequent elastic wave full waveform inversion, for specific description, reference can be made to the related description of S102 above, which will not be repeated here. Further, the same noise attenuation parameter is applied to the Z component and R component, and a step-by-step wave field separation method of multi-channel vector median filtering is adopted to separate the uplink and downlink P wave and SV wave of the Z component and R component, for specific description, reference can be made to the related description of S103 above, which will not be repeated here. Then, the amplitude compensation is performed on the uplink and downlink P wave and SV wave separated from the Z component and R component, and the compensated wave field components are stacked, for specific description, reference can be made to the related description of S104 above, which will not be repeated here. Different source energies are normalized and corrected, for specific description, reference can be made to the related description of S105 above, which will not be repeated here. At the same time, the first arrival of the downlink P wave is extracted from the Z component, and the first arrival of the downlink S wave is extracted from the R component to determine the initial velocity model, and then the initial density model is obtained, and the target initial model containing the P wave velocity model, the S wave velocity model and the initial density model is determined, for specific description, reference can be made to the related description of S201-S203 above, which will not be repeated here. Then, the model parameters of the target initial model are iteratively updated through the residual value between the simulated wave field data and the actual wave field data, for specific description, reference can be made to the related description of S204-207 above, which will not be repeated here. Finally, the initial velocity model is verified based on the corridor stacking method and the prestack depth migration method, for specific description, reference can be made to the related description of S208-S209 above, which will not be repeated here.

[0159] It should be noted that the method of the embodiments of the present application can be executed by a single device, for example, a computer or a server, etc. The method of the embodiments can also be applied to a distributed scenario, and be completed by multiple devices cooperating with each other. In the case of such a distributed scenario, one of the multiple devices can only execute one or more steps in the method of the embodiments of the present application, and the multiple devices can interact with each other to complete the method.

[0160] It should be noted that some embodiments of the present application have been described above. Other embodiments are within the scope of the appended claims. In some cases, the actions or steps recited in the claims can be performed in a different order and still achieve desirable results. Additionally, the processes depicted in the figures do not necessarily require the particular order shown or sequential order in order to achieve the desired results. In certain implementations, multitasking and parallel processing can be advantageous.

[0161] Based on the same inventive concept, the present application also provides an elastic wave full waveform inversion device based on a vertical seismic profile, corresponding to any of the above-mentioned embodiment methods.

[0162] Reference Figure 16 The elastic wave full waveform inversion device based on a vertical seismic profile includes an acquisition module 1601, a processing module 1602, and a correction module 1603.

[0163] The acquisition module 1601 is configured to acquire original shot gather data, wherein the original shot gather data includes a first original component, a second original component, and a third original component obtained based on three-component orientation.

[0164] The processing module 1602 is configured to perform horizontal component rotation on the first original component and the second original component to obtain a first radial component and a tangential component.

[0165] The processing module 1602 is further configured to perform wave field separation on the third original component and the first radial component based on multi-channel vector median filtering to obtain a first downgoing wave, a second downgoing wave, a first upgoing wave, and a second upgoing wave.

[0166] The processing module 1602 is further configured to perform amplitude compensation on the first downgoing wave, the second downgoing wave, the first upgoing wave, and the second upgoing wave to obtain a fourth original component and a second radial component.

[0167] The correction module 1603 is configured to perform consistency correction on the fourth original component and the second radial component based on a root mean square value to obtain a target original component and a target radial component.

[0168] In some embodiments, the device further comprises a determining module 1604 and an updating module 1605.

[0169] The determining module 1604 is configured to determine an initial velocity model according to the target original component and the target radial component; the initial velocity comprises a P-wave velocity model and a S-wave velocity model.

[0170] The determining module 1604 is further configured to determine an initial density model according to the P-wave velocity model.

[0171] The determining module 1604 is further configured to determine a target initial model according to the initial velocity model and the initial density model.

[0172] The determining module 1604 is further configured to determine first wave field data according to the target initial model.

[0173] The obtaining module 1601 is configured to obtain second wave field data; the second wave field data is actual wave field data.

[0174] The determining module 1604 is further configured to determine a residual value according to the first wave field data and the second wave field data.

[0175] The residual value satisfies the following expression:

[0176]

[0177] Wherein, E(c) is the residual value, x s is a source position, x r is a receiver position, t is time, u cal (x s ,x r ,t) is the first wave field data, u obs (x s ,x r ,t) is the second wave field data.

[0178] The updating module 1605 is configured to iteratively update a model parameter of the target initial model based on the residual value.

[0179] The processing module 1602 is further configured to stop iteration in response to the target initial model reaching a preset condition; the preset condition comprises at least one of the following: an iteration number reaching a first threshold value and the residual value being less than or equal to a second threshold value.

[0180] In some embodiments, the determining module 1604 is specifically configured to: acquire a first arrival time of the target original component and the target radial component; and determine the initial velocity model based on the first arrival time.

[0181] In some embodiments, the updating module 1605 is specifically configured to: acquire model vectors of the n th iteration and the n+1 th iteration;

[0182] determine a step length of the n th iteration based on the model vectors of the n th iteration and the n+1 th iteration and the residual value;

[0183] The step length of the n th iteration satisfies the following expression:

[0184] p n =m n+1 -m n ;

[0185]

[0186] wherein m n and m n+1 are the model vectors of the n th iteration and the n+1 th iteration, p n is a model update vector of the n+1 th iteration, q n is a gradient difference vector of the n+1 th iteration, a n is the step length of the n th iteration, is a transpose of p n-1 , ||q n-1 ||2 is a two norm of q n-1 ;

[0187] update the model parameters of the target initial model based on the steepest gradient descent method according to the step length of the n th iteration and the residual value;

[0188] The model vector of the n+1 th iteration satisfies the following expression:

[0189]

[0190] In some embodiments, the elastic wave full waveform inversion device based on the vertical seismic profile further comprises a verifying module 1606.

[0191] The verifying module 1606 is configured to verify the initial velocity model based on a corridor stacking method and a prestack depth migration method.

[0192] In some embodiments, the third original component and the first radial component are denoised based on a same noise attenuation parameter.

[0193] In some embodiments, the amplitude compensation includes spherical spreading compensation and absorption attenuation compensation.

[0194] For the convenience of description, the above apparatus is described in various modules in terms of functions. Of course, in the implementation of the present application, the functions of the modules can be implemented in one or more software and / or hardware.

[0195] The apparatus of the above embodiments is used to implement the corresponding elastic wave full waveform inversion method based on vertical seismic profile in any of the preceding embodiments, and has the beneficial effects of the corresponding method embodiments, which will not be described here.

[0196] Based on the same inventive concept, the present application also provides an electronic device corresponding to the method of any of the above embodiments, comprising a memory, a processor and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the elastic wave full waveform inversion method based on vertical seismic profile according to any of the above embodiments.

[0197] Figure 17 A more specific hardware structure of an electronic device according to the present embodiment is shown, which can include a processor 1010, a memory 1020, an input / output interface 1030, a communication interface 1040 and a bus 1060. The processor 1010, the memory 1020, the input / output interface 1030 and the communication interface 1040 are connected to each other through the bus 1060 for communication within the device.

[0198] The processor 1010 can be implemented in the form of a general-purpose CPU (Central Processing Unit), a microprocessor, an ASIC (Application Specific Integrated Circuit) or one or more integrated circuits, etc., for executing related programs to implement the technical solutions provided by the present embodiment.

[0199] The memory 1020 can be implemented in the form of a ROM (Read Only Memory), a RAM (Random Access Memory), a static storage device, a dynamic storage device, etc. The memory 1020 can store an operating system and other application programs, and when the technical solutions provided by the present embodiment are implemented by software or firmware, the related program codes are stored in the memory 1020 and executed by the processor 1010.

[0200] The input / output interface 1030 is configured to connect an input / output module to realize information input and output. The input / output module can be configured in the device (not shown in the figure) or externally connected to the device to provide corresponding functions. The input device can include a keyboard, a mouse, a touch screen, a microphone, various sensors, etc., and the output device can include a display, a speaker, a vibrator, an indicator light, etc.

[0201] The communication interface 1040 is configured to connect a communication module (not shown in the figure) to realize communication interaction between the device and other devices. The communication module can realize communication through a wired manner (such as a USB, a network cable, etc.) or a wireless manner (such as a mobile network, WIFI, Bluetooth, etc.).

[0202] The bus 1050 includes a channel to transmit information between various components (such as the processor 1010, the memory 1020, the input / output interface 1030, and the communication interface 1040) of the device.

[0203] It should be noted that although the above device only shows the processor 1010, the memory 1020, the input / output interface 1030, the communication interface 1040, and the bus 1050, in the specific implementation process, the device can also include other components necessary for normal operation. In addition, those skilled in the art can understand that the above device can also only include components necessary for implementing the embodiments of the present specification, and does not necessarily include all the components shown in the figure.

[0204] The electronic device of the above embodiments is used to implement the corresponding elastic wave full waveform inversion method based on a vertical seismic profile in any of the preceding embodiments, and has the beneficial effects of the corresponding method embodiments, which are not described here again.

[0205] Based on the same inventive concept, corresponding to the method of any of the above embodiments, the present application also provides a non-transitory computer readable storage medium storing computer instructions for causing the computer to execute the elastic wave full waveform inversion method based on a vertical seismic profile as described in any of the above embodiments.

[0206] The computer readable medium of the embodiments can include permanent and non-permanent, removable and non-removable media, and can be implemented by any method or technology to store information. The information can be computer readable instructions, data structures, program modules or other data. Examples of computer storage media include, but are not limited to, phase change memory (PRAM), static random access memory (SRAM), dynamic random access memory (DRAM), other types of random access memory (RAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other memory technologies, compact disc read-only memory (CD-ROM), digital versatile disc (DVD) or other optical storage, magnetic cassette, magnetic tape, magnetic disk storage or other magnetic storage devices, or any other non-transmission medium that can be used to store information accessible by a computing device.

[0207] The computer instructions stored in the storage medium of the above embodiments are used to make the computer execute the elastic wave full waveform inversion method based on the vertical seismic profile as described in any of the above embodiments, and have the beneficial effects of the corresponding method embodiments, which are not repeated here.

[0208] Based on the same inventive concept, corresponding to the elastic wave full waveform inversion method based on the vertical seismic profile as described in any of the above embodiments, the present disclosure also provides a computer program product comprising a computer program. In some embodiments, the computer program is executable by one or more processors to cause the processors to perform the elastic wave full waveform inversion method based on the vertical seismic profile. The processor performing the corresponding step can belong to the corresponding execution subject corresponding to each step in each embodiment of the elastic wave full waveform inversion method based on the vertical seismic profile.

[0209] The computer program product of the above embodiments is used to make the processor execute the elastic wave full waveform inversion method based on the vertical seismic profile as described in any of the above embodiments, and has the beneficial effects of the corresponding method embodiments, which are not repeated here.

[0210] Those skilled in the art should understand that the discussion of any of the above embodiments is only exemplary and is not intended to imply that the scope of the present application (including claims) is limited to these examples; the above embodiments or technical features between different embodiments can also be combined, and there are many other changes of the aspects of the embodiments of the present application as described above. In order to be brief, they are not provided in detail.

[0211] Additionally, to simplify the description and discussion, and so as not to obscure the embodiments of the application being presented, the well-known functions or constructions of integrated circuit (IC) chips and other components can or can not be shown in the figures and will be omitted as not to unnecessarily obscure the embodiments of the application being presented. Moreover, the devices can be shown in block diagram form in order to avoid unnecessary obscurity of the present embodiments, and this also acknowledges the fact that the details in regard to the implementation of such block diagram devices are highly dependent on the platform within which the present embodiments are to be implemented (i.e., such details should be well within the purview of one of ordinary skill in the art). Where specific details are set forth in order to describe an illustrative embodiment of the application, it will be apparent to one of ordinary skill in the art that the embodiment of the application can be practiced without these specific details. In other instances, detailed descriptions of well-known methods, devices, and materials can be omitted so as not to obscure the description of the present embodiments of the application. It is intended that the specific embodiments disclosed herein are presented by way of example only and that the present application is not limited by the embodiments presented herein.

[0212] Although the present application has been described in connection with certain specific embodiments thereof, many modifications, changes, variations and substitutions will be apparent to those of ordinary skill in the art. For example, other memory architectures (e.g., dynamic RAM (DRAM)) can use the embodiments discussed.

[0213] It is therefore intended that the present application cover all such modifications, changes, variations and substitutions that fall within the broad scope of the appended claims. Accordingly, any one or more of the features, functions, structures, or other aspects of the embodiments described herein can be combined in any suitable manner to form additional embodiments, which are also within the scope of the present application. Thus, various additional embodiments of the present application are also contemplated. Therefore, the foregoing description is not intended to be limiting. Hence, any non-included alternatives, modifications, or variations should be understood as within the scope of the present application.

Claims

1. A method for elastic wave full waveform inversion based on vertical seismic profiles, characterized in that, The method includes: Obtain the raw shot gather data; the raw shot gather data includes a first raw component, a second raw component, and a third raw component obtained based on three-component orientation; The first original component and the second original component are rotated horizontally to obtain the first radial component and the tangential component. Based on multi-channel vector median filtering, the third original component and the first radial component are separated into wave fields to obtain a first down-going wave, a second down-going wave, a first up-going wave and a second up-going wave. The amplitudes of the first downwave, the second downwave, the first upwave, and the second upwave are compensated and superimposed to obtain the fourth original component and the second radial component. The root mean square values ​​of the fourth original component and the second radial component are calculated; Based on the root mean square value, the fourth original component and the second radial component are subjected to consistency correction to obtain the target original component and the target radial component. An initial velocity model is determined based on the original component of the target and the radial component of the target; the initial velocity model includes a longitudinal wave velocity model and a transverse wave velocity model. Based on the described longitudinal wave velocity model, determine the initial density model; Based on the initial velocity model and the initial density model, determine the target initial model; Based on the initial target model, determine the first wavefield data; Acquire the second wavefield data; the second wavefield data is the actual wavefield data. The residual value is determined based on the first wavefield data and the second wavefield data; The residual value satisfies the following expression: ； in, The residual value is... Location of the epicenter. For receiver location, For time, For the first wave field data, This refers to the second wavefield data; Based on the residual value, the model parameters of the target initial model are iteratively updated; The iteration stops when the target initial model reaches a preset condition; the preset condition includes at least one of the following: the number of iterations reaches a first threshold and the residual value is less than or equal to a second threshold.

2. The method according to claim 1, characterized in that, The step of determining the initial velocity model based on the original component of the target and the radial component of the target includes: Obtain the initial arrival times of the target's original component and the target's radial component; The initial velocity model is determined based on the initial arrival time.

3. The method according to claim 1, characterized in that, The step of iteratively updating the model parameters of the target initial model based on the residual value includes: Obtain the model vectors for the nth and (n+1)th iterations; Based on the model vectors of the nth and (n+1)th iterations, and the residual value, the step size of the nth iteration is determined. The step size of the nth iteration satisfies the following expression: ； ； ； in, and The model vectors for the nth and (n+1)th iterations are... This is the model update vector for the (n+1)th iteration. Let this be the gradient difference vector for the (n+1)th iteration. Let n be the step size of the nth iteration. for transpose, for The second norm; Based on the step size of the nth iteration and the residual value, the model parameters of the target initial model are iteratively updated using the steepest gradient descent method. The model vector for the (n+1)th iteration satisfies the following expression: 。 4. The method according to claim 3, characterized in that, The method further includes: The initial velocity model was verified based on the corridor overlay method and the pre-stack depth offset method.

5. The method according to claim 1, characterized in that, The third original component and the first radial component are denoised based on the same noise attenuation parameter.

6. The method according to claim 1, characterized in that, The amplitude compensation includes spherical diffusion compensation and absorption attenuation compensation.

7. An elastic wave full waveform inversion device based on a vertical seismic profile, characterized in that, The device includes: an acquisition module, a processing module, and a correction module; The acquisition module is used to acquire raw gun gathering data; the raw gun gathering data includes a first raw component, a second raw component, and a third raw component obtained based on three-component orientation; The processing module is used to perform horizontal component rotation on the first original component and the second original component to obtain a first radial component and a tangential component. The processing module is further configured to perform wavefield separation on the third original component and the first radial component based on multi-channel vector median filtering to obtain a first down-going wave, a second down-going wave, a first up-going wave and a second up-going wave; The processing module is also used to compensate the amplitudes of the first downwave, the second downwave, the first upwave, and the second upwave, and to superimpose them to obtain a fourth original component and a second radial component. The correction module is used to calculate the root mean square value of the fourth original component and the second radial component; Based on the root mean square value, the fourth original component and the second radial component are subjected to consistency correction to obtain the target original component and the target radial component. An initial velocity model is determined based on the original component of the target and the radial component of the target; the initial velocity model includes a longitudinal wave velocity model and a transverse wave velocity model. Based on the described longitudinal wave velocity model, determine the initial density model; Based on the initial velocity model and the initial density model, determine the target initial model; Based on the initial target model, determine the first wavefield data; Acquire the second wavefield data; the second wavefield data is the actual wavefield data. The residual value is determined based on the first wavefield data and the second wavefield data; The residual value satisfies the following expression: ； in, The residual value is... Location of the epicenter. For receiver location, For time, For the first wave field data, This refers to the second wavefield data; Based on the residual value, the model parameters of the target initial model are iteratively updated; The iteration stops when the target initial model reaches a preset condition; the preset condition includes at least one of the following: the number of iterations reaches a first threshold and the residual value is less than or equal to a second threshold.

8. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor, when executing the program, implements the method as claimed in any one of claims 1 to 6.

9. A non-transitory computer-readable storage medium storing computer instructions for causing a computer to perform the method of any one of claims 1 to 6.

Images