阅读说明：本技术 Vsp地震资料各向异性建模方法 (Method for modeling anisotropy of VSP seismic data ) 是由 左建军 杨宏伟 谷玉田 孙卫国 吕小伟 孔庆丰 苗永康 于 2019-07-09 设计创作，主要内容包括：本发明提供一种VSP地震资料各向异性建模方法,该VSP地震资料各向异性建模方法包括：步骤1,读取地震观测记录,构建地震观测系统；步骤2,建立参数初始模型；步骤3,将震源波场正向延拓,求取记录残差；步骤4,将震源波场逆向延拓,求取梯度；步骤5,求取迭代步长,更新模型参数；步骤6,判断是否满足收敛条件,在满足收敛条件时,输出模型。该VSP地震资料各向异性建模方法针对VSP地震资料进行的全波场波形反演,能够在于地面记录不同的方向上接收地震波,使得全波形反演的过程中对于陡倾角的地质构造有很好的恢复效果,最终的反演结果,能够改善复杂地质构造的成像与解释结果。(The invention provides VSP seismic data anisotropic modeling methods which comprise the steps of 1, reading seismic observation records and constructing a seismic observation system, 2, establishing a parameter initial model, 3, extending a seismic source wave field in a forward direction to obtain a record residual error, 4, extending the seismic source wave field in a reverse direction to obtain a gradient, 5, obtaining an iteration step length and updating model parameters, and 6, judging whether a convergence condition is met or not, and outputting a model when the convergence condition is met.)

The method for modeling the anisotropy of the VSP seismic data is characterized by comprising the following steps:

step 1, reading earthquake observation records and constructing an earthquake observation system;

step 2, establishing a parameter initial model;

step 3, extending the forward direction of the seismic source wave field, and solving a recording residual error;

step 4, reversely extending the seismic source wave field to obtain a gradient;

step 5, calculating iteration step length and updating model parameters;

and 6, judging whether the convergence condition is met, and outputting the model when the convergence condition is met.

2. The method for anisotropically modeling VSP seismic data according to claim 1, wherein in step 1, an observation system of the VSP seismic record is constructed by reading the position coordinates of the seismic observation record and the wave field information; the method comprises the steps of constructing a grid space required by a finite difference method in an algorithm, and using a staggered grid technology, wherein the staggered grid technology is to open two sets of mutually staggered grids simultaneously, and the two sets of mutually staggered grids represent the velocity and the stress component of particle vibration respectively.

3. The method of anisotropically modeling VSP seismic data according to claim 1, wherein in step 2, a full waveform inversion is performed, the task of the full waveform inversion being to obtain a finer model from the initial parametric model.

4. The method of claim 1 where in step 3 the numerical solution of the seismic wavefield at each time instants is solved using a finite difference solution based on existing parameters of the subsurface medium and the pseudo-longitudinal wave equation in the transversely isotropic medium.

5. The method of claim 1, wherein in step 3, while solving the wavefield numerical solution according to the observation system in step 1,

recording wave field information at the detection point; comparing the seismic simulation record obtained according to the forward continuation of the seismic source wave field with the observed seismic data to obtain a residual error; and in the process of obtaining the wave field residual error, the wave field residual error in the least square sense is used, all seismic record residual errors are summed to obtain a target function in the least square sense, and the target function is an important index for measuring the difference between the parameter model and the real underground medium in the inversion process.

The transverse isotropic medium pseudo longitudinal wave equation is as follows:

formula (1) is a pseudo longitudinal wave (qP wave) of transverse isotropic medium for forward continuation of seismic source wave field, also called as QP wave (quasi-P waves) of VTI (vertical transverse transformation isotope) medium, since the polarization characteristic of seismic wave changes in anisotropic medium, the amplitude of longitudinal wave no longer corresponds to the propagation direction , so it is defined as pseudo longitudinal wave (qP wave), in formula sigma_VAnd σ_HRepresenting vertical and horizontal stress components of a seismic wavefield; v_xAnd V_zRepresenting a horizontal velocity component and a vertical velocity component of the wavefield, respectively; ρ represents the density of the subsurface medium, V_pThe axial velocity of the VTI medium is expressed, and epsilon and delta represent Thomsen parameters in the anisotropic medium so as to characterize the directional difference in the anisotropic medium; x represents the ground level orientation coordinate and z represents the subsurface medium depth. t represents the travel time of the seismic wave. Solving by using a finite difference method, and setting an initial state, namely t is 0 and the wave field value is 0; and (4) carrying out forward continuation according to the derivation recurrence format of the formula (1).

The recursive format of finite differences is as follows:

in the formula, Δ t and Δ x represent time-discrete and space-discrete steps, respectively. The superscript 1 represents the wavefield variable at the instant of interest, and the superscript 0 represents the wavefield variable at the current instant.

6. The method of claim 1, wherein in step 4, the updating gradient is determined by performing inverse estimation of the wave field using a adjoint state method: the satisfied equation for the adjoint wavefield: m^TAnd λ ═ s, and in the backward pushing process, performing backward pushing by using a finite difference method by using the following equation:

gradient of objective function

The gradient formula is obtained based on the adjoint state method as follows:

sigma in the formula_V' and σ_H' represents the vertical and horizontal stress components of the adjoint wavefield; v_x' and V_z' represents a horizontal velocity component and a vertical velocity component of the adjoint wavefield, respectively; ρ represents the density of the subsurface medium, v represents the axial velocity of the VTI medium, and ε and δ represent anisotropy parameters; x represents the ground level orientation coordinate and z represents the subsurface medium depth. t represents the travel time of the seismic wave.

According to the gradient formula, in the process of recording residual reverse thrust, combining all components of the forward wave field of the initial model to obtain the gradient.

7. The method of claim 1, wherein in step 5, the gradient obtained in step 4 is normalized , that is, the absolute value of each value in the gradient is less than 1, the probing step size of the full waveform inversion is controlled in a range due to linear approximation, so that the whole inversion process is stable, and v is calculated according to the computed step size and the gradient normalized _p，δ，ε(v_pRepresenting the longitudinal wave velocity of the VTI medium in the axial direction, and delta and epsilon are Thomsen parameters of the anisotropic VTI medium) to obtain a new parameter model.

8. The method of anisotropic modeling of VSP seismic data according to claim 1, wherein in step 6, new subsurface medium parameters are obtained through step 5, forward continuation of the seismic wavefield is performed and the seismic record is solved, the residual is solved and new objective function values are obtained; if the residual error of the seismic record is smaller than a specified value, stopping iteration and outputting the existing parameter model; and if the seismic record residual is larger than the specified value, returning to the step 3, and continuing to iterate to make the objective function keep convergence.

Technical Field

The invention relates to the technical field of geophysics, in particular to an anisotropic modeling method for VSP seismic data.

Background

In general seismic exploration, an artificial seismic source is excited on the ground surface, seismic waves are excited by the seismic source through a method of receiving seismic data through detectors distributed on the ground surface, and then series of processing procedures such as superposition, denoising, migration and inversion are performed on the seismic data to obtain a final seismic processing result.

In the oil and gas exploration field, the anisotropy of some complex rock strata or geological structures can reach more than 50 percent, such as shale, salt dome flanks, rock formations and the like, in the data processing process, the anisotropy is considered, the inversion effect can be effectively improved, more bases are provided for reservoir prediction and seismic data explanation, the exploration precision is improved, in the case of considering the anisotropy, the VSP technology can provide seismic records in the vertical direction, compared with the conventional ground reception, the inversion effect of wave field information in different directions on the anisotropic medium is particularly obvious, novel VSP data anisotropy modeling methods are invented for solving the technical problems.

Disclosure of Invention

The invention aims to provide VSP seismic data anisotropy modeling methods for solving the processing result deviation caused by the earth medium anisotropy and providing an accurate underground parameter model aiming at the VSP seismic data.

The object of the invention can be achieved by the following technical measures: the method for modeling the anisotropy of the VSP seismic data comprises the following steps: step 1, reading earthquake observation records and constructing an earthquake observation system; step 2, establishing a parameter initial model; step 3, extending the forward direction of the seismic source wave field, and solving a recording residual error; step 4, reversely extending the seismic source wave field to obtain a gradient; step 5, calculating iteration step length and updating model parameters; and 6, judging whether the convergence condition is met, and outputting the model when the convergence condition is met.

The object of the invention can also be achieved by the following technical measures:

in the step 1, an observation system of the VSP seismic record is constructed by reading the position coordinate and wave field information of the seismic observation record; the method comprises the steps of constructing a grid space required by a finite difference method in an algorithm, and using a staggered grid technology, wherein the staggered grid technology is to open two sets of mutually staggered grids simultaneously, and the two sets of mutually staggered grids represent the velocity and the stress component of particle vibration respectively.

In step 2, a full waveform inversion is performed, which is tasked with obtaining a finer model from the initial parametric model.

In step 3, a numerical solution of the seismic wavefield at each time instants is solved using a finite difference solution based on the existing parameters of the subsurface medium and the pseudo-longitudinal wave equation in the transversely isotropic medium.

In step 3, according to the observation system in the step 1, the wave field numerical solution is solved, and meanwhile, the wave field information at the detection point is recorded; comparing the seismic simulation record obtained according to the forward continuation of the seismic source wave field with the observed seismic data to obtain a residual error; in the process of obtaining the wave field residual error, the wave field residual error under the least square meaning is used, all seismic record residual errors are summed to obtain an objective function under the least square meaning, the objective function is an important index for measuring the difference between a parameter model and a real underground medium in the inversion process,

formula (1) is a transverse isotropic medium pseudo-longitudinal wave formula for forward continuation of seismic source wave field, and is defined as pseudo-longitudinal wave because the polarization characteristic of seismic wave is changed in anisotropic medium and the amplitude of longitudinal wave is no longer consistent with propagation direction , and in the formula, sigma is_VAnd σ_HRepresenting vertical and horizontal stress components of a seismic wavefield; v_xAnd V_zRepresenting a horizontal velocity component and a vertical velocity component of the wavefield, respectively; ρ represents the density of the subsurface medium, V_pThe axial velocity of the VTI medium is expressed, and epsilon and delta represent Thomsen parameters in the anisotropic medium so as to characterize the directional difference in the anisotropic medium; x represents the ground horizontal direction coordinate, and z represents the depth of the underground medium; t represents the travel time of the seismic wave; solving by using a finite difference method, and setting an initial state, namely t is 0 and the wave field value is 0; deducing a recurrence format according to a formula (1) to carry out forward continuation;

the recursive format of finite differences is as follows:

in the formula, Δ t and Δ x respectively represent the step length of time dispersion and space dispersion; the superscript 1 represents the wavefield variable at the instant of interest, and the superscript 0 represents the wavefield variable at the current instant.

In step 4, the wave field is reversely estimated using the adjoint state method to obtain an update gradient: the satisfied equation for the adjoint wavefield: m^TAnd λ ═ s', in the backward pushing process, performing backward pushing by using a finite difference method by using the following equation:

gradient of objective function

Where M represents the equation matrix of the adjoint wavefield and T represents the transpose of the matrix; λ represents the adjoint wavefield, s' represents the source of the adjoint wavefield, typically the seismic record residual; e represents an inverted objective function; sigma 'in formula'_VAnd σ'_HRepresenting vertical and horizontal stress components of the adjoint wavefield; v_x' and V_z' represents a horizontal velocity component and a vertical velocity component of the adjoint wavefield, respectively; ρ represents the density of the subsurface medium, v represents the axial velocity of the VTI medium, and ε and δ represent anisotropy parameters; x represents the ground horizontal direction coordinate, and z represents the depth of the underground medium; t represents the travel time of the seismic wave;

the gradient formula is obtained based on the adjoint state method as follows:

sigma 'in formula'_VAnd σ'_HRepresenting vertical and horizontal stress components of the adjoint wavefield; v_x' and V_z' represents a horizontal velocity component and a vertical velocity component of the adjoint wavefield, respectively; ρ represents the density of the subsurface medium, v represents the axial velocity of the VTI medium, and ε and δ represent anisotropy parameters; x represents the ground horizontal direction coordinate, and z represents the depth of the underground medium; t represents the travel time of the seismic wave;

according to the gradient formula, in the process of recording residual reverse thrust, combining all components of the forward wave field of the initial model to obtain the gradient.

In step 5, the gradient obtained in step 4 is subjected to normalization treatment, namely the absolute value of each numerical value in the gradient is less than 1, the tentative step length of full waveform inversion is controlled within a fixed range due to linear approximation, so that the whole inversion process is kept stable, and v is subjected to normalization treatment according to the calculated step length and the gradient_pDelta, epsilon three parameters are updated, where v_pRepresenting the longitudinal wave velocity of the transverse isotropic medium in the axial direction, and delta and epsilon are Thomsen parameters of the anisotropic VTI medium, and a new parameter model is obtained.

In step 6, obtaining new underground medium parameters through step 5, carrying out forward continuation of a seismic wave field and solving a seismic record, and solving a residual error and obtaining a new objective function value; if the residual error of the seismic record is smaller than a specified value, stopping iteration and outputting the existing parameter model; and if the seismic record residual is larger than the specified value, returning to the step 3, and continuing to iterate to make the objective function keep convergence.

According to the VSP seismic data anisotropic modeling method, seismic waves can be received in different directions recorded on the ground by performing full-wave-field waveform inversion on VSP seismic data, so that a good recovery effect is achieved on a steep dip angle geological structure in the full-wave-shape inversion process, and the imaging and interpretation results of a complex geological structure can be improved by the final inversion result. The full-wave-field waveform inversion method has the advantages that the full-wave-field waveform inversion method is more accurate, the anisotropic influence of the underground medium can be overcome, VSP seismic data and ground data are processed, and an accurate underground medium parameter model is obtained.

Drawings

FIG. 1 is a flow chart of an embodiment of a method of the invention for anisotropic modeling of VSP seismic data;

FIG. 2 is a schematic view of a dimple pattern in an embodiment of the present invention;

FIG. 3 is a schematic illustration of a dimple model VSP seismic record (multiple shots) in an embodiment of the present invention;

FIG. 4 is a diagram illustrating results of full waveform inversion of a dimple model VSP in embodiment of the present invention.

Detailed Description

In order to make the aforementioned and other objects, features and advantages of the present invention comprehensible, preferred embodiments accompanied with figures are described in detail below.

As shown in fig. 1, fig. 1 is a flow chart of the method for modeling anisotropy of VSP seismic data according to the invention.

(1) And reading the seismic observation records to construct a seismic observation system.

The method comprises the steps of establishing a grid space required by a finite difference method in an algorithm, using a staggered grid technology, wherein the staggered grid technology is used for simultaneously opening two sets of mutually staggered grids which respectively represent the velocity and stress components of particle vibration, the staggered grids can effectively improve the space precision and improve the stability of the algorithm, the VSP seismic record, figure 3 represents the seismic record, is vibration information recorded by a wave detection point in a well, the vertical axis represents time, and the horizontal axis represents seismic channels, figure 3 represents multi-shot records in an anisotropic medium, and represents the change of well in-well vibration caused by a single shot recorded by pieces of information on the seismic record along with time when the horizontal axis is unchanged.

(2) And establishing a parameter initial model.

Full waveform inversion is used as inversion methods, initial model parameters are essential to the inversion process, the initial model is established by means of velocity spectrum, seismic tomography and the like, the full waveform inversion is performed by obtaining a finer model through the initial parameter model, FIG. 2 is a velocity parameter of a depression model which is common in seismic exploration, the vertical axis of the velocity parameter represents the depth from the earth surface, and the horizontal axis of the velocity parameter represents coordinates on the earth surface.

(3) And (5) carrying out forward continuation on the seismic source wave field, and solving a record residual error.

The invention relates to a forward continuation process of a seismic wave field, in particular to a seismic wave field forward evolution process, which is characterized in that a finite difference solution is used for solving the numerical value of the seismic wave field at every moments according to the existing parameters (speed parameters, anisotropic parameters and density) of an underground medium and a qP wave equation (pseudo longitudinal wave equation in a transverse isotropic medium) in a VTI medium, wave field information at a demodulation point is recorded while the numerical value of the wave field is solved according to an observation system in the step (1), seismic simulation records obtained according to the forward continuation of the seismic wave field are compared with observed seismic data, and a residual error is solved.

Equation (1) is the qP wave (quasi-P waves) equation for the transversely isotropic medium pseudo-longitudinal wave, also known as VTI (vertical Transversely isotropic) medium, used for forward continuation of the seismic source wavefield. Due to deflection of seismic waves in anisotropic mediaThe vibration characteristics are changed, and the amplitude of the longitudinal wave no longer coincides with the propagation direction , and is therefore defined as a pseudo longitudinal wave (qP wave), where σ_VAnd σ_HRepresenting vertical and horizontal stress components of a seismic wavefield; v_xAnd V_zRepresenting a horizontal velocity component and a vertical velocity component of the wavefield, respectively; ρ represents the density of the subsurface medium, V_pThe axial velocity of the VTI medium is expressed, and epsilon and delta represent Thomsen parameters in the anisotropic medium so as to characterize the directional difference in the anisotropic medium; x represents the ground level orientation coordinate and z represents the subsurface medium depth. t represents the travel time of the seismic wave. Solving by using a finite difference method, and setting an initial state, namely t is 0 and the wave field value is 0; and (4) carrying out forward continuation according to the derivation recurrence format of the formula (1).

The recursive format of finite differences is as follows:

in the formula, Δ t and Δ x represent time-discrete and space-discrete steps, respectively. The superscript 1 represents the wavefield variable at the instant of interest, and the superscript 0 represents the wavefield variable at the current instant.

(4) And (5) reversely extending the wave field of the seismic source to obtain a gradient.

Defining a least squares objective function as:

adding constraint conditions:

wherein λ is an adjoint wave field, λ ═ v'_xv'_zσ'_xσ'_z]^TAnd p is the calculated wavefield,

representing an observed seismic record. Mp-s-0 represents the anisotropy equation for the wavefield. E denotes the inverted objective function. By derivation of the adjoint wavefield λ, the adjoint wavefield equation is obtained:

thus, the satisfied equation for the adjoint wavefield: m^Tλ ═ s', i.e.:

gradient of objective function

Where M represents the equation matrix for the adjoint wavefield and T represents the transpose of the matrix. s represents the source of the wavefield. s' represents the source of the adjoint wavefield. Sigma 'in formula'_VAnd σ'_HRepresenting vertical and horizontal stress components of the adjoint wavefield; v_x' and V_z' represents a horizontal velocity component and a vertical velocity component of the adjoint wavefield, respectively; ρ represents the density of the subsurface medium, v represents the axial velocity of the VTI medium, and ε and δ represent anisotropy parameters; x represents the ground level orientation coordinate and z represents the subsurface medium depth. t represents the travel time of the seismic wave.

Sigma 'in formula'_VAnd σ'_HRepresenting vertical and horizontal stress components of the adjoint wavefield; v_x' and V_z' represents a horizontal velocity component and a vertical velocity component of the adjoint wavefield, respectively; ρ represents the density of the subsurface medium, v represents the axial velocity of the VTI medium, and ε and δ represent anisotropy parameters; x represents the ground level orientation coordinate and z represents the subsurface medium depth. t represents the travel time of the seismic wave. According to the gradient formula, the gradient can be obtained by combining all components of the forward wave field of the initial model in the process of recording the residual reverse thrust.

(5) And (5) calculating an iteration step length and updating the model parameters.

The change of the target function along with the step length can be regarded as a quadratic function, the optimal step length suitable for each iteration is solved by using a three-point method, namely the minimum value point of the target function is solved under the condition of quadratic approximation, the gradient obtained in the step (4) is subjected to treatment, namely the absolute value of each values in the gradient is less than 1, the step length of the full waveform inversion needs to be controlled within a fixed range (5 percent) due to linear approximation, so that the whole inversion process is kept stable, and v is subjected to linear inversion according to the calculated step length and the gradient which is obtained in the step _pAnd (d), updating three parameters of delta and epsilon (VTI medium axial speed and Thomsen parameters delta and epsilon) to obtain a new parameter model.

(6) Judging whether the convergence condition is satisfied, and outputting the model

And (5) acquiring new underground medium parameters, carrying out forward continuation on the seismic wave field, solving the seismic record, solving the residual error and acquiring a new objective function value. And if the residual error of the seismic record is less than a specified value, stopping iteration and outputting the existing parameter model. If the seismic record residual is greater than the specified value, iteration is continued to keep the objective function converged. FIG. 4 shows velocity parameters obtained by inversion of VSP seismic records. Since the receiving direction of the VSP seismic records is different from the ground receiving, the ideal effect is achieved on the recovery of the steep dip.