Dual-speed based full waveform inversion method and device for VTI medium
Patent Information
- Application Number
- CN202211086204.2
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2022-09-06
- Publication Date
- 2026-09-08
- Estimated Expiration
- 2042-09-06
AI Technical Summary
[0004]本发明实施例提供一种基于双速度的VTI介质全波形反演方法及设备,用以解决现有方法中参数藕合性较强,反演精度不足的问题
[0039]The dual-velocity VTI medium full waveform inversion method and apparatus provided in this invention uses initial vertical P-wave velocity, initial horizontal P-wave velocity, vertical-to-horizontal velocity transition parameters, medium density, and source wavelet as input data. It then performs forward modeling of the wave equation using the first-order velocity-stress equation under the dual-velocity mode of the VTI medium to obtain the first simulated seismic data. When the error between the first simulated seismic data and the observed seismic data does not meet the preset termination iteration condition, the error is used as input data, and the wavefield is backpropagated using the first-order velocity-stress equation under the dual-velocity mode of the VTI medium to obtain the second simulated seismic data. Based on the first and second simulated seismic data, the gradients of the vertical and horizontal P-wave velocities are determined respectively. The initial vertical and horizontal P-wave velocities are iteratively updated until the error meets the preset termination iteration condition. This yields a vertical P-wave velocity with abundant high wavenumber components and a horizontal P-wave velocity with abundant low wavenumber components, effectively reducing parameter coupling and improving inversion accuracy. Vertical P-wave velocities with abundant high wavenumber components can clearly characterize structural features and structural interfaces, and can be combined with interpretation results for reservoir prediction. Horizontal P-wave velocities with abundant low wavenumber components can be applied to seismic imaging to improve imaging accuracy and reduce well-seismic errors and exploration risks.
Smart Images

Figure CN117665937B_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of seismic data processing and interpretation technology, specifically to a method and device for full waveform inversion of transverse isotropy (VTI) media based on dual velocities. Background Technology
[0002] The pervasive anisotropy of the Earth's medium causes its elastic characteristics to vary with direction. This results in seismic waves exhibiting directional characteristics in propagation velocity, polarization direction, and amplitude attenuation, leading to phenomena such as inter-wave coupling and shear wave splitting. The study of seismic anisotropy has become a crucial topic in seismology and a challenge for both theoretical and applied seismological research. It represents a step forward in theoretical research towards understanding the wave theory of the actual Earth's medium. Conducting seismic anisotropy studies has both theoretical significance and practical value for the exploration and development of complex oil and gas reservoirs.
[0003] The conventional method for full waveform inversion in anisotropic media currently uses one velocity and two anisotropic media parameters for inversion. For example, the vertical P-wave velocity, anisotropic intensity parameter, and a transition parameter from vertical to horizontal velocity can be used for inversion. However, the coupling between the P-wave velocity and anisotropic parameters obtained by existing methods is relatively strong, and the inversion accuracy is insufficient. Summary of the Invention
[0004] This invention provides a dual-velocity VTI medium full waveform inversion method and device to solve the problem of strong parameter coupling and insufficient inversion accuracy in existing methods.
[0005] In a first aspect, embodiments of the present invention provide a method for full waveform inversion of VTI media based on dual velocity, comprising:
[0006] Step 1: Construct the initial model, which includes the initial vertical P-wave velocity, the initial horizontal P-wave velocity, the transition parameters from vertical to horizontal velocity, the medium density, and the source wavelet;
[0007] Step 2: Using the initial vertical P-wave velocity, initial horizontal P-wave velocity, vertical-to-horizontal velocity transition parameters, medium density, and source wavelet as input data, the wave equation forward modeling is performed using the first-order velocity-stress equation in the VTI medium dual-velocity mode to obtain the first simulated seismic data of the forward modeling simulation.
[0008] Step 3: Calculate the error between the first simulated seismic data and the observed seismic data. If the error meets the preset termination iteration condition, terminate the iteration and output the initial vertical P-wave velocity and the initial horizontal P-wave velocity.
[0009] Step 4: Using the error as input data, the wave field is backpropagated using the first-order velocity-stress equation in the dual-velocity mode of the VTI medium to obtain the second simulated seismic data;
[0010] Step 5: Determine the gradients of the vertical P-wave velocity and the horizontal P-wave velocity based on the first and second simulated seismic data, respectively.
[0011] Step 6: Update the initial vertical P-wave velocity according to the iteration step size and gradient of the vertical P-wave velocity to obtain the new vertical P-wave velocity; update the initial horizontal P-wave velocity according to the iteration step size and gradient of the horizontal P-wave velocity to obtain the new horizontal P-wave velocity.
[0012] Step 7: Use the new vertical P-wave velocity as the initial vertical P-wave velocity and the new horizontal P-wave velocity as the initial horizontal P-wave velocity. Repeat steps 2 to 7 until the preset termination iteration condition is met.
[0013] In one embodiment, the first-order velocity-stress equation in the dual-velocity mode of the VTI medium is determined according to the following expression:
[0014]
[0015] Among them, v p0 The vertical P-wave velocity, v, represents the velocity of a seismic wave propagating in a medium. h The horizontal P-wave velocity represents the propagation velocity of seismic waves in the medium, δ represents the transition parameter from vertical velocity to horizontal velocity, and v x ,v y ,v z Let ρ represent the displacement velocity of a medium particle in the x, y, and z directions, t represent time, f represent the source wavelet, and σ represent the displacement velocity of the medium particle in the x, y, and z directions, respectively. H σ represents the normal stress in the horizontal direction. V This represents the normal stress in the vertical direction.
[0016] In one embodiment, the first simulated seismic data includes horizontal normal stress and vertical normal stress, and the second simulated seismic data includes the conjugate of the horizontal normal stress and the conjugate of the vertical normal stress. The gradient of the vertical P-wave velocity is determined according to the following expression:
[0017]
[0018] in, This represents the conjugate of the horizontal normal stress. The conjugate of the normal stress in the vertical direction, grad vp0 This represents the gradient of the vertical longitudinal wave velocity.
[0019] In one embodiment, the gradient of the horizontal P-wave velocity is determined according to the following expression:
[0020]
[0021] Among them, grad vh This represents the gradient of the horizontal longitudinal wave velocity.
[0022] In one embodiment, the vertical P-wave velocity, the horizontal P-wave velocity, and the transition parameter from vertical to horizontal velocity satisfy the following frequency domain wave equation:
[0023]
[0024] Among them, v p0 The vertical P-wave velocity, v, represents the velocity of a seismic wave propagating in a medium. h δ represents the horizontal P-wave velocity of a seismic wave propagating in a medium, and δ represents the transition parameter from vertical velocity to horizontal velocity. Let z, x, and y represent the frequency domain wavenumbers in the z, x, and y directions, respectively, ω represent the angular frequency, and P represent the wave field.
[0025] In one embodiment, the wave field includes a background wave field and a perturbation wave field, the perturbation wave field being determined according to the following expression:
[0026]
[0027] Where s(ω) represents the frequency domain source, x s Denotes the excitation point, x r G(x) represents the receiving point. s (x, ω) represents the seismic wave originating from the excitation point x. s Green's function to point x, G(x) r (x, ω) represents the distance of the seismic wave from point x to the receiving point x. r Green's function, K vp0 Indicates v p0 Sensitive kernel function, K vh Indicates v h Sensitive kernel function, K δ The sensitive kernel function of δ, ρ represents the medium density, and v0 represents v p0 Background velocity, Δ h Δ represents the second partial derivative in the horizontal direction. v P represents the second partial derivative in the vertical direction. s This represents a perturbation wave field.
[0028] In one embodiment, the sensitive kernel function corresponding to the vertical P-wave velocity, the horizontal P-wave velocity, and the vertical-to-horizontal velocity transition parameter is determined according to the following expression:
[0029]
[0030] Among them, G s It is G(x) s The abbreviation for (x, ω) indicates that the seismic wave originates from the excitation point x. s Green's function to point x; G r It is G(x) r The abbreviation for (x, ω) represents the distance of a seismic wave from point x to the receiving point x. r Green's function.
[0031] In one embodiment, the radiation modes corresponding to the vertical P-wave velocity, horizontal P-wave velocity, and the vertical-to-horizontal velocity transition parameter are determined according to the following expression:
[0032]
[0033] Where, χ vp0 Indicates v p0 The radiation mode, χ vh Indicates v h The radiation mode, χ δ p represents the radiation mode of δ. sh p represents the sinusoidal parameter representing the seismic wave exit angle at the excitation point. rh p represents the sine parameter representing the seismic wave emission angle at the receiving point. sz p represents the cosine parameter of the seismic wave exit angle at the excitation point. rz The cosine parameter represents the seismic wave emission angle at the receiving point.
[0034] In a second aspect, embodiments of the present invention provide an electronic device, comprising:
[0035] At least one processor and memory;
[0036] The memory stores instructions that the computer executes;
[0037] At least one processor executes computer execution instructions stored in memory, causing the at least one processor to perform the dual-speed VTI medium full waveform inversion method as described in any of the first aspects.
[0038] Thirdly, embodiments of the present invention provide a computer-readable storage medium storing computer-executable instructions, which, when executed by a processor, are used to implement the dual-speed VTI medium full waveform inversion method as described in any of the first aspects.
[0039] The dual-velocity VTI medium full waveform inversion method and apparatus provided in this invention uses initial vertical P-wave velocity, initial horizontal P-wave velocity, vertical-to-horizontal velocity transition parameters, medium density, and source wavelet as input data. It then performs forward modeling of the wave equation using the first-order velocity-stress equation under the dual-velocity mode of the VTI medium to obtain the first simulated seismic data. When the error between the first simulated seismic data and the observed seismic data does not meet the preset termination iteration condition, the error is used as input data, and the wavefield is backpropagated using the first-order velocity-stress equation under the dual-velocity mode of the VTI medium to obtain the second simulated seismic data. Based on the first and second simulated seismic data, the gradients of the vertical and horizontal P-wave velocities are determined respectively. The initial vertical and horizontal P-wave velocities are iteratively updated until the error meets the preset termination iteration condition. This yields a vertical P-wave velocity with abundant high wavenumber components and a horizontal P-wave velocity with abundant low wavenumber components, effectively reducing parameter coupling and improving inversion accuracy. Vertical P-wave velocities with abundant high wavenumber components can clearly characterize structural features and structural interfaces, and can be combined with interpretation results for reservoir prediction. Horizontal P-wave velocities with abundant low wavenumber components can be applied to seismic imaging to improve imaging accuracy and reduce well-seismic errors and exploration risks. Attached Figure Description
[0040] The accompanying drawings, which are incorporated in and form part of this specification, illustrate embodiments consistent with the invention and, together with the description, serve to explain the principles of the invention.
[0041] Figure 1 A flowchart of a VTI medium full waveform inversion method based on dual velocity provided in an embodiment of the present invention;
[0042] Figure 2 A schematic diagram of the initial vertical longitudinal wave velocity provided in an embodiment of the present invention;
[0043] Figure 3 A schematic diagram of the initial horizontal longitudinal wave velocity provided in an embodiment of the present invention;
[0044] Figure 4 A schematic diagram of the actual vertical longitudinal wave velocity provided in an embodiment of the present invention;
[0045] Figure 5 A schematic diagram of the actual horizontal longitudinal wave velocity provided in an embodiment of the present invention;
[0046] Figure 6 This is a schematic diagram of the vertical longitudinal wave velocity inversion result provided in an embodiment of the present invention;
[0047] Figure 7This is a schematic diagram of the horizontal longitudinal wave velocity inversion result provided in an embodiment of the present invention;
[0048] Figure 8 A schematic diagram of a radiation mode provided in an embodiment of the present invention;
[0049] Figure 9 This is a schematic diagram of the structure of an electronic device provided in an embodiment of the present invention.
[0050] The accompanying drawings have illustrated specific embodiments of the invention, which will be described in more detail below. These drawings and descriptions are not intended to limit the scope of the invention in any way, but rather to illustrate the concept of the invention to those skilled in the art through reference to particular embodiments. Detailed Implementation
[0051] The present invention will now be described in further detail with reference to specific embodiments and accompanying drawings. Similar elements in different embodiments are referred to by associated similar element reference numerals. In the following embodiments, many details are described to facilitate a better understanding of this application. However, those skilled in the art will readily recognize that some features may be omitted in different situations, or may be replaced by other elements, materials, or methods. In some cases, certain operations related to this application are not shown or described in the specification. This is to avoid obscuring the core parts of this application with excessive description. For those skilled in the art, detailed description of these related operations is not necessary; they can fully understand the related operations based on the description in the specification and general technical knowledge in the art.
[0052] Furthermore, the features, operations, or characteristics described in the specification can be combined in any suitable manner to form various embodiments. At the same time, the steps or actions in the method description can be rearranged or adjusted in a manner obvious to those skilled in the art. Therefore, the various orders in the specification and drawings are only for the clear description of a particular embodiment and do not imply a necessary order, unless otherwise stated that a particular order must be followed.
[0053] The serial numbers assigned to components in this document, such as "first" and "second," are used only to distinguish the described objects and have no sequential or technical meaning. The terms "connection" and "linkage" used in this application, unless otherwise specified, include both direct and indirect connections (linkages).
[0054] Full waveform inversion helps assist seismic imaging and reservoir prediction. VTI medium is a common anisotropic medium, and full waveform inversion of VTI medium has important practical application value. Usually, VTI pseudo-acoustic anisotropic medium is inverted using one velocity and two anisotropic medium parameters, such as (1) vp0 ,ε,δ;(2)v h ,δ,ε;(3)v h ,η,ε;(4)v n ,η,δ;(5)v p0 , η, δ, etc., parameterized forms. Where v p0 The vertical P-wave velocity represents the velocity of the seismic wave propagating in the medium, ε represents the anisotropy intensity parameter, δ represents the transition parameter from vertical velocity to horizontal velocity, and v h The horizontal longitudinal wave velocity v represents the velocity of a seismic wave propagating in a medium. n The dynamic correction velocity is represented by η, which indicates the quantity connecting the horizontal P-wave velocity and the dynamic correction velocity. However, the anisotropic parameters in existing parameterization forms are strongly coupled, making decoupling difficult in seismic inversion. This leads to low reliability of the inversion results, hindering high-precision seismic imaging and reservoir prediction. To address these issues, this application abandons existing methods and breaks with conventional thinking by using a dual-velocity approach—emphasizing both vertical and horizontal P-wave velocities—for inversion. This reduces coupling, improves inversion accuracy, and provides more reliable data support for further seismic imaging and reservoir prediction. Specific embodiments will be described in detail below.
[0055] Example 1
[0056] Figure 1 This is a flowchart illustrating a dual-velocity VTI medium full waveform inversion method according to an embodiment of the present invention. Figure 1 As shown, the dual-velocity VTI medium full waveform inversion method provided in this embodiment may include:
[0057] S101. Construct the initial model, which includes the initial vertical P-wave velocity, the initial horizontal P-wave velocity, the transition parameters from vertical to horizontal velocity, the medium density, and the source wavelet.
[0058] In this embodiment, the initial model can be constructed, for example, using ray tomography. Alternatively, the initial model can be constructed based on observed seismic data.
[0059] S102. Using the initial vertical P-wave velocity, initial horizontal P-wave velocity, vertical-to-horizontal velocity transition parameters, medium density, and source wavelet as input data, the wave equation forward modeling is performed using the first-order velocity-stress equation in the VTI medium dual-velocity mode to obtain the first simulated seismic data of the forward modeling simulation.
[0060] The first-order velocity-stress equation in the dual-velocity mode of the VTI medium in this embodiment reflects the relationship between the vertical P-wave velocity, the horizontal P-wave velocity, and the horizontal and vertical normal stresses in the VTI medium. By performing forward modeling of the wave equation based on the first-order velocity-stress equation in the dual-velocity mode of the VTI medium, the propagating wave field can be obtained. The first simulated seismic data in the forward modeling can include the horizontal and vertical normal stresses.
[0061] S103. Calculate the error between the first simulated earthquake data and the observed earthquake data.
[0062] The full waveform inversion takes minimizing the wavefield error between observed seismic data and simulated seismic data as the objective functional. The initial given parameter model is continuously updated through parameter iteration, so that the model continuously approaches the true parameter model.
[0063] In this embodiment, after obtaining the first simulated seismic data through forward modeling using the wave equation, the error between the first simulated seismic data and the observed seismic data is calculated. The observed seismic data in this embodiment can be seismic data obtained from a real model, or it can be seismic data obtained from actual field measurements.
[0064] S104. Determine whether the error meets the preset termination iteration condition. If it does, proceed to step S105; otherwise, proceed to step S106.
[0065] The goal of full waveform inversion is to minimize error. Therefore, in this embodiment, the preset termination iteration condition can be that the error is less than a preset threshold.
[0066] S105, terminate the iteration and output the initial vertical P-wave velocity and the initial horizontal P-wave velocity.
[0067] When the error meets the preset termination iteration condition, the iteration terminates and the initial vertical P-wave velocity and initial horizontal P-wave velocity at this point are output. In this embodiment, the final obtained vertical P-wave velocity has a richer high-wavenumber velocity component, and the final obtained horizontal P-wave velocity has a richer low-wavenumber component. The two velocities correspond to high-wavenumber and low-wavenumber components, respectively, effectively reducing coupling phenomena in anisotropic parameter inversion. The effect is significant and stable, and practical. The vertical P-wave velocity with rich high-wavenumber components can clearly characterize structural features and structural interfaces, and can be combined with interpretation results for reservoir prediction. The horizontal P-wave velocity with rich low-wavenumber components can be applied to seismic imaging, improving imaging accuracy and reducing well-seismic errors and exploration risks.
[0068] S106. Using the error as input data, the wave field is backpropagated using the first-order velocity-stress equation in the dual-velocity mode of the VTI medium to obtain the second simulated seismic data.
[0069] If the error does not meet the preset termination iteration condition, the wavefield error can be used as input data, and the same operator as in the forward modeling of seismic waves can be used for wavefield backpropagation, i.e., the first-order velocity-stress equation in the dual-velocity mode of the VTI medium can be used to obtain the second simulated seismic data. The second simulated seismic data obtained by backpropagation includes the conjugate of the horizontal normal stress and the conjugate of the vertical normal stress.
[0070] S107. Determine the gradient of the vertical P-wave velocity and the gradient of the horizontal P-wave velocity based on the first simulated seismic data and the second simulated seismic data, respectively.
[0071] Based on the first simulated seismic data representing the forward propagation wavefield and the second simulated data representing the reverse propagation wavefield, the direction and rate of update of the vertical P-wave velocity can be determined, i.e., the gradient of the vertical P-wave velocity; the direction and rate of update of the horizontal P-wave velocity can also be determined, i.e., the gradient of the horizontal P-wave velocity.
[0072] S108. Update the initial vertical P-wave velocity according to the iteration step size and gradient of the vertical P-wave velocity to obtain a new vertical P-wave velocity; update the initial horizontal P-wave velocity according to the iteration step size and gradient of the horizontal P-wave velocity to obtain a new horizontal P-wave velocity.
[0073] The initial values of the iteration step size for both the vertical and horizontal P-wave velocities can be preset and can be constants. After determining the gradient and iteration step size, the iterations can be performed according to the iterative formula to obtain new vertical and horizontal P-wave velocities.
[0074] For example, the vertical P-wave velocity can be updated using the following iterative formula:
[0075]
[0076] Where, α p0 The iteration step size for the vertical P-wave velocity, grad vp0 The gradient representing the vertical P-wave velocity. This represents the initial vertical P-wave velocity during this iteration process. This represents the new vertical longitudinal wave velocity obtained during this iteration.
[0077] For example, the horizontal P-wave velocity can be updated using the following iterative formula:
[0078]
[0079] Where, α h The iteration step size, grad, represents the horizontal longitudinal wave velocity. hThe gradient representing the horizontal longitudinal wave velocity. This represents the initial horizontal P-wave velocity during this iteration process. This represents the new horizontal longitudinal wave velocity obtained during this iteration.
[0080] In this embodiment, different iteration step sizes are used for the vertical P-wave velocity and the horizontal P-wave velocity, which makes it easier to adopt matching iteration step sizes for different P-wave velocities and helps to improve the inversion accuracy.
[0081] S109. Using the new vertical P-wave velocity as the initial vertical P-wave velocity and the new horizontal P-wave velocity as the initial horizontal P-wave velocity, iterate until the preset termination iteration condition is met.
[0082] The dual-velocity VTI medium full-waveform inversion method provided in this embodiment uses the initial vertical P-wave velocity, initial horizontal P-wave velocity, vertical-to-horizontal velocity transition parameters, medium density, and source wavelet as input data. It then performs forward modeling of the wave equation using the first-order velocity-stress equation under the VTI medium dual-velocity mode to obtain the first simulated seismic data. When the error between the first simulated seismic data and the observed seismic data does not meet the preset termination iteration condition, the error is used as input data, and the wavefield is backpropagated using the first-order velocity-stress equation under the VTI medium dual-velocity mode to obtain the second simulated seismic data. Based on the first and second simulated seismic data, the gradients of the vertical and horizontal P-wave velocities are determined respectively. The initial vertical and horizontal P-wave velocities are iteratively updated until the error meets the preset termination iteration condition. This yields a vertical P-wave velocity with abundant high-wavenumber components and a horizontal P-wave velocity with abundant low-wavenumber components, effectively reducing parameter coupling and improving inversion accuracy. Vertical P-wave velocities with abundant high wavenumber components can clearly characterize structural features and structural interfaces, and can be combined with interpretation results for reservoir prediction. Horizontal P-wave velocities with abundant low wavenumber components can be applied to seismic imaging to improve imaging accuracy and reduce well-seismic errors and exploration risks.
[0083] Example 2
[0084] The wave field vector in an anisotropic medium can be defined as:
[0085] w = [v x ,v y ,v z ,σ xx ,σ yy ,σ zz ,σ xy ,σ xz ,σ yz ] T
[0086] Where w represents the wave field vector, v x ,v y ,v z σ represents the displacement velocity of the medium particle in the x, y, and z directions, respectively. xx ,σ yy ,σ zz Let σ represent the normal stresses in the x, y, and z directions, respectively. xy ,σ xz ,σ yz These represent the shear stresses in the xoy, xoz, and yoz planes, respectively.
[0087] Under the assumption of acoustic medium
[0088] w = [v x ,v y ,v z ,σ xx ,σ yy ,σ zz ] T
[0089] In VTI medium, the horizontal normal stress σ H =σ xx +σ yy Vertical normal stress σ V =σ zz .
[0090] In this mode, the first-order velocity-stress equation in the VTI medium dual-velocity mode, used for forward modeling of the wave equation and backward propagation of the wave field, can be expressed as follows:
[0091]
[0092] Among them, v p0 The vertical P-wave velocity, v, represents the velocity of a seismic wave propagating in a medium. h The horizontal P-wave velocity represents the propagation velocity of seismic waves in the medium, δ represents the transition parameter from vertical velocity to horizontal velocity, and v x ,v y ,v z Let ρ represent the displacement velocity of a medium particle in the x, y, and z directions, respectively, and let t represent time. H σ represents the normal stress in the horizontal direction. V Indicates the normal stress in the vertical direction. This indicates the partial derivative. It should be noted that when the above equation is used for forward modeling of the wave equation, f represents the source wavelet; when the above equation is used for backward propagation of the wave field, f represents the wave field error.
[0093] The first simulated seismic data includes horizontal and vertical normal stresses, while the second simulated seismic data includes the conjugate of horizontal and vertical normal stresses.
[0094] Based on the first-order velocity-stress equation in the dual-velocity mode of the VTI medium described above, the gradient of the vertical P-wave velocity can be determined using the following expression:
[0095]
[0096] in, This represents the conjugate of the horizontal normal stress. The conjugate of the normal stress in the vertical direction, grad vp0 This represents the gradient of the vertical P-wave velocity, and * indicates the conjugate of the corresponding parameter.
[0097] Based on the first-order velocity-stress equation in the dual-velocity mode of the VTI medium described above, the gradient of the horizontal P-wave velocity can be determined using the following expression:
[0098]
[0099] Among them, grad vh This represents the gradient of the horizontal longitudinal wave velocity.
[0100] Based on the first-order velocity-stress equation in the dual-velocity mode of the VTI medium described above, the gradient of the transition parameter from vertical velocity to horizontal velocity can be determined by the following expression:
[0101]
[0102] Among them, grad δ This represents the gradient of the transition parameter from vertical velocity to horizontal velocity.
[0103] It should be noted that in the gradient expression above:
[0104] The effectiveness of the method provided in this embodiment will be further illustrated below through test results. Please refer to the appendix. Figure 2-7 , Figure 2 and Figure 3 These are schematic diagrams showing the initial vertical P-wave velocity and the initial horizontal P-wave velocity, respectively. Figure 4 and Figure 5 These are schematic diagrams showing the true vertical P-wave velocity and the true horizontal P-wave velocity, respectively. Figure 6 and Figure 7 The diagrams show the vertical P-wave velocity inversion results and the horizontal P-wave velocity inversion results obtained by performing full waveform inversion using the dual-velocity VTI medium full waveform inversion method provided by this invention. Figure 2 and Figure 3 It can be done separately by... Figure 4 and Figure 5 The actual speed shown is obtained after smoothing. (This is achieved through...) Figure 2-7 It can be seen that the vertical P-wave velocity obtained by full-waveform inversion using the method provided in this invention has a richer high-wavenumber velocity component, which can clearly characterize structural features and provide a clear structural interface; the obtained horizontal P-wave velocity has a richer low-wavenumber component. During the inversion process, the two velocities correspond to high-wavenumber and low-wavenumber components respectively, reducing coupling phenomena in anisotropic parameter inversion, demonstrating significant effectiveness, stability, and practicality. The vertical P-wave velocity with rich high-wavenumber velocity components can be combined with interpretation results for reservoir prediction, while the horizontal P-wave velocity with rich low-wavenumber components can be applied to seismic imaging, improving imaging accuracy and reducing well-seismic errors and exploration risks.
[0105] Example 3
[0106] This application abandons existing methods and breaks with conventional thinking, based on dual velocities, namely, using vertical P-wave velocity and horizontal P-wave velocity for inversion. Building upon any of the above embodiments, the VTI medium full waveform inversion method based on dual velocities provided in this embodiment satisfies the following frequency domain wave equation for the vertical P-wave velocity, horizontal P-wave velocity, and the transition parameter from vertical to horizontal velocity:
[0107]
[0108] Among them, v p0 The vertical P-wave velocity, v, represents the velocity of a seismic wave propagating in a medium. h δ represents the horizontal P-wave velocity of a seismic wave propagating in a medium, and δ represents the transition parameter from vertical velocity to horizontal velocity. Let represent the wavenumbers in the frequency domains of the z, x, and y directions, respectively; ω represent the angular frequency; and P represent the wave field (including the background wave field and the perturbation wave field). The wave field P can be decomposed into the background wave field P0 and the perturbation wave field P... s That is: P = P0 + P s Under the Born approximation, the perturbation wave field P s The expression for the analytical solution is shown below:
[0109]
[0110] Where s(ω) represents the frequency domain source, x s Denotes the excitation point, x r G(x) represents the receiving point. s (x, ω) represents the seismic wave originating from the excitation point x. s Green's function to point x, G(x) r (x, ω) represents the distance of the seismic wave from point x to the receiving point x.r Green's function, K vp0 Indicates v p0 Sensitive kernel function, K vh Indicates v h Sensitive kernel function, K δ The sensitive kernel function of δ, ρ represents the medium density, and v0 represents v p0 Background velocity, Δ h Δ represents the second partial derivative in the horizontal direction. v P represents the second partial derivative in the vertical direction. s This represents a perturbation wave field.
[0111] According to the perturbation wave field P s The analytical solution yields the sensitive kernel functions corresponding to the vertical P-wave velocity, horizontal P-wave velocity, and the transition parameter from vertical to horizontal velocity. Specifically, these functions can be determined using the following expression:
[0112]
[0113] Among them, G s It is G(x) s The abbreviation for (x, ω) indicates that the seismic wave originates from the excitation point x. s Green's function to point x; G r It is G(x) r The abbreviation for (x, ω) represents the distance of a seismic wave from point x to the receiving point x. r Green's function.
[0114] According to the perturbation wave field P s The analytical solution can also yield the radiation modes corresponding to the vertical P-wave velocity, horizontal P-wave velocity, and the transition parameters from vertical to horizontal velocity. Specifically, these can be determined using the following expressions:
[0115]
[0116] Where, χ vp0 Indicates v p0 The radiation mode, χ vh Indicates v h The radiation mode, χ δ p represents the radiation mode of δ. sh p represents the sinusoidal parameter representing the seismic wave exit angle at the excitation point. rh p represents the sine parameter representing the seismic wave emission angle at the receiving point. sz p represents the cosine parameter of the seismic wave exit angle at the excitation point. rz p represents the cosine parameter of the seismic wave emission angle at the receiving point. sh =p rh =sinθ, p sz =prz =cosθ, where θ is the angle of incidence.
[0117] Please refer to Figure 8 The diagram shows a radiation pattern, where the lightest color and largest amplitude curve is χ. vp0 The curve with the darkest color is χ. vh The curve with the smallest amplitude is χ. δ The angle shown in the diagram is the angle between the incident wave and the reflected wave, which is twice the angle of incidence. Figure 8 It can be seen that the vertical longitudinal wave velocity v p0 It mainly causes disturbances to small-angle seismic data (incident angle < 45 degrees), and the horizontal P-wave velocity v h The main disturbance caused is large-angle seismic data. Compared to the two velocities, the δ parameter contributes less to the seismic data disturbance. Horizontal and vertical P-wave velocities represent large and small offset seismic data, respectively, exhibiting low coupling, which is beneficial for obtaining reliable anisotropy parameters. Therefore, performing full-waveform inversion based on horizontal and vertical P-wave velocities can not only reduce parameter coupling but also improve inversion accuracy.
[0118] Example 4
[0119] Based on any of the above embodiments, to further improve the inversion accuracy and speed, the VTI medium full waveform inversion method based on dual velocities provided in this embodiment, after determining the gradient of the vertical P-wave velocity, can also determine the update amount of the vertical P-wave velocity based on the iteration step size and the gradient of the vertical P-wave velocity, and optimize the iteration step size of the vertical P-wave velocity based on the update amount of the vertical P-wave velocity. After determining the gradient of the horizontal P-wave velocity, can also determine the update amount of the horizontal P-wave velocity based on the iteration step size and the gradient of the horizontal P-wave velocity, and optimize the iteration step size of the horizontal P-wave velocity based on the update amount of the horizontal P-wave velocity. Then, in the next iteration update process, the optimized iteration step size of the vertical P-wave velocity and the optimized iteration step size of the horizontal P-wave velocity are used for velocity update. Optionally, the LBFGS algorithm can be used to optimize the iteration step size of the vertical P-wave velocity and the iteration step size of the horizontal P-wave velocity.
[0120] Example 5
[0121] This invention also provides an electronic device, please refer to [link to relevant documentation]. Figure 9 As shown, the embodiments of the present invention are only used as examples. Figure 9 The examples are provided for illustration only and do not imply that the invention is limited to these examples. Figure 9 This is a schematic diagram of the structure of an electronic device provided according to an embodiment of the present invention. Figure 9As shown, the electronic device 20 provided in this embodiment may include: a memory 201, a processor 202, and a bus 203. The bus 203 is used to connect the various components.
[0122] The memory 201 stores a computer program, which, when executed by the processor 202, can implement the technical solutions of any of the above method embodiments.
[0123] The memory 201 and processor 202 are electrically connected directly or indirectly to enable data transmission or interaction. For example, these components can be electrically connected to each other via one or more communication buses or signal lines, such as bus 203. The memory 201 stores a computer program that implements a dual-speed VTI medium full waveform inversion method, including at least one software function module that can be stored in the memory 201 in the form of software or firmware. The processor 202 executes various functional applications and data processing by running the software program and modules stored in the memory 201.
[0124] The memory 201 may be, but is not limited to, Random Access Memory (RAM), Read Only Memory (ROM), Programmable Read-Only Memory (PROM), Erasable Programmable Read-Only Memory (EPROM), Electrically Erasable Programmable Read-Only Memory (EEPROM), etc. The memory 201 stores programs, and the processor 202 executes the programs after receiving execution instructions. Furthermore, the software programs and modules within the memory 201 may also include an operating system, which may include various software components and / or drivers for managing system tasks (such as memory management, storage device control, power management, etc.) and can communicate with various hardware or software components to provide an operating environment for other software components.
[0125] Processor 202 can be an integrated circuit chip with signal processing capabilities. The aforementioned processor 202 can be a general-purpose processor, including a Central Processing Unit (CPU), a Network Processor (NP), etc. It can implement or execute the methods, steps, and logic block diagrams disclosed in the embodiments of this invention. The general-purpose processor can be a microprocessor or any conventional processor. It is understood that... Figure 9 The structure shown is for illustrative purposes only and may include more... Figure 9 The more or fewer components shown, or having the same Figure 9 The different configurations shown. Figure 9 The components shown can be implemented in hardware and / or software.
[0126] This invention also provides a computer-readable storage medium storing a computer program thereon, which is executed by a processor to implement the technical solutions of any of the above method embodiments.
[0127] The various embodiments in this disclosure are described in a progressive manner. The same or similar parts between the various embodiments can be referred to each other. Each embodiment focuses on describing the differences from other embodiments.
[0128] The scope of protection of this disclosure is not limited to the embodiments described above. Obviously, those skilled in the art can make various modifications and variations to this disclosure without departing from its scope and spirit. If such modifications and variations fall within the scope of the claims of this disclosure and their equivalents, then the intent of this disclosure also includes such modifications and variations.
Claims
1. A method for full waveform inversion of VTI media based on dual velocity, characterized in that, include: Step 1: Construct an initial model, which includes the initial vertical P-wave velocity, the initial horizontal P-wave velocity, the transition parameters from vertical to horizontal velocity, the medium density, and the source wavelet; Step 2: Using the initial vertical P-wave velocity, the initial horizontal P-wave velocity, the vertical-to-horizontal velocity transition parameter, the medium density, and the source wavelet as input data, perform forward modeling of the wave equation using the first-order velocity-stress equation in the VTI medium dual-velocity mode to obtain the first simulated seismic data of the forward modeling simulation. Step 3: Calculate the error between the first simulated seismic data and the observed seismic data. If the error meets the preset termination iteration condition, terminate the iteration and output the initial vertical P-wave velocity and the initial horizontal P-wave velocity. Step 4: Using the error as input data, the wave field is backpropagated using the first-order velocity-stress equation in the dual-velocity mode of the VTI medium to obtain the second simulated seismic data; Step 5: Determine the gradient of the vertical P-wave velocity and the gradient of the horizontal P-wave velocity based on the first simulated seismic data and the second simulated seismic data, respectively; Step 6: Update the initial vertical P-wave velocity according to the iteration step size and gradient of the vertical P-wave velocity to obtain a new vertical P-wave velocity; update the initial horizontal P-wave velocity according to the iteration step size and gradient of the horizontal P-wave velocity to obtain a new horizontal P-wave velocity. Step 7: Using the new vertical P-wave velocity as the initial vertical P-wave velocity and the new horizontal P-wave velocity as the initial horizontal P-wave velocity, repeat steps 2 to 7 until the preset termination iteration condition is met.
2. The method according to claim 1, characterized in that, The first-order velocity-stress equation for the VTI medium dual-velocity mode is determined according to the following expression: Among them, v p0 The vertical P-wave velocity, v, represents the velocity of a seismic wave propagating in a medium. h The horizontal P-wave velocity represents the propagation velocity of seismic waves in the medium, δ represents the transition parameter from vertical velocity to horizontal velocity, and v x ,v y ,v z Let ρ represent the displacement velocity of a medium particle in the x, y, and z directions, t represent time, f represent the source wavelet, and σ represent the displacement velocity of the medium particle in the x, y, and z directions, respectively. H σ represents the normal stress in the horizontal direction. V This represents the normal stress in the vertical direction.
3. The method according to claim 2, characterized in that, The first simulated seismic data includes horizontal normal stress and vertical normal stress, and the second simulated seismic data includes the conjugate of horizontal normal stress and the conjugate of vertical normal stress. The gradient of the vertical P-wave velocity is determined according to the following expression: in, This represents the conjugate of the horizontal normal stress. The conjugate of the normal stress in the vertical direction, grad vp0 This represents the gradient of the vertical longitudinal wave velocity.
4. The method according to claim 3, characterized in that, The gradient of the horizontal longitudinal wave velocity is determined according to the following expression: Among them, grad vh This represents the gradient of the horizontal longitudinal wave velocity.
5. The method according to any one of claims 1-4, characterized in that, The vertical P-wave velocity, the horizontal P-wave velocity, and the transition parameter from vertical to horizontal velocity satisfy the following frequency domain wave equation: Among them, v p0 The vertical P-wave velocity, v, represents the velocity of a seismic wave propagating in a medium. h δ represents the horizontal P-wave velocity of a seismic wave propagating in a medium, and δ represents the transition parameter from vertical velocity to horizontal velocity. Let z, x, and y represent the frequency domain wavenumbers in the z, x, and y directions, respectively, ω represent the angular frequency, and P represent the wave field.
6. The method according to claim 5, characterized in that, The wave field includes a background wave field and a perturbation wave field, and the perturbation wave field is determined according to the following expression: Where s(ω) represents the frequency domain source, x s Denotes the excitation point, x r G(x) represents the receiving point. s (x, ω) represents the seismic wave originating from the excitation point x. s Green's function to point x, G(x) r (x, ω) represents the distance of the seismic wave from point x to the receiving point x. r Green's function, K vp0 Indicates v p0 Sensitive kernel function, K vh Indicates v h Sensitive kernel function, K δ The sensitive kernel function of δ, ρ represents the medium density, and v0 represents v p0 Background velocity, Δ h Δ represents the second partial derivative in the horizontal direction. v P represents the second partial derivative in the vertical direction. s This represents a perturbation wave field.
7. The method according to claim 6, characterized in that, The sensitive kernel functions corresponding to the vertical P-wave velocity, the horizontal P-wave velocity, and the vertical-to-horizontal velocity transition parameter are determined according to the following expression: Among them, G s It is G(x) s The abbreviation for (x, ω) indicates that the seismic wave originates from the excitation point x. s Green's function to point x; G r It is G(x) r The abbreviation for (x, ω) represents the distance of a seismic wave from point x to the receiving point x. r Green's function.
8. The method according to claim 6, characterized in that, The radiation modes corresponding to the vertical P-wave velocity, the horizontal P-wave velocity, and the vertical-to-horizontal velocity transition parameter are determined according to the following expression: Where, χ vp0 Indicates v p0 The radiation mode, χ vh Indicates v h The radiation mode, χ δ p represents the radiation mode of δ. sh p represents the sinusoidal parameter representing the seismic wave exit angle at the excitation point. rh p represents the sine parameter representing the seismic wave emission angle at the receiving point. sz p represents the cosine parameter of the seismic wave exit angle at the excitation point. rz The cosine parameter represents the seismic wave emission angle at the receiving point.
9. An electronic device, characterized in that, include: At least one processor and memory; The memory stores computer-executed instructions; The at least one processor executes computer execution instructions stored in the memory, causing the at least one processor to perform the dual-speed VTI medium full waveform inversion method as described in any one of claims 1-8.
10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores computer-executable instructions, which, when executed by a processor, are used to implement the dual-speed VTI medium full waveform inversion method as described in any one of claims 1-8.
Citation Information
Patent Citations
Seismic anisotropy parameter full waveform inversion method and device
CN103713315A
Full waveform inversion method and system based on logging data constraints
CN106324678A