WO2016008105A1

WO2016008105A1 - Post-stack wave impedance inversion method based on cauchy distribution

Info

Publication number: WO2016008105A1
Application number: PCT/CN2014/082265
Authority: WO
Inventors: 杨顺伟
Original assignee: 杨顺伟
Priority date: 2014-07-15
Filing date: 2014-07-15
Publication date: 2016-01-21
Also published as: CN104769458A

Abstract

A post-stack wave impedance inversion method based on a Cauchy distribution, comprising: employing a conventional seismic exploration method to collect seismic data, and conducting conventional processing to obtain post-stack seismic data; conducting horizon picking on the post-stack seismic data to obtain horizon data, and checking, correcting, interpolating and smoothing the horizon of a determined target layer; employing a conventional logging method to obtain the time difference curve and density curve of a logged sound wave; according to the post-stack seismic data, horizon data and known well stratification data, calibrating a time difference curve and density curve of the sound wave of a depth domain as a curve of a time domain, simultaneously generating wave impedance curve data in the time domain in a well, and extracting a seismic wavelet; generating an initial wave impedance model; respectively normalizing wavelet data; acquiring a reflection coefficient sequence of a channel via inversion; deducing, via a definition of a relative wave impedance at an n-th sampling point, to obtain a relationship between the wave impedance and the reflection coefficient, and repeating the above process on all seismic channels to obtain wave impedance inversion results for all channels.

Description

Post-stack wave impedance inversion method based on Cauchy distribution

The invention relates to geophysical exploration technology, belonging to the reservoir prediction inversion technology class, and is a post-stack wave impedance inversion method based on Cauchy distribution.

Background technique

Seismic exploration is to artificially excite seismic waves, record the seismic response of seismic waves on the surface or underground with single-component or multi-component sensors, study their propagation patterns in the formation, and identify the underground geological structures through seismic data processing and inversion methods. Lithological characteristics, and then search for geophysical exploration methods for mineral resources such as oil and gas. Seismic exploration begins with understanding the subsurface structural form, developing direct application of seismic information to determine lithology, analyzing lithofacies, quantifying the physical parameters of rock formations, and finding oil and gas displays. Earthquake inversion technology is the product of this development process.

The basic purpose of seismic inversion is to use the propagation law of seismic waves in underground media, and to estimate the spatial distribution of underground rock stratum structure and physical parameters through data collection, processing and interpretation processes, which provides an important basis for exploration and development. In seismic inversion studies, there are various parameters of inversion methods such as wave impedance, velocity, density, porosity, permeability, Poisson's ratio, and so on. Since the wave impedance information is a bridge connecting geology and geophysics, the amount of calculated data after stacking is relatively small, and it is convenient and effective in actual production. Therefore, wave impedance inversion has a special position in seismic inversion, earthquake Inversion usually refers to wave impedance inversion.

Conventional seismic wave impedance inversion refers to seismic special processing technology that uses seismic materials to invert formation/rock wave impedance. Compared with statistical multi-parameter pattern recognition prediction reservoir oil and gas, neural network prediction formation parameters, amplitude fitting prediction reservoir thickness and other statistical methods, wave impedance inversion has a clear physical meaning, it is reservoir lithology prediction, oil The deterministic method of reservoir feature description has achieved remarkable results in practical applications. Geological effects. Most of the current inversion methods are model-based methods. These methods generally establish an initial model based on logging and geological data, and iteratively obtains lithology parameters by generalized linear inversion. Due to the nonlinearity of the problem, in addition to requiring fine wavelets, the initial model is required to be close to the real model in order to achieve reliable results, that is, the inversion results strongly depend on the choice of the initial model. In addition to such methods, globally optimized inversion methods (such as genetic algorithms and simulated annealing algorithms) overcome the drawbacks of the model-based approach to the initial model, but because of the globally optimal inversion results, The inversion is very slow. Summary of the invention

The object of the present invention is to provide a post-stack impedance inversion method based on Cauchy distribution.

The invention is achieved by the following technical solutions:

1) Collecting seismic data by conventional seismic exploration methods, and performing conventional processing on seismic data to obtain post-stack seismic data;

2) performing post-stack seismic data on the layered data to obtain layer data, verifying and correcting the determined target layer level, and interpolating and smoothing;

3) Using conventional logging methods to obtain logging data, and obtaining time-lapse curves and density curves of logging acoustic waves;

4) According to the post-stack seismic data, the horizon data and the known drilling stratification data, the acoustic time difference curve and the density curve of the depth domain are calibrated into a time domain curve, and the time domain wave impedance curve data in the well is generated, and Extracting seismic wavelets from the well-side seismic trace;

The calibration is to simulate the well-side seismic record by using the logging curve and the seismic wavelet, and realize the calibration and mapping of the logging layer to the seismic horizon, thereby obtaining the time-depth relationship curve, and thus the depth-depth relationship can be the depth domain. The log curve is converted to a time domain curve.

5) Using the horizon data of step 2) and the time domain impedance curve obtained in step 4) Initial wave impedance model;

6) Normalize the seismic data and the wavelet data obtained in step 4) and normalize them to the range [-1, 1];

7) A seismic data read in step 1), a data of the initial wave impedance model generated in step 5) and the wavelet data extracted in step 4) are substituted into the following formula, and the inversion is obtained by inversion. Sequence of reflection coefficients: r = (G ^T G + +

(G ^T d + pC ^T ξ) (1) where is the post-stack seismic data, Ν is the total number of sampling points of the seismic data; r = [^r ₂ ,...,r _N ] ^T is the sequence of reflection coefficients; G is N XN dimensional wavelet convolution matrix, the upper angle T represents the transposition of the matrix; μ is the sparse constraint factor, which controls the sparsity of the reflection coefficient; the diagonal element of the matrix Q is: „ represents the η row of the matrix Q The value of the column element, η is the row and column number of the Q matrix,

The matrix Q except the diagonal elements are all zero, ^ represents the standard deviation of the noise, R„ is the initial reflection coefficient at the position of the nth sample point calculated by the initial wave impedance model; The degree of dependence of the results on the initial model; C is the integral operator matrix, and its discrete form is expressed as:

1 0 ··· 0

1 1 0 :

C =

: '·. '·. 0

1 1 1 Equation (1) The upper-angle-1 is the inverse of the matrix; 丄l _n ^ =∑. is the nth sampling point

2 two 1

The impedance value of the wave, /. In order to invert the initial wave impedance value corresponding to the first sampling point in the time window, /„ is the initial wave impedance value at the nth sampling point in the inversion, In is the natural logarithmic symbol, which is the inversion time window. The reflection coefficient value of the 1st water sampling point indicates the request for the ^ from the 1st sampling point to the nth sampling point. 8) by the definition of the relative wave impedance 丄ln = ^. at the nth sampling point in step 7),

twenty one

The relationship between the wave impedance and the reflection coefficient can be obtained by derivation:

Where: /„ is the wave impedance value of the nth sampling point in the inversion window, / is the initial wave impedance value corresponding to the first sampling point in the inversion window, which is obtained by inversion in step 7) The reflection coefficient value of the first sampling point, e represents the bottom of the natural logarithm, and represents the summation calculation for the ^ from the first sampling point to the nth sampling point.

By converting the series of reflection coefficients obtained by the inversion of step 7) into a wave impedance sequence by formula (2), the wave impedance inversion result of the channel can be obtained;

9) Repeat steps 7) through 8) for all seismic traces to obtain the wave impedance inversion results for all channels.

The invention has the following characteristics, and the main performances are as follows:

(1) By assuming that the prior probability of the reflection coefficient obeys the Cauchy distribution, the Cauchy distribution belongs to the long tail distribution, and the distribution at the peak is narrower than the Gaussian distribution, and the velocity approaching zero is also slow, resulting in a small amount. Non-zero values and a large number of zero values to achieve sparse pulse inversion;

(2) The accuracy and stability of the inversion results can be controlled by adding wave impedance model constraints;

(3) Sparse pulse inversion is a seismic-based inversion method. The resolution, signal-to-noise ratio and reliability of the inversion results mainly depend on the quality of the seismic data itself. The seismic noise is sensitive to the inversion results and affects Large, so seismic data for sparse pulse inversion should have features such as wider frequency bands, lower noise, relative amplitude retention, and accurate imaging. Logging data, especially sonic logging and density logging data, are the comparison criteria and interpretation basis for seismic lateral prediction. They should be carefully edited and corrected before the inversion process to correctly reflect the physical characteristics of the rock formation; (4) The wave impedance result obtained by the sparse pulse inversion has a high degree of coincidence with the well curve at the well position, and the calculation efficiency is high.

The inversion method of the present invention not only has the characteristics of the general recursive inversion method, that is, the inversion result is faithful to the seismic data and can reflect the lateral variation of the reservoir. Moreover, the introduction of geological and logging data into the inversion constraints during the iterative process increases some of the low frequency and high frequency components and broadens the inversion frequency band to some extent. This method is less dependent on the initial model, and the inversion results are more unique and less prone to artifacts.

DRAWINGS

Figure 1 Application example of the present invention 1 Post-stack seismic profile of model inversion results;

2 is a wave impedance profile of the model inversion result of the application example 1 of the present invention;

Figure 3 is a comparison of the inversion wave impedance results of the present invention in the Dagang Port area with the Jason software inversion results. detailed description

The method provided by the present invention is a post-stack wave impedance inversion method based on Cauchy distribution, and the operation efficiency is very high.

The invention is achieved by the following technical solutions:

The calibration is to simulate well logging by using well logs and seismic wavelets to achieve logging stratification. To the calibration and mapping of the seismic horizon, the time-depth relationship curve is obtained, and the time-depth relationship can convert the logging curve in the depth domain into the time domain curve.

5) using the horizon data of step 2) and the time domain impedance curve obtained in step 4) to generate an initial wave impedance model;

7) A seismic data read in step 1), a data of the initial wave impedance model generated in step 5), and a wavelet data extracted in step 4) are substituted into the following formula, and the inversion is obtained by inversion. Sequence of reflection coefficients: r = (G ^T G + +

(G ^T d + pC ^T ξ) ( 1 ) where is the post-stack seismic data, Ν is the total number of sampling points of the seismic data; r = [^ r ₂ , ..., r _N ] ^T is the sequence of reflection coefficients; N XN dimensional wavelet convolution matrix, the upper angle T represents the transposition of the matrix; μ is the sparse constraint factor, which controls the sparsity of the reflection coefficient; the diagonal element of the matrix Q is: „ represents the η row of the matrix Q The value of the column element, η is the row and column number of the Q matrix,

In the formula (1), the upper mark-1 is the inverse of the matrix; =丄1 ₁₁ = ^ is the 11th sample point

2 2 1 pair of wave impedance values, /. In order to invert the initial wave impedance value corresponding to the first sampling point in the time window, /„ is the initial wave impedance value at the nth sampling point in the inversion, In is the natural logarithmic symbol, which is the inversion time window. 1st The reflection coefficient value of each sampling point represents the summation calculation of ^ from the first sampling point to the nth sampling point.

8) by the definition of the relative wave impedance 丄 ln = ^. at the nth sampling point in step 7),

twenty one

2Y r,

I _n = e (2) where: / „ is the wave impedance value of the nth sampling point in the inversion window, / is the initial wave impedance value corresponding to the first sampling point in the inversion window, which is the step 7) The first sample point reflection coefficient value obtained by the inversion, e represents the bottom of the natural logarithm, and represents the summation calculation for the ^ from the 1st sample point to the nth sample point.

The effects of the present invention will be described below by way of specific examples.

Invention Application Example 1 :

Using this patented method, a two-dimensional model data (from Jason Corporation) was tested. The vertical direction of the model is three diamond sand bodies. The fluids contained in each sand body from top to bottom are gas, oil and water, respectively. Figure 1 shows the noise-free post-stack seismic profile of the model. Based on the previous steps 1) -9), the sparse pulse inversion based on the Cauchy distribution is obtained, and the wave impedance profile of Fig. 2 is obtained. It can be seen from Fig. 2 that the inversion result is in good agreement with the well curve at the well position, and the impedance difference inside the sand body can be clearly reflected from the figure, indicating that the inversion result is very accurate. The inversion time window is 1250 milliseconds, the sampling interval is 4 milliseconds, a total of 200 channels, and the inversion on the HP2 single machine by the present invention takes about 110 seconds.

Invention Application Example 2:

Using this patented method to invert the Dagangkoukou area, according to the previous steps 1) -9) to perform the sparse pulse inversion based on the Cauchy distribution, Figure 3 (a) is the wave impedance profile obtained by the present invention, Figure 3 (b) For the wave impedance profile obtained by inversion using the foreign software Jason system, it can be seen from the figure that the inversion result of the present invention is generally comparable to the Jason system inversion result.

The inversion target layer time window is about 900 milliseconds, the sampling interval is 4 milliseconds, and a total of 300 channels. It takes about 340 seconds to perform the inversion on the HP2 stand-alone unit with the present invention.

The theoretical model and the actual data show that the Cauchy distribution can not only restore the sparseness of the reflection coefficient sequence, but also achieve a balance between improving the resolution of seismic data and reducing the suppression of weak reflection information. The inversion results obtained by this method are not only accurate, but also highly accurate, and the calculation efficiency is very high. This result can be used for accurate reservoir prediction.

Claims

claims

1. A post-stack wave impedance inversion method based on Cauchy distribution. The steps are as follows:

1) Use conventional seismic exploration methods to collect seismic data, and perform conventional processing on the seismic data to obtain post-stack seismic data;

2) Perform layer picking on post-stack seismic data to obtain layer data, and perform inspection, correction, interpolation and smoothing on the determined target layer layers;

3) Use conventional logging methods to obtain logging data, and obtain the logging acoustic transit time curve and density curve;

4) Based on the post-stack seismic data, layer data and known drilling layering data, calibrate the acoustic travel time difference curve and density curve in the depth domain into curves in the time domain, and at the same time generate the time domain wave impedance curve data in the well, and Extract seismic wavelets from seismic traces next to the well;

The calibration described in step 4) is to use well logging curves and seismic wavelets to simulate wellside seismic records to achieve calibration and mapping of well log layers to seismic layers, thereby obtaining the time-depth relationship curve, and from this time-depth relationship Depth domain logging curves can be converted into time domain curves.

5) Use the layer data in step 2) and the time domain wave impedance curve obtained in step 4) to generate an initial wave impedance model;

6) Normalize the seismic data and the wavelet data obtained in step 4) respectively, and normalize them to the range [-1, 1];

7) Substitute the seismic data read in step 1), the initial wave impedance model data generated in step 5) and the wavelet data extracted in step 4) into the following formula, and obtain the trace of the trace through inversion. Reflection coefficient sequence:

r = (G ^T G + +

(G ^T d + pC ^T ξ) (1) where is the post-stack seismic data, Ν is the total number of sampling points of the seismic data; r = [^ r ₂ , ..., r _N ] ^T is the reflection coefficient sequence; G is N XN dimensional wavelet convolution matrix, the superscript T represents the moment The transpose of the matrix; μ is the sparsity constraint factor, which controls the sparsity of the reflection coefficient; the diagonal elements of the matrix Q are:

„Represents the value of the nth row and nth column element of the matrix Q, n is the row and column number of the Q matrix,

Except for the diagonal elements of the matrix Q, all other elements are zero, ^ represents the standard deviation of the noise, R„ is the initial reflection coefficient at the nth sampling point calculated by the initial wave impedance model; is the model constraint factor, controlling the inverse The degree of dependence of the simulation results on the initial model; C is the integral operator matrix, and its discrete form is expressed as:

In formula (1), the superscript -1 is the inversion of the matrix; =丄1 ₁₁ = ^ is the value at the 11th sampling point

twenty one

Wave impedance value, /. is the initial wave impedance value corresponding to the first sampling point within the inversion time window, /„ is the initial wave impedance value at the nth sampling point within the inversion time window, In is the natural logarithm sign, and is the initial wave impedance value within the inversion time window Reflection coefficient value of the 1st water sampling point, ^ represents the calculation of ^ from the 1st sampling point to the nth sampling point

8) Through the definition of relative wave impedance 丄l _n = ^. at the nth sampling point in step 7)

twenty one

The relationship between wave impedance and reflection coefficient is obtained by derivation:

2∑r,

(2) Among them: /„ is the wave impedance value of the nth sampling point in the inversion time window, /. is the initial wave impedance value corresponding to the first sampling point in the inversion time window, is the inversion value in step 7) The reflection coefficient value of the first sampling point obtained by deducing, e represents the base of the natural logarithm, and represents the summation calculation of ^ from the first sampling point to the nth sampling point; Use formula (2) to convert the reflection coefficient series of a channel obtained by the inversion in step 7) into a wave impedance series, and then the wave impedance inversion result of the channel can be obtained;

9) Repeat steps 7) to 8) for all seismic channels to obtain the wave impedance inversion results of all channels.