Elastic wave full waveform inversion method based on vertical seismic section

Abstract

The invention provides an elastic wave full waveform inversion method based on a vertical seismic section. The method comprises the following steps: acquiring original shot gather data; performing horizontal component rotation on the first original component and the second original component to obtain a first radial component and a tangential component; 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; performing amplitude compensation on the first downlink wave, the second downlink wave, the first uplink wave and the second uplink wave, and performing superposition 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 the root-mean-square value to obtain a target original component and a target radial component. According to the method, the accuracy and the reliability of subsequent elastic wave full waveform inversion are improved through the fine preprocessing step.

Description

Technical Field

[0001] The present application relates to the technical field of seismic exploration, and in particular to a method for full waveform inversion of elastic waves based on vertical seismic profiles. Background Art

[0002] Elastic full waveform inversion (EFWI) technology can provide high-fidelity full-wavefield information, including upgoing and downgoing shear and longitudinal waves. This feature enables the technology to achieve accurate velocity modeling, thereby improving the accuracy of depth migration.

[0003] However, EFWI has very strict requirements for data preprocessing. 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, EFWI is also very sensitive to the initial velocity model and requires an accurate initial velocity model to avoid deviations in the inversion results. Summary of the invention

[0004] In view of this, the purpose of this application is to propose an elastic wave full waveform inversion method based on vertical seismic profiles.

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

[0006] 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;

[0007] Perform horizontal component rotation on 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, the third original component and the first radial component are subjected to wave field separation to obtain a first downgoing wave, a second downgoing wave, a first upgoing wave, and a second upgoing wave;

[0009] Compensate the amplitudes of the first downgoing wave, the second downgoing wave, the first upgoing wave, and the second upgoing wave, and superimpose them to obtain a fourth original component and a second radial component;

[0010] Based on the 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 shear wave velocity model;

[0012] Determining an initial density model according to the longitudinal wave velocity model;

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

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

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

[0016] determining a residual value according to the first wavefield data and the second wavefield data;

[0017] The residual value satisfies the following expression:

[0018]

[0019] Where, E(c) is the residual value, x s is the earthquake source location, x r is the receiver position, t is the 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 updating the model parameters of the target initial model based on the residual value;

[0021] In response to the target initial model reaching a preset condition, the iteration is stopped; 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.

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

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

[0024] Based on the first arrival time, the initial velocity model is determined.

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

[0026] Get the model vectors for the nth and n+1th iterations;

[0027] Determine the step size of the nth iteration based on the model vectors of the nth and n+1th iterations, and the residual value;

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

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

[0030]

[0031]

[0032] Among them, m n and m n+1 is the model vector of the nth and n+1th iterations, p n is the model update vector for the n+1th iteration, q n is the gradient difference vector of the n+1th iteration, α n is the step size of the nth iteration, For p n-1 The transpose of , ‖q n-1 ‖ 2 for q n-1 The second norm of ;

[0033] Iteratively updating the model parameters of the target initial model based on the steepest gradient descent method according to the step size of the nth iteration and the residual value;

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

[0035]

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

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

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

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

[0040] The acquisition module is used 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 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;

[0042] The processing module is further used 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;

[0043] The processing module is further used to compensate the amplitudes of the first downgoing wave, the second downgoing wave, the first upgoing wave and the second upgoing wave, and to superpose them to obtain a fourth original component and a second radial component;

[0044] The correction module is used 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 vertical seismic profile also includes a determination module and an update module.

[0046] The determination module is used to determine an initial velocity model according to the target original component and the target radial component; the initial velocity includes a longitudinal wave velocity model and a shear wave velocity model;

[0047] The determination module is also used to determine an initial density model according to the longitudinal wave velocity model;

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

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

[0050] The acquisition module is used 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 wavefield data and the second wavefield data;

[0052] The residual value satisfies the following expression:

[0053]

[0054] Where, E(c) is the residual value, x s is the earthquake source location, x r is the receiver position, t is the 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 used to iteratively update the model parameters of the target initial model based on the residual value;

[0056] The processing module is also used to 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 and the residual value is less than or equal to a second threshold.

[0057] Optionally, the determination module is specifically used to obtain the 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 used to obtain the model vectors of the nth and n+1th iterations;

[0059] Determine the step size of the nth iteration based on the model vectors of the nth and n+1th iterations and the residual value;

[0060] The step size of the nth iteration satisfies the following expression:

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

[0062]

[0063]

[0064] Among them, m n and m n+1 is the model vector of the nth and n+1th iterations, p n is the model update vector for the n+1th iteration, q n is the gradient difference vector of the n+1th iteration, α n is the step size of the nth iteration, For p n-1 The transpose of ||q n-1 || 2 for q n-1 The second norm of ;

[0065] Iteratively updating the model parameters of the target initial model based on the steepest gradient descent method according to the step size of the nth iteration and the residual value;

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

[0067]

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

[0069] The verification module is used to verify the initial velocity model based on a corridor stacking method and a prestack depth migration method.

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

[0071] Optionally, the amplitude compensation includes spherical diffusion compensation and absorption attenuation compensation.

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

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

[0074] As can be seen from the above, the elastic wave full waveform inversion method based on vertical seismic profile provided by the present application replaces the three-component directional processing of the original shot data by rotating the horizontal component, converting the seismic wave component from the original coordinate system to the radial and tangential coordinate systems that are more convenient for analysis, which provides convenience for subsequent wave field separation and reduces the error that may be caused by coordinate system mismatch. Furthermore, the third original component and the first radial component are separated by multi-channel vector median filtering technology, which can effectively distinguish different types of wave vectors, including downgoing waves and upgoing waves, and help to more clearly identify and analyze the propagation path and characteristics of seismic waves. Further, the separated downgoing waves and upgoing waves are amplitude compensated and superimposed to obtain the fourth original component and the second radial component. The amplitude compensation can correct the amplitude attenuation caused by factors such as propagation distance and medium attenuation, and the superposition 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 corrected for consistency 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 arise during data acquisition and processing, and ensures the consistency between different components. The above-mentioned sophisticated preprocessing steps provide high-quality input data for the subsequent elastic wave full waveform inversion, which helps to reduce the impact of noise. At the same time, it retains the vector characteristics of the target original component and the target radial component, improves the accuracy and reliability of the inversion, and thus avoids deviations in the inversion results. BRIEF DESCRIPTION OF THE DRAWINGS

[0075] In order to more clearly illustrate the technical solutions in the present application or related technologies, the drawings required for use in the embodiments or related technical descriptions are briefly introduced below. Obviously, the drawings described below are only embodiments of the present application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0076] Figure 1 A schematic flow chart of a method for full waveform inversion of elastic waves based on vertical seismic profiles provided in an embodiment of the present application;

[0077] Figure 2 A schematic diagram of an original shot gather record provided in an embodiment of the present application;

[0078] Figure 3 A schematic diagram of a shot gathering record after horizontal rotation provided in an embodiment of the present application;

[0079] Figure 4 A schematic diagram of a filtered shot gather record and filtered noise provided in an embodiment of the present application;

[0080] Figure 5 A schematic diagram of an R component uplink and downlink wave field recording provided in an embodiment of the present application;

[0081] Figure 6 A schematic diagram of uplink and downlink wave field recording of a Z component provided in an embodiment of the present application;

[0082] Figure 7 A schematic diagram of source energy consistency correction provided in an embodiment of the present application;

[0083] Figure 8 A schematic flow chart of another elastic wave full waveform inversion method based on vertical seismic profile provided in an embodiment of the present application;

[0084] Fig. 9 A schematic flow chart of another elastic wave full waveform inversion method based on vertical seismic profile provided in an embodiment of the present application;

[0085] Fig.10 A schematic diagram of a target initial model provided in an embodiment of the present application;

[0086] Fig.11 A schematic diagram of an elastic wave full waveform inversion model provided in an embodiment of the present application;

[0087] Fig.12 A schematic diagram of a forward recording of an elastic wave full waveform inversion model provided in an embodiment of the present application;

[0088] Fig.13 A schematic flow chart of another elastic wave full waveform inversion method based on vertical seismic profile provided in an embodiment of the present application;

[0089] Fig.14 A schematic diagram of a pre-stack depth migration profile and corridor stack calibration provided in an embodiment of the present application;

[0090] Fig.15 A schematic flow chart of another elastic wave full waveform inversion method based on vertical seismic profile provided in an embodiment of the present application;

[0091] Fig.16 A schematic diagram of the composition of an elastic wave full waveform inversion device based on a vertical seismic profile provided in an embodiment of the present application;

[0092] Fig.17 A schematic diagram of the composition of an electronic device provided in an embodiment of the present application. DETAILED DESCRIPTION

[0093] In order to make the purpose, technical solutions and advantages of the present application more clearly understood, the present application is further described in detail below in combination with specific embodiments and with reference to the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present application, rather than all of the embodiments. Based on the embodiments in the present application, all other embodiments obtained by ordinary technicians in the field without creative work are within the scope of protection of the present application.

[0094] It should be noted that, in the embodiments of the present application, words such as "exemplarily" or "for example" are used to indicate examples, illustrations or descriptions. Any embodiment or design described as "exemplarily" or "for example" in the embodiments of the present application should not be interpreted as being more preferred or more advantageous than other embodiments or designs. Specifically, the use of words such as "exemplarily" or "for example" is intended to present related concepts in a concrete way. Unless otherwise defined, the technical terms or scientific terms used in the embodiments of the present application should have the usual meaning understood by people with ordinary skills in the field to which the present application belongs. SUMMARY OF THE INVENTION

[0096] Among the related technologies, 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, cracks and lithology changes. Compared with traditional seismic exploration methods, EFWI can better reveal the geological details of the underground, thus helping to improve the success rate of exploration and reduce development risks. However, EFWI has very strict requirements for data preprocessing. The preprocessing technology must be able to maintain the vector characteristics of the seismic wave amplitude, which means that in the processing process, not only the amplitude information of the wave must be retained, but also the polarization direction of the wave. If these vector characteristics are damaged during the preprocessing process, it will directly affect the accuracy of the inversion results.

[0097] The inventors of the present application have discovered that although the first-arrival time inversion of P-wave layer velocity based on vertical seismic profiling (VSP) technology is already quite mature, it faces challenges in processing the first-arrival of shear waves, because the first-arrival of shear waves is often masked by other wave field signals, making it difficult to extract the first-arrival of shear waves with high precision. Furthermore, noise interference is another important factor affecting the inversion stability. Many existing VSP pre-stack processing technologies find it difficult to effectively maintain the polarization relationship between the vertical component and the radial component of the elastic wave field. This defect seriously weakens the reliability and accuracy of the inversion results.

[0098] In order to solve the above problems, the present application provides a full waveform inversion method of elastic waves based on vertical seismic profiles, which replaces the three-component directional processing of the original shot data by rotating the horizontal component, converts the seismic wave components from the original coordinate system to the radial and tangential coordinate systems that are more convenient for analysis, and provides convenience for subsequent wave field separation, while reducing the errors that may be caused by coordinate system mismatch. Furthermore, the third original component and the first radial component are separated by multi-channel vector median filtering technology, which can effectively distinguish different types of wave vectors, including downgoing waves and upgoing waves, and help to more clearly identify and analyze the propagation path and characteristics of seismic waves. Further, the separated downgoing waves and upgoing waves are amplitude compensated and superimposed to obtain the fourth original component and the second radial component. The amplitude compensation can correct the amplitude attenuation caused by factors such as propagation distance and medium attenuation, and the superposition 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 corrected for consistency 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 arise during data acquisition and processing, and ensures the consistency between different components. The above-mentioned sophisticated preprocessing steps provide high-quality input data for the subsequent elastic wave full waveform inversion, which helps to reduce the impact of noise. At the same time, it retains the vector characteristics of the target original component and the target radial component, improves the accuracy and reliability of the inversion, and thus avoids deviations in the inversion results.

[0099] After introducing the basic principles of the present application, various non-limiting implementation methods of the present application are described in detail below.

[0100] Figure 1 The present invention provides a flowchart of a method for full waveform inversion of elastic waves based on vertical seismic profiles. Figure 1 As shown, the elastic wave full waveform inversion method based on vertical seismic profile provided by the present application specifically includes the following steps:

[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 data are received by the detector after being excited by the earthquake source and propagated through the underground medium. These data record information such as the amplitude, frequency, and phase of the seismic waves at different times and locations. Among them, the three-component detector can measure the vibration of the seismic waves 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 seismic waves in three mutually perpendicular directions, and the three components are respectively the X component (the first original component), the Y component (the second original component) and the Z component (the third original component).

[0105] S102, performing horizontal component rotation on the first original component and the second original component to obtain a first radial component and a tangential component.

[0106] In some embodiments, Figure 2 The X component (a) and Y component (b) of the original shot data are shown. The X component and Y component are converted through the horizontal component, as shown in Figure 3 As shown, the R component (first radial component) and the T component (tangential component) are obtained, among which the longitudinal wave (P wave) and the vertical shear wave (SV wave) are mainly projected on the R component, while the horizontal shear wave (SH wave) is mainly polarized on the T component.

[0107] It should be noted that horizontal component rotation is the process of rotating the two mutually perpendicular horizontal components (X component and Y component) in the seismic data from the observation coordinate system to the source-detector 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 in the position of the surface shot point and the fixed downhole detector, the X component and Y component of the detector both receive SH waves and SV waves. Through horizontal rotation, these wave fields can be separated into two components, R and T, making it easier to identify and analyze the characteristics of seismic waves.

[0108] S103 . Based on multi-channel vector median filtering, perform wave field separation 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, during the wave field separation process, the same noise attenuation parameter is applied to the third original component and the first radial component to maintain the polarization relationship between the amplitudes of the Z and R components of the P wave and the SV wave, such as Figure 4 As shown, including R component noise (a), Z component noise (b), R component (c) and Z component (d), it can be seen that no strong energy up and down longitudinal and transverse waves are seen in the residual record, which shows that the filtering is amplitude-preserving.

[0110] In some embodiments, wavefield separation decomposes seismic wave data into wavefield components in different directions, which is helpful for subsequent data analysis and processing. A multi-channel vector median filtering step-by-step wavefield separation method is used. According to the characteristics that the first arrival of the downgoing wave of the VSP wavefield can be well distinguished, while the phase axis of the upgoing wave is often difficult to identify, step-by-step separation is performed according to the energy of different types of wavefields from strong to weak. 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 R component wavefield records after separation are: (a) downgoing PP wave, (b) downgoing PS wave, (c) upgoing PP wave and (d) upgoing PS wave. Figure 5 As shown in Figure 2, the wave field records of the Z component after separation are: (a) downgoing PP wave, (b) downgoing PS wave, (c) upgoing PP wave and (d) upgoing PS wave. Figure 6 shown.

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

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

[0113] S104, compensating the amplitudes of the first downgoing wave, the second downgoing wave, the first upgoing wave, and the second upgoing wave, and superimposing them to obtain a fourth original component and a second radial component.

[0114] In some embodiments, amplitude compensation is performed on the up and down P waves and SV waves separated from the Z and R components, including spherical diffusion compensation and absorption attenuation compensation. In addition, for the Z and R components after wave field separation in the same mode, the amplitude compensation uses the same parameters. After amplitude compensation, the up and down P waves and SV waves have similar waveforms in different time windows to ensure that a wave field with polarization-preserving characteristics is obtained.

[0115] It should be noted that spherical diffusion compensation is a compensation method that takes into account the attenuation of seismic waves due to spherical diffusion during propagation, and absorption attenuation compensation is a compensation method that takes into account the attenuation of seismic waves due to medium absorption during propagation.

[0116] In some embodiments, an amplitude compensation factor is calculated based on the propagation distance of the seismic wave and the medium characteristics, and the amplitudes of the downgoing wave and the upgoing wave are compensated to correct the amplitude attenuation caused by the propagation distance and the medium attenuation. The compensated wave field components are superimposed, that is, the upgoing and downgoing P waves and SV waves separated in the Z component are superimposed to obtain the fourth original component that retains the vector characteristics, and the upgoing and downgoing P waves and SV waves separated in the R component are superimposed to obtain the second radial component that retains the vector characteristics.

[0117] S105 . Based on the root mean square value, perform consistency correction on the fourth original component and the second radial component to obtain a target original component and a target radial component.

[0118] In some embodiments, the root mean square (RMS) is a statistic used to measure the degree of data fluctuation or dispersion, which represents the square root of the average value of the sum of squares of the data. By calculating the root mean square value of the fourth original component and the second radial component, the average is obtained, and the Figure 7 As shown, the fourth original component and the second radial component are corrected for consistency according to the root mean square value to ensure that the energy of the fourth original component and the second radial component are 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 vertical seismic profile provided by the present application replaces the three-component directional processing of the original shot data by rotating the horizontal component, converts the seismic wave component from the original coordinate system to the radial and tangential coordinate systems that are more convenient for analysis, provides convenience for subsequent wave field separation, and reduces the error that may be caused by the mismatch of the coordinate system. Furthermore, the third original component and the first radial component are separated by the multi-channel vector median filtering technology, which can effectively distinguish different types of wave vectors, including downgoing waves and upgoing waves, and help to more clearly identify and analyze the propagation path and characteristics of seismic waves. Further, the separated downgoing waves and upgoing waves are amplitude compensated and superimposed to obtain the fourth original component and the second radial component. The amplitude compensation can correct the amplitude attenuation caused by factors such as propagation distance and medium attenuation, and the superposition 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 corrected for consistency 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 arise during data acquisition and processing, and ensures the consistency between different components. The above-mentioned sophisticated preprocessing steps provide high-quality input data for the subsequent elastic wave full waveform inversion, which helps to reduce the impact of noise. At the same time, it retains the vector characteristics of the target original component and the target radial component, improves the accuracy and reliability of the inversion, and thus avoids deviations in the inversion results.

[0120] In some embodiments, Figure 8 As shown, the elastic wave full waveform inversion method based on vertical seismic profile provided in the embodiment of the present application also includes the following S201-S208:

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

[0122] The initial velocity includes a longitudinal wave velocity model and a shear wave velocity model.

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

[0124] In some embodiments, Fig. 9 As shown, S201 can be specifically implemented as S2011-S2012 as follows:

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

[0126] In some embodiments, the first arrival time of the longitudinal wave is extracted from the downward longitudinal wave of the Z component, and the first arrival time of the shear wave is extracted from the downward shear wave of the R component. The first arrival time refers to the first time when the seismic wave arrives at the detector from the source. For the longitudinal wave (P wave) and the shear wave (S wave), the first arrival time can provide basic information about the velocity of the underground medium.

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

[0128] In some embodiments, the first arrival time of the downward P-wave of the Z component and the first arrival time of the downward S-wave of the R component are used to determine the P-wave and S-wave velocity models of the underground medium through waveform inversion techniques such as ray tracing or wave equation methods.

[0129] S202. Determine an initial density model according to the longitudinal wave velocity model.

[0130] In some embodiments, an empirical formula or a rock physics model (Gardner) is used to estimate the density distribution of the medium according to the P-wave velocity to obtain an initial density model.

[0131] S203: Determine 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 subsurface medium model that combines an initial velocity model and an initial density model, such as Fig.10 As shown, it provides a basis for the subsequent waveform inversion.

[0133] S204: Determine first wave field data according to the target initial model.

[0134] In some embodiments, the first wavefield data is seismic wavefield data obtained by simulating a target initial model. The first wavefield data is generated by seismic forward modeling, such as a finite difference method or a spectral element method.

[0135] S205: Acquire second wave field data.

[0136] In some embodiments, the second wavefield data is actually received seismic wavefield data, that is, wavefield data obtained based on a target original component and a target radial component.

[0137] S206: Determine a residual value according to the first wavefield data and the second wavefield data.

[0138] In some embodiments, the residual value is the difference between the simulated wavefield data and the actual wavefield 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 wavefield data. The residual value satisfies the following expression:

[0139]

[0140] Among them, E(c) is the residual value, x s is the earthquake source location, x r is the receiver position, t is the 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.

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

[0142] In some embodiments, based on the residual value, an iterative algorithm, such as the steepest gradient 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 capability and convergence speed.

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

[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] Among them, m n and m n+1 is the model vector for the nth and n+1th iterations, p n is the model update vector for the n+1th iteration, q n is the gradient difference vector of the n+1th iteration, α n is the step size of the nth iteration, For p n-1 The transpose of ||q n-1 || 2 for q n-1 The second norm of ;

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

[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 a preset condition, stop iteration.

[0152] In some embodiments, the preset conditions include 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, an upper limit on the number of iterations and a lower limit on the residual value are set as conditions for stopping iteration, and during the iteration process, a real-time check is made to see whether the stopping conditions are met, and if so, the iteration is stopped.

[0153] In some embodiments, Fig.11 As shown in Figure 1, it includes the longitudinal wave velocity model (a), the shear wave velocity model (b), and the density model (c). Fig.11 The inversion results show that the update amount of the shallow layer and the wellbore model is large, and the update amount of the deep layer is relatively small, which is related to the maximum well source distance recorded by the shot gather. Fig.12 As shown, it includes 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 full waveform inversion of elastic waves. By comparing the forward records of the initial model with the forward records of the updated model by full waveform inversion of elastic waves, it can be seen intuitively that the reflection waves in the records of the initial model are less and less obvious. After full waveform inversion of elastic waves, the seismic records can better approximate the seismic records of the real model. Fig.12 c, d and Figure 4 Comparing c and d, the final elastic wave forward record and the actual record have a high degree of similarity above 1000ms; the wave group relationship on the forward record basically corresponds to the actual record, but the multiple waves on the actual record still have a relatively obvious viscoelastic effect. Therefore, in order to accelerate convergence, the stop condition of EFWI is set to reach the preset minimum residual of more than 1000ms.

[0154] In some embodiments, Fig.13 As shown, the elastic wave full waveform inversion method based on vertical seismic profile provided in the embodiment of the present application also includes the following S209:

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

[0156] In some embodiments, a corridor stacked seismic profile is obtained based on zero-bias VSP data and the corridor stacking method, which helps to intuitively calibrate the layers on the seismic profile and serves as a bridge between the seismic profile and the geological layers; the prestack depth migration method is used to convert seismic data from the time domain to the depth domain and improve the resolution and signal-to-noise ratio of the seismic image. After the corridor stacked profile is converted to the depth domain, the accuracy of the velocity model is verified by comparing it with the prestack depth migration profile.

[0157] In some embodiments, the P- and S-wave layer velocity models obtained by inversion are subjected to Kirchhoff prestack depth migration of PP and PS waves of Walkaway VSP to verify the reliability of elastic wave full waveform inversion velocity modeling. Fig.14 As shown in the figure, it includes PP wave depth domain profile (a), PP wave corridor (b), and PS wave depth domain profile (c). The PP wave prestack depth migration profile is calibrated by the corridor stacking profile. It can be seen that the wave group relationship between the two is almost the same, and the seismic reflection phase axis has a good consistency, indicating that the P wave layer velocity modeling is accurate. In addition, the formation occurrence reflected by the PP wave and PS wave prestack depth migration is relatively close, which further illustrates the accuracy of the elastic wave full waveform inversion results and the reliability of the inversion process.

[0158] like Fig.15 As shown, an embodiment of the present application provides a flow chart 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 the shot collection records of the R component and the T component, and the Z component and the R component obtained by the three-component orientation are used for the subsequent elastic wave full waveform inversion. For a specific description, please refer to the relevant description of S102 above, which will not be repeated in this application. Then, the same noise attenuation parameter is applied to the Z component and the 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 waves and SV waves of the Z and R components. For a specific description, please refer to the relevant description of S103 above, which will not be repeated in this application. Then, the uplink and downlink P waves and SV waves separated from the Z and R components are amplitude compensated, and the compensated wave field components are superimposed. For a specific description, please refer to the relevant description of S104 above, which will not be repeated in this application. The co-seismic source energy is normalized and corrected. For the specific description, please refer to the relevant description of S105 above, which will not be described in detail in this application. At the same time, the initial arrival of the P-wave is extracted from the downward P-wave of the Z component and the initial arrival of the S-wave is extracted from the downward S-wave of the R component to determine the initial velocity model, and then the initial density model is obtained, and the target initial model including the P-wave velocity model, the S-wave velocity model and the initial density model is determined. For the specific description, please refer to the relevant description of S201-S203 above, which will not be described in detail in this application. Then, the model parameters of the target initial model are iteratively updated through the residual values ​​between the simulated wave field data and the actual wave field data. For the specific description, please refer to the relevant description of S204-207 above, which will not be described in detail in this application. Finally, based on the corridor stacking method and the pre-stack depth migration method, the initial velocity model is verified. For the specific description, please refer to the relevant description of S208-S209 above, which will not be described in detail in this application.

[0159] It should be noted that the method of the embodiment of the present application can be performed by a single device, such as a computer or server. The method of this embodiment can also be applied to a distributed scenario and completed by multiple devices cooperating with each other. In the case of such a distributed scenario, one of the multiple devices can only perform one or more steps in the method of the embodiment of the present application, and the multiple devices will interact with each other to complete the described method.

[0160] It should be noted that the above describes some embodiments of the present application. Other embodiments are within the scope of the appended claims. In some cases, the actions or steps recorded in the claims can be performed in an order different from that in the above embodiments and still achieve the desired results. In addition, the processes depicted in the accompanying drawings do not necessarily require the specific order or continuous order shown to achieve the desired results. In some embodiments, multitasking and parallel processing are also possible or may be advantageous.

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

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

[0163] The acquisition module 1601 is used 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;

[0164] The processing module 1602 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;

[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 used to compensate the amplitudes of the first downgoing wave, the second downgoing wave, the first upgoing wave and the second upgoing wave, and superimpose them to obtain a fourth original component and a second radial component;

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

[0168] In some embodiments, the elastic wave full waveform inversion device based on vertical seismic profile also includes a determination module 1604 and an update module 1605 .

[0169] The determination module 1604 is used to determine an initial velocity model according to the target original component and the target radial component; the initial velocity includes a longitudinal wave velocity model and a shear wave velocity model;

[0170] The determination module 1604 is also used to determine an initial density model according to the longitudinal wave velocity model;

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

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

[0173] The acquisition module 1601 is used to acquire 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 wavefield data and the second wavefield data;

[0175] The residual value satisfies the following expression:

[0176]

[0177] Where, E(c) is the residual value, x s is the earthquake source location, x r is the receiver position, t is the 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 used to iteratively update the model parameters of the target initial model based on the residual value;

[0179] The processing module 1602 is also used to 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 and the residual value is less than or equal to a second threshold.

[0180] In some embodiments, the determination module 1604 is specifically used to obtain the 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 used to obtain the model vectors of the nth and n+1th iterations;

[0182] Determine the step size of the nth iteration based on the model vectors of the nth and n+1th iterations, and the residual value;

[0183] The step size of the nth iteration satisfies the following expression:

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

[0185]

[0186] Among them, m n and m n+1 is the model vector of the nth and n+1th iterations, p n is the model update vector for the n+1th iteration, q n is the gradient difference vector of the n+1th iteration, α n is the step size of the nth iteration, For p n-1 The transpose of ||q n-1 || 2 for q n-1 The second norm of ;

[0187] Iteratively updating the model parameters of the target initial model based on the steepest gradient descent method according to the step size of the nth iteration and the residual value;

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

[0189]

[0190] In some embodiments, the elastic wave full waveform inversion device based on vertical seismic profile also includes a verification module 1606.

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

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

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

[0194] For the convenience of description, the above device is described in terms of functions divided into various modules. Of course, when implementing the present application, the functions of each module can be implemented in the same or multiple software and / or hardware.

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

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

[0197] Fig.17 A more specific schematic diagram of the hardware structure of an electronic device provided in this embodiment is shown, and the device may 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 in the device.

[0198] The processor 1010 can be implemented by a general-purpose CPU (Central Processing Unit), a microprocessor, an application-specific integrated circuit (ASIC), or one or more integrated circuits, and is used to execute relevant programs to implement the technical solutions provided in the embodiments of this specification.

[0199] The memory 1020 may be implemented in the form of ROM (Read Only Memory), RAM (Random Access Memory), static storage device, dynamic storage device, etc. The memory 1020 may store an operating system and other application programs. When the technical solutions provided in the embodiments of this specification are implemented by software or firmware, the relevant program codes are stored in the memory 1020 and are called and executed by the processor 1010.

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

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

[0202] The bus 1050 includes a path that transmits information between the various components of the device (eg, the processor 1010, the memory 1020, the input / output interface 1030, and the communication interface 1040).

[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 may also include other components necessary for normal operation. In addition, it can be understood by those skilled in the art that the above device may also only include the 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 embodiment is used to implement the corresponding elastic wave full waveform inversion method based on vertical seismic profile in any of the above embodiments, and has the beneficial effects of the corresponding method embodiment, which will not be repeated here.

[0205] Based on the same inventive concept, corresponding to any of the above-mentioned embodiment methods, the present application also provides a non-transitory computer-readable storage medium, wherein the non-transitory computer-readable storage medium stores computer instructions, and the computer instructions are used to enable the computer to execute the elastic wave full waveform inversion method based on vertical seismic profile as described in any of the above embodiments.

[0206] The computer-readable medium of this embodiment includes permanent and non-permanent, removable and non-removable media, and information storage can be implemented by any method or technology. Information can be computer-readable instructions, data structures, modules of programs, 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 technology, read-only compact disk read-only memory (CD-ROM), digital versatile disk (DVD) or other optical storage, magnetic cassettes, magnetic tape magnetic disk storage or other magnetic storage devices or any other non-transmission media that can be used to store information that can be accessed by a computing device.

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

[0208] Based on the same inventive concept, corresponding to the elastic wave full waveform inversion method based on vertical seismic profile described in any of the above embodiments, the present disclosure also provides a computer program product, which includes a computer program. In some embodiments, the computer program is executable by one or more processors so that the processor executes the elastic wave full waveform inversion method based on vertical seismic profile. Corresponding to the execution subject corresponding to each step in each embodiment of the elastic wave full waveform inversion method based on vertical seismic profile, the processor that executes the corresponding step may belong to the corresponding execution subject.

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

[0210] Those skilled in the art should understand that the discussion of any of the above embodiments is merely illustrative and is not intended to imply that the scope of the present application (including the claims) is limited to these examples. In line with the concept of the present application, the technical features in the above embodiments or different embodiments may also be combined, and there are many other variations of different aspects of the embodiments of the present application as described above, which are not provided in detail for the sake of simplicity.

[0211] In addition, to simplify the description and discussion, and in order not to make the embodiments of the present application difficult to understand, the known power supply / ground connection with the integrated circuit (IC) chip and other components may or may not be shown in the provided drawings. In addition, the device can be shown in the form of a block diagram to avoid making the embodiments of the present application difficult to understand, and this also takes into account the fact that the details of the implementation of these block diagram devices are highly dependent on the platform to be implemented in the embodiments of the present application (that is, these details should be fully within the scope of understanding of those skilled in the art). In the case of elaborating specific details (e.g., circuits) to describe exemplary embodiments of the present application, it is obvious to those skilled in the art that the embodiments of the present application can be implemented without these specific details or when these specific details are changed. Therefore, these descriptions should be considered to be illustrative rather than restrictive.

[0212] Although the present application has been described in conjunction with specific embodiments of the present application, many replacements, modifications and variations of these embodiments will be apparent to those skilled in the art from the foregoing description. For example, other memory architectures (e.g., dynamic RAM (DRAM)) may use the embodiments discussed.

[0213] The embodiments of the present application are intended to cover all such substitutions, modifications and variations that fall within the broad scope of the appended claims. Therefore, any omissions, modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the embodiments of the present application should be included in the scope of protection of the present application.

Claims

1. A method for full waveform inversion of elastic waves based on vertical seismic profiles, characterized in that: The method comprises: 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; Perform horizontal component rotation on the first original component and the second original component to obtain a first radial component and a tangential component; Based on multi-channel vector median filtering, the third original component and the first radial component are subjected to wave field separation to obtain a first downgoing wave, a second downgoing wave, a first upgoing wave, and a second upgoing wave; Compensate the amplitudes of the first downgoing wave, the second downgoing wave, the first upgoing wave, and the second upgoing wave, and superimpose them to obtain a fourth original component and a second radial component; Based on the 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.

2. The method according to claim 1, characterized in that The method further comprises: Determine an initial velocity model according to the target original component and the target radial component; the initial velocity includes a longitudinal wave velocity model and a shear wave velocity model; Determining an initial density model according to the longitudinal wave velocity model; Determining a target initial model according to the initial velocity model and the initial density model; Determining first wave field data according to the target initial model; Acquire second wave field data; the second wave field data is actual wave field data; determining a residual value according to the first wavefield data and the second wavefield data; The residual value satisfies the following expression: Where, E(c) is the residual value, x s is the earthquake source location, x r is the receiver position, t is the 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; Iteratively updating the model parameters of the target initial model based on the residual value; In response to the target initial model reaching a preset condition, the iteration is stopped; 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.

3. The method according to claim 2, characterized in that The determining of the initial velocity model according to the target original component and the target radial component comprises: Obtaining the first arrival time of the target original component and the target radial component; Based on the first arrival time, the initial velocity model is determined.

4. The method according to claim 2, characterized in that: The iterative updating of the model parameters of the target initial model based on the residual value comprises: Get the model vectors for the nth and n+1th iterations; Determine the step size of the nth iteration based on the model vectors of the nth and n+1th iterations, and the residual value; The step size of the nth iteration satisfies the following expression: p n =m n+1 -m n ; Among them, m n and m n+1 is the model vector of the nth and n+1th iterations, p n is the model update vector for the n+1th iteration, q n is the gradient difference vector of the n+1th iteration, α n is the step size of the nth iteration, For p n-1 The transpose of ||q n-1 ||2 is q n-1 The second norm of ; Iteratively updating the model parameters of the target initial model based on the steepest gradient descent method according to the step size of the nth iteration and the residual value; The model vector of the n+1th iteration satisfies the following expression:

5. The method according to claim 4, characterized in that The method further comprises: The initial velocity model is verified based on the corridor stacking method and the prestack depth migration method.

6. 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.

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

8. An elastic wave full waveform inversion device based on vertical seismic profile, characterized in that: The device comprises: an acquisition module, a processing module and a correction module; The acquisition module is used 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; 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 used 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; The processing module is further used to compensate the amplitudes of the first downgoing wave, the second downgoing wave, the first upgoing wave and the second upgoing wave, and to superpose them to obtain a fourth original component and a second radial component; The correction module is used 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.

9. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor implements the method according to any one of claims 1 to 7 when executing the program. 10 . A non-transitory computer-readable storage medium storing computer instructions, wherein the computer instructions are used to cause a computer to execute the method according to claim 1 .

Images