CN111694051A - Gaussian beam-based viscoacoustic medium seismic wave forward modeling method - Google Patents

Abstract

The invention provides a forward modeling method for viscoacoustic medium seismic waves based on Gaussian beams, and relates to the technical field of seismic wave numerical simulation. The method first reads the forward modeling velocity model, the Q value model and the parameter files of the two models, and determines the source function; then traces the central ray in different directions from the shot point, and calculates the Gaussian beam corresponding to each ray; The window is divided into units, and the Gaussian beam corresponding to each ray is calculated after tracing the central ray from the center of the window in different directions; the ray beam pair is selected from the shot point and the center point of the window, and the reflectivity of the forward velocity model is mapped to the local plane wave ; Perform inverse tilt stacking of the Fourier transform of the local plane wave to obtain the seismic records of the receiving channels within the window; finally superimpose the seismic records corresponding to all the receiving channels to obtain the final single-shot seismic records. On the premise of ensuring the calculation accuracy, this method improves the calculation efficiency of the forward modeling of the primary scattered wave.

Description

A forward modeling method for seismic waves in viscoacoustic media based on Gaussian beams

technical field

The invention relates to the technical field of seismic wave numerical simulation, in particular to a forward modeling method for viscoacoustic medium seismic waves based on Gaussian beams.

Background technique

Seismic wavefield numerical simulation is an important method throughout the acquisition, processing and interpretation of seismic data. It is of great significance to study the fast and high-precision seismic wave forward modeling method. The seismic wave simulation method based on acoustic medium ignores the influence of the viscosity of the medium on the propagation of seismic waves, but this property will significantly affect the dynamic characteristics of seismic waves, including energy attenuation, frequency band narrowing, and phase distortion.

In the 2019 02 issue of "Acta Geophysics", Yue Yubo and others published the "Bom forward modeling of Gaussian beams of primary scattering wavefields in acoustic media", and introduced a method for the forward modeling of Gaussian beams in acoustic media. This method is based on the linearized seismic inversion theory, uses the Green's function based on the representation of Gaussian beams, and uses inverse tilt stacking to convert the local plane wave into the scattered wave field at the receiving point. The computational efficiency of the algorithm. Yue Yubo et al. simulated the complex layered model and the Marmousi model (a public seismic acoustic wave model) through the Gaussian beam Born forward modeling method of the primary scattering wave field in the acoustic medium, and the forward modeling results proved the correctness and effectiveness of the method. sex.

"Journal of Geophysics", 2017-04, published Wu Yu et al. "Seismic Forward Simulation and Reverse Time Migration in Viscoacoustic Media Based on Fractional Laplacian Decoupling", and introduced a kind of dispersion effect in the medium. A forward modeling method for seismic waves based on the fractional Laplacian viscoacoustic wave equation decoupled from attenuation effects. This method employs staggered grid finite differences to approximate time derivatives, an improved pseudospectral method to calculate spatial derivatives, and PML absorbing boundaries to remove boundary reflections. The test on the sag model shows that the seismic wave simulation method in viscoacoustic medium has good forward modeling effect.

"Petroleum Geophysical Exploration", 2019 03, published Cai Ruiqian et al.'s "Finite Difference Forward Model of Viscoacoustic Wave Equation Variable Mechanisms", and introduced a variable relaxation mechanism of the viscoacoustic wave equation proposed under the generalized standard linear body model with limited number of relaxation mechanisms The differential seismic wavefield forward modeling method uses different numbers of relaxation mechanisms and different simulation precisions in different regions of the model to achieve the purpose of unifying computational efficiency and simulation precision. And the seismic wave simulation calculation of the layered model and the SEAM model is carried out by the variable mechanical finite difference forward modeling method of the visco-acoustic wave equation, and the experimental results have obtained good results.

It can be seen from the above examples that the existing positive seismic wave forward modeling method based on Gaussian beams is aimed at acoustic media rather than viscoacoustic media. Although the finite difference forward modeling method for viscoacoustic media can obtain good seismic wave simulation results, the calculation Efficiency is also lower.

SUMMARY OF THE INVENTION

The technical problem to be solved by the present invention is to provide a forward modeling method for seismic waves in viscoacoustic medium based on Gaussian beams, aiming at the deficiencies of the above-mentioned prior art.

In order to solve the above-mentioned technical problems, the technical scheme adopted by the present invention is: a forward modeling method for viscoacoustic medium seismic waves based on Gaussian beams, comprising the following steps:

Step 1: Read in the forward modeling velocity model, the Q value model and the parameter files of the two models, and determine the source function; the parameter files of the two models include the horizontal and vertical grid numbers, grid numbers of the forward modeling velocity model and the Q value model. Spacing, shot position, main frequency of the source, reference frequency, maximum frequency, initial beam width, position of receiving channels, number of receiving channels, distance between receiving channels, sampling distance and number of samples of each seismic data;

Step 2: Trace the central ray in different directions from the shot point, and calculate the Gaussian beam corresponding to each ray;

For visco-acoustic media, after tracing the central ray from the shot point in different directions, the Gaussian beam corresponding to each ray is as follows:

Among them, u _GB (x, x _s , p _S , ω) represents the Gaussian beam corresponding to each ray after the shot point traces the central ray in different directions, x represents the grid point coordinates in the forward velocity model, and x _s is the shot point point position coordinates, p _S and A _GB (x, x _s ) are the slowness and amplitude of the Gaussian beam propagating from the shot point x _s to the grid point x, respectively, i is the imaginary unit, τ _R (x, x _s ) , τ _I (x, x _s ), τ _O (x, x _s ) represent the real travel time, imaginary travel time and Q value integral of the Gaussian beam propagating from the shot point x _s to the grid point x, respectively, ω is the angular frequency, ω _r is the reference angular frequency, Q is used to characterize the strength of the anisotropy of the forward model, and t is the travel time of the central ray;

Step 3: Divide the receiving trace of the received seismic wave by window, trace the central ray in different directions from the center of the window, and calculate the Gaussian beam corresponding to each ray;

After tracing the central ray in different directions from the center of the window, the Gaussian beam corresponding to each ray is given by the following formula:

Among them, u _GB (x, x _L , p _L , ω) is the Gaussian beam corresponding to each ray after tracing the central ray from the center of the window along different directions, x _L is the position of the center of the window, τ _R (x, x _L ) , τ _I (x, x _L ), τ _Q (x, x _L ) represent the real travel time, imaginary travel time and Q value integral of the Gaussian beam propagating from the window center x _L to the grid point x, respectively; p _L and A _GB ( x, x _L ) are the slowness and amplitude of the Gaussian beam propagating from the window center x _L to the grid point x, respectively;

Step 4: Select the ray beam pair from the shot point and the center point of the window, and use the wavelet banking method considering the Q value to map the reflectivity of the forward velocity model to a local plane wave;

The expression for mapping the reflectivity of the forward modeling velocity model to a local plane wave using the wavelet banking method considering the Q value is:

为正演速度模型中散射点的集合，δ()为雷克子波，*为褶积符号，S(ω)表示地震子波，s^H(t，τ′_I，τ′_Q)为的s(t，τ′_I，τ′_Q)的希伯特变换，τ′_R、τ′_I、τ′_Q分别为高斯束从炮点到窗中心点的实走时、虚走时以及Q值积分，w₀为射线束的初始波束宽度；Among them, U(L, x _s , p _sx , P _Lx , t) is the local plane wave of reflectivity mapping, L represents different window centers, and p _sx and p _Lx are the Gaussian beams propagated from the shot point and the window center point to the horizontal component of slowness at grid point x, c ₀ (x) and c ₁ (x) are the background velocity and perturbation velocity of the forward velocity model, respectively, _AR is the real part of the product of ray beam pair amplitudes, and A _I is the imaginary part of the product of beam pair amplitudes,

is the set of scattering points in the forward modeling velocity model, δ() is the Rake wavelet, * is the convolution symbol, S(ω) is the seismic wavelet, and s ^H (t, τ′ _I , τ′ _Q ) is the s (t, τ′ _I , τ′ _Q ) Hiebert transform, τ′ _R , τ′ _I , τ′ _Q are the real travel time, imaginary travel time and Q value integral of the Gaussian beam from the shot point to the window center point, respectively, w ₀ is the initial beam width of the ray beam;

τ′_Q对应的最大值

为中心射线上离散点记录的最大值，针对τ′_I和τ′_Q建立时间采样序列，如下公式所示：The maximum value corresponding to τ′ _I

The maximum value corresponding to τ′ _Q

is the maximum recorded at discrete points on the center ray, and a time sampling sequence is established for τ′ _I and τ′ _Q , as shown in the following formula:

分别为针对τ′_I和τ′_Q的采样间隔；Among them, N and K are the total number of sampling points for τ′ _I and τ′ _Q , respectively,

are the sampling intervals for τ′ _I and τ′ _Q , respectively;

Then s(t, τ′ _I , τ′ _Q ) is shown in the following formula at each time sampling sequence point:

进行双线性插值，进一步求得s(t，τ′_I，τ′_Q)如下公式所示：By sampling the values at the time series points

Perform bilinear interpolation to further obtain s(t, τ′ _I , τ′ _Q ) as shown in the following formula:

均为权重系数，具体表达式为：in,

are weight coefficients, and the specific expression is:

Step 5: Perform inverse tilt stacking on the Fourier transform of the local plane wave to obtain the seismic records of the receiving channel in the window, as shown in the following formula:

where u is the seismic record of the receiver trace in the window, x _r is the receiver trace in the window, p _sz and p _Lz are the vertical components of the slowness of the Gaussian beam propagating from the shot point and the center point of the window to the grid point x, respectively, ΔL is the distance between the centers of adjacent windows, U(L, x _s , p _sx , p _Lx , ω) is the Fourier transform of U(L, x _s , p _sx , p _Lx , t);

Step 6: Superimpose the seismic records corresponding to all the receiving channels to obtain the final single-shot seismic records.

The beneficial effect of adopting the above technical solution is that: a Gaussian beam-based viscoacoustic medium seismic wave forward modeling method provided by the present invention uses the Green's function represented by the Gaussian beam, and adopts the wavelet bank method considering the Q value to synthesize the local area. It can realize forward modeling of seismic waves for visco-acoustic media, and improve the calculation efficiency of the forward modeling of primary scattered waves on the premise of ensuring the calculation accuracy.

Description of drawings

1 is a flowchart of a Gaussian beam-based viscoacoustic medium seismic wave forward modeling method provided by an embodiment of the present invention;

2 is a velocity distribution diagram of a complex layered viscoacoustic model provided by an embodiment of the present invention;

3 is a Q value distribution diagram of a complex layered model provided by an embodiment of the present invention;

4 is a simulation result diagram of a complex layered viscoacoustic model finite difference method provided by an embodiment of the present invention;

FIG. 5 is a simulation result diagram of a Gaussian beam Bon forward method for a primary scattering field of a complex layered viscoacoustic model provided for an embodiment of the present invention.

Detailed ways

The specific embodiments of the present invention will be described in further detail below with reference to the accompanying drawings and embodiments. The following examples are intended to illustrate the present invention, but not to limit the scope of the present invention.

This embodiment is based on a complex layered viscoacoustic model, and uses the Gaussian beam-based viscoacoustic medium seismic wave forward modeling method of the present invention to perform forward modeling simulation of seismic waves in a viscoacoustic medium.

In this embodiment, a Gaussian beam-based seismic wave forward modeling method for viscoacoustic medium, as shown in FIG. 1 , includes the following steps:

Step 1: Read in the forward modeling velocity model, the Q value model and the parameter files of the two models, and determine the source function; the parameter files of the two models include the horizontal and vertical grid numbers, grid spacing, The main frequency of the source of the shot location, the reference frequency ω _r , the maximum frequency, the initial beam width w ₀ of the ray beam, the position of the receiving channel, the number of receiving channels, the distance between the receiving channels, the sampling interval of each seismic data and the number of samples;

The velocity distribution and Q value distribution of the complex layered viscoacoustic model in this embodiment are shown in Figures 2 and 3. In the figures, x represents the lateral distance, and z represents the depth. There are 1000 grid points in the lateral direction of the model. The distance is 10 meters; there are 550 grid points in the longitudinal direction, and the grid spacing is 5 meters; the shot point is located in the middle of the top layer of the model; a total of 241 receiving channels, with a channel spacing of 20 meters, are evenly on both sides of the shot point.

Step 2: Trace the central ray in different directions from the shot point, and calculate the Gaussian beam corresponding to each ray;

For visco-acoustic media, after tracing the central ray from the shot point in different directions, the Gaussian beam corresponding to each ray is as follows:

Among them, u _GB (x, x _s , p _S , ω) represents the Gaussian beam corresponding to each ray after the shot point traces the central ray in different directions, x represents the grid point coordinates in the forward velocity model, and x _s is the shot point Point position coordinates, p _S and A _GB (x, x _s ) are the slowness and amplitude of Gaussian beam propagation from shot point x _s to grid point x, respectively, i is an imaginary unit, τ _R (x, x _s ), τ _I (x, x _s ), τ _Q (x, x _s ) represent the real travel time, imaginary travel time and Q value integral of the Gaussian beam propagating from the shot point x _s to the grid point x, respectively, ω is the angular frequency, ω _r is the reference angular frequency, Q is used to characterize the strength of the anisotropy of the forward velocity model, and t is the travel time of the central ray;

Step 3: Divide the receiving trace of the received seismic wave by window, trace the central ray in different directions from the center of the window, and calculate the Gaussian beam corresponding to each ray;

After tracing the central ray in different directions from the center of the window, the Gaussian beam corresponding to each ray is given by the following formula:

Among them, u _GB (x, x _L , p _L , ω) is the Gaussian beam corresponding to each ray after tracing the central ray from the center of the window in different directions, x _L is the position of the center of the window, p _L and A _GB (x, x _L ) are the slowness and amplitude of the Gaussian beam propagating from the window center x _L to the grid point x, respectively, τ _R (x, x _L ), τ _I (x, x _L ), τ _Q (x, x _L ) represent the real travel time, imaginary travel time and Q value integral of the Gaussian beam propagating from the window center x _L to the grid point x, respectively;

Step 4: Select the ray beam pair from the shot point and the center point of the window, and use the wavelet banking method considering the Q value to map the reflectivity of the forward velocity model to a local plane wave;

The expression for mapping the reflectivity of the forward modeling velocity model to a local plane wave using the wavelet banking method considering the Q value is:

为正演速度模型中散射点的集合，δ()为雷克子波，*为褶积符号，S(ω)表示地震子波，s^H(t，τ′_I，τ′_Q)为的s(t，τ′_I，τ′_Q)的希伯特变换，τ′_R、τ′_I、τ′_Q分别为高斯束从炮点到窗中心点的实走时、虚走时以及Q值积分，w₀为射线束的初始波束宽度；Among them, U(L, x _s , p _sx , P _Lx , t) is the local plane wave of reflectivity mapping, L represents different window centers, and p _sx and p _Lx are the Gaussian beams propagated from the shot point and the window center point to the horizontal component of slowness at grid point x, c ₀ (x) and c ₁ (x) are the background velocity and perturbation velocity of the forward velocity model, respectively, _AR is the real part of the product of ray beam pair amplitudes, and A _I is the imaginary part of the product of beam pair amplitudes,

is the set of scattering points in the forward modeling velocity model, δ() is the Rake wavelet, * is the convolution symbol, S(ω) is the seismic wavelet, and s ^H (t, τ′ _I , τ′ _Q ) is the s (t, τ′ _I , τ′ _Q ) Hiebert transform, τ′ _R , τ′ _I , τ′ _Q are the real travel time, imaginary travel time and Q value integral of the Gaussian beam from the shot point to the window center point, respectively, w ₀ is the initial beam width of the ray beam;

τ′_Q对应的最大值

为中心射线上离散点记录的最大值，针对τ′_I和τ′_Q建立时间采样序列，如下公式所示：The maximum value corresponding to τ′ _I

The maximum value corresponding to τ′ _Q

is the maximum recorded at discrete points on the center ray, and a time sampling sequence is established for τ′ _I and τ′ _Q , as shown in the following formula:

分别为针对τ′_I和τ′_Q的采样间隔；Among them, N and K are the total number of sampling points for τ′ _I and τ′ _Q , respectively,

are the sampling intervals for τ′ _I and τ′ _Q , respectively;

Then s(t, τ′ _I , τ′ _Q ) is shown in the following formula at each time sampling sequence point:

进行双线性插值，进一步求得s(t，τ′_I，τ′_Q)如下公式所示：By sampling the values at the time series points

Perform bilinear interpolation to further obtain s(t, τ′ _I , τ′ _Q ) as shown in the following formula:

均为权重系数，具体表达式为：in,

are weight coefficients, and the specific expression is:

Step 5: Perform inverse tilt stacking on the Fourier transform of the local plane wave to obtain the seismic records of the receiving channel in the window, as shown in the following formula:

where u(x _r , x _s , ω) is the seismic record of the receiver trace in the window, x _r is the receiver trace in the window, p _sz , p _Lz are the Gaussian beams propagated from the shot point and the center point of the window to the grid, respectively Vertical component of slowness at point x, ΔL distance between adjacent window centers, U(L, x _s , p _sx , p _Lx , ω) is the Fu of U(L, x _s , p _sx , P _Lx , t) Lie transform;

Step 6: Superimpose the seismic records corresponding to all the receiving channels to obtain the final single-shot seismic records.

In this embodiment, based on the complex layered visco-acoustic model, the seismic wave simulation result using the finite difference method is shown in Figure 4, while the seismic wave simulation result using the Gaussian beam-based viscoacoustic medium seismic wave forward modeling method of the present invention is shown in Figure 5 As shown in the two figures, offset is the offset distance, which is used to represent the distance between the receiving channel and the shot point, and T is the receiving time, which represents the time it takes for the receiving channel to receive the seismic wave. It can be seen from the two forward modeling results that the primary scattering seismic wavefield obtained by the method of the present invention and the results obtained by the finite difference method are almost the same, but under the same calculation conditions, the calculation time used by the finite difference method is 321 seconds, while the method of the present invention only takes 5.3 seconds, which improves the computing efficiency by nearly 60 times.

Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, but not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that it can still be The technical solutions described in the foregoing embodiments are modified, or some or all of the technical features thereof are equivalently replaced; and these modifications or replacements do not make the essence of the corresponding technical solutions depart from the scope defined by the claims of the present invention.

Claims (6)

1. a viscoacoustic medium seismic wave forward modeling method based on Gaussian beam, is characterized in that: comprise the following steps: Step 1: Read in the forward modeling velocity model, the Q value model and the parameter files of the two models, and determine the source function; Step 2: Trace the central ray in different directions from the shot point, and calculate the Gaussian beam corresponding to each ray; Step 3: Divide the receiving trace of the received seismic wave by window, trace the central ray in different directions from the center of the window, and calculate the Gaussian beam corresponding to each ray; Step 4: Select the ray beam pair from the shot point and the center point of the window, and use the wavelet banking method considering the Q value to map the reflectivity of the forward velocity model to a local plane wave; Step 5: Perform inverse tilt stacking on the Fourier transform of the local plane wave to obtain the seismic records of the receiving channel in the window; Step 6: Superimpose the seismic records corresponding to all the receiving channels to obtain the final single-shot seismic records. 2. a kind of viscoacoustic medium seismic wave forward modeling method based on Gaussian beam according to claim 1, is characterized in that: the parameter files of the two models described in step 1 comprise the horizontal and vertical grids of forward modeling velocity model and Q value model number, grid spacing, shot location, main frequency of the source, reference frequency, maximum frequency, initial beam width, position of receiving channels, number of receiving channels, distance between receiving channels, sampling interval and number of samples of each seismic data. 3. a kind of viscoacoustic medium seismic wave forward modeling method based on Gaussian beam according to claim 2, it is characterized in that: for viscoacoustic medium, after the center ray is traced along different directions from shot point, the Gaussian beam corresponding to each ray As shown in the following formula:

Among them, u _GB (x, x _s , p _S , ω) represents the Gaussian beam corresponding to each ray after the shot point traces the central ray in different directions, x represents the grid point coordinates in the forward velocity model, and x _s is the shot point point position coordinates, p _S and A _GB (x, x _s ) are the slowness and amplitude of the Gaussian beam propagating from the shot point x _s to the grid point x, respectively, i is the imaginary unit, τ _R (x, x _s ) , τ _I (x, x _s ), τ _Q (x, x _s ) represent the real travel time, imaginary travel time and Q value integral of the Gaussian beam propagating from the shot point x _s to the grid point x, respectively, ω is the angular frequency, ω _r is the reference angular frequency, Q is used to characterize the strength of the anisotropy of the forward model, and t is the travel time of the central ray. 4. A method for seismic wave forward modeling of viscoacoustic medium based on Gaussian beam according to claim 3, characterized in that: after tracing the central ray from the center of the window along different directions in step 3, the Gaussian beam corresponding to each ray As shown in the following formula:

Among them, u _GB (x, x _L , p _L , ω) is the Gaussian beam corresponding to each ray after tracing the central ray from the center of the window along different directions, x _L is the position of the center of the window, τ _R (x, x _L ) , τ _I (x, x _L ), τ _Q (x, x _L ) represent the real travel time, imaginary travel time and Q value integral of the Gaussian beam propagating from the window center x _L to the grid point x, respectively; p _L and A _GB ( x, x _L ) are the slowness and amplitude, respectively, of the Gaussian beam propagating from the window center x _L to the grid point x. 5. a kind of viscoacoustic medium seismic wave forward modeling method based on Gaussian beam according to claim 4, is characterized in that: described in step 4, use the wavelet bank method that considers Q value to map the reflectivity of forward modeling velocity model as: The expression for a local plane wave is:

where U(L, x _s , p _sx , p _Lx , t) is the local plane wave of reflectivity mapping, L represents different window centers, and p _sx and p _Lx are the Gaussian beams propagating from the shot point and the window center point to the horizontal component of slowness at grid point x, c ₀ (x) and c ₁ (x) are the background velocity and perturbation velocity of the forward velocity model, respectively, _AR is the real part of the product of ray beam pair amplitudes, and A _I is the imaginary part of the product of beam pair amplitudes,

is the set of scattering points in the forward modeling velocity model, δ() is the Rake wavelet, * is the convolution symbol, S(ω) is the seismic wavelet, and s ^H (t, τ′ _I , τ′ _Q ) is the s (t, τ′ _I , τ′ _Q ) Hiebert transform, τ′ _R , τ′ _I , τ′ _Q are the real travel time, imaginary travel time and Q value integral of the Gaussian beam from the shot point to the window center point, respectively, w ₀ is the initial beam width of the ray beam; The maximum value corresponding to τ′ _I

The maximum value corresponding to τ′ _Q

is the maximum recorded at discrete points on the center ray, and a time sampling sequence is established for τ′ _I and τ′ _Q , as shown in the following formula:

Among them, N and K are the total number of sampling points for τ′ _I and τ′ _Q , respectively,

are the sampling intervals for τ′ _I and τ′ _Q , respectively; Then s(t, τ′ _I , τ′ _Q ) is shown in the following formula at each time sampling sequence point:

By sampling the values at the time series points

Perform bilinear interpolation to further obtain s(t, τ′ _I , τ′ _Q ) as shown in the following formula:

in,

are weight coefficients, and the specific expression is:

6. A kind of viscoacoustic medium seismic wave forward modeling method based on Gaussian beam according to claim 5, is characterized in that: described in step 5, the Fourier transform of the local plane wave is carried out inverse tilt stacking to obtain the earthquake of the receiving channel in the window The record is shown in the following formula:

where u is the seismic record of the receiver trace in the window, x _r is the receiver trace in the window, p _sz and p _Lz are the vertical components of the slowness of the Gaussian beam propagating from the shot point and the center point of the window to the grid point x, respectively, ΔL is the distance between adjacent window centers, and U(L, x _s , p _sx , p _Lx , ω) is the Fourier transform of U(L, x _s , p _sx , p _Lx , t).

Images