Inverse filtering processing for improving longitudinal resolution

According to the propagation theory of seismic waves, seismic waves propagate underground in the form of seismic wavelets in viscoelastic media, and the seismic records of reflected waves received on the ground are the convolution of the reflection coefficient of strata and seismic wavelets. Therefore, the strata are equivalent to a filter, which makes the reflection coefficient sequence become a seismic record composed of wavelets, reducing the vertical resolution of seismic exploration. The purpose of inverse filtering is to design an inverse filter, and then filter the seismic records to eliminate the effect of stratum filtering and improve the longitudinal resolution of seismic records.

3.3.3.1 seismic wavelet

seismic wavelet is a signal with definite starting time and finite energy and finite duration, and it is the basic unit of seismic wave in seismic records. It is generally believed that the seismic wave generated when the earthquake source is excited is only a sharp pulse with a very short duration. As the sharp pulse propagates in viscoelastic media, the high-frequency component of the sharp pulse quickly attenuates, and the waveform increases accordingly, becoming a seismic wavelet with a limited frequency band width and a certain duration. Seismic waves propagate underground in the form of seismic wavelets.

3.3.3.1.1 Mathematical model of seismic wavelet

In practice, seismic wavelet is a very complicated problem, because seismic wavelet is related to formation lithology, which itself is a complex body. However, for the convenience of research, it is still necessary to simulate the seismic wavelet. At present, it is generally believed that the mathematical model of seismic wavelet proposed by Lake is widely representative, which is called Lake wavelet. The mathematical model of the minimum phase seismic wavelet is

b(t)=e-αt2sin2π? T (3.3-25)

where: Is the main frequency of wavelet; α=2? 2ln(M) is the wavelet attenuation coefficient; M = | m1/m2 | is the ratio of the maximum wave peak value m1 to the maximum wave trough m2. The wavelet shape calculated by formula (3.3-25) is shown in Figure 3-16.

figure 3-16 schematic diagram of lake wavelet

3.3.3.2 phase problem of wavelet

if the spectrum of seismic wavelet b(t) is calculated by Fourier transform as B(ω), then

b (ω) = | b (ω) | EJ φ (ω) (3.3-26. Like any wave function, the characteristics of this wave function can be described by its amplitude spectrum and phase spectrum. For complex and changeable wavelets, the most frequent changes are the attenuation form and duration of the waveform. Therefore, generally, seismic wavelets have a relatively stable amplitude spectrum, but a relatively large phase spectrum. If

ψ(ω)=-φ(ω) (3.3-27)

, ψ (ω) is called the phase delay spectrum. The wavelets with the same amplitude spectrum can be classified according to their different phase delay spectra.

the phase delay of seismic wavelet is also a relative concept, for example, there are two terms (only two values) that wavelet B1 (n) = {1,.5} and B2 (n) = {.5,1}, which can be written as general formula B (n) = {A,a1}. Get z transform

B(z)=a+a1z (3.3-28)

Then use z = e-j Ω δ to get amplitude spectrum and phase spectrum

Principle, method and explanation of seismic exploration

Substitute a=1, a1=.5 and a=.5, a1=1. It can be seen that the two sub-waves b1(n) and b2(n) have the same amplitude spectrum, but the phase delay spectrum is different. Because Φ 1 (ω) has less phase delay than Φ 2 (ω), it is called the minimum delay phase. However, the phase delay of φ2(ω) is larger than that of φ1(ω), which is called the maximum delay phase.

The minimum and maximum phasor wavelets can be described as finding the zero α of the wavelet z transformation equation (a+a1z)= in the z domain (Z plane), that is, there are

principles, methods and explanations of seismic exploration

When | α | > 1, wavelet b(n) is called the minimum phasor wavelets.

| α | < 1, wavelet b(n) is called the largest phase wavelet.

| α | = 1, wavelet b(n) is called mixed phase wavelet.

the phase classification of binomial wavelets can be extended to the classification of n wavelets, that is, let b(n)={b, b1, …, bn}, and its z transformation is

B(z) = B+B1z+B2Z2+…+BNZn (3.3-33)

to find b (z). If all zeros are in the unit circle (| α I | < 1, I = 1,2, ..., n), then b(n) is the largest phase wavelet. If zeros exist both inside and outside the unit circle, then b(n) is a mixed phase wavelet. It can be seen that the phase of the wavelet can be completely judged according to the position of the zero point of the wavelet Z-transform polynomial in the Z-plane. Therefore, the Z-plane can be divided into two areas by the unit circle, with the minimum phase area outside the unit circle and the maximum phase area inside the unit circle, as shown in Figure 3-17. The waveform and energy characteristics of the other three kinds of phase wavelet are: the energy of the smallest phase wavelet is mainly concentrated in the front, the energy of the largest phase wavelet is mainly concentrated in the back, and the energy of the mixed phase wavelet is mainly concentrated in the middle. Three kinds of phase wavelet are shown in Figure 3-18. In practice, seismic wavelets are mainly minimum phase wavelet and mixed phase wavelet.

Phase delay property of zero point position indication wavelet on Z plane in Figure 3-17

3.3.3 Principle and method of inverse filtering

As mentioned above, seismic record is the convolution of strata reflection coefficient sequence r(t) and seismic wavelet b(t), that is,

x(t)=r(t)? B(t)(3.3-34)

Due to the problem of wavelet, the high-resolution reflection coefficient pulse sequence becomes a low-resolution seismic record, and b(t) is equivalent to the formation filter factor. In order to improve the resolution, an inverse filter can be designed, and the inverse filter factor is a(t), and it is required that a(t) and b(t) satisfy the following relationship

a(t)? B (t) = Δ (t) (3.3-35)

Figure 3-18 m+1 signal classification and characteristics of counter-signal

Using a(t) to filter seismic record x(t)

x(t)? a(t)=r(t)? b(t)? a(t)=r(t)? δ(t)=r(t)(3.3-36)

The result is a reflection coefficient sequence. The above is the basic principle of inverse filtering.

the core of inverse filtering is to determine the inverse filtering factor a(t). due to the uncertainty of seismic wavelet and the existence of noise interference in seismic records, it is very difficult or even impossible to determine the accurate a(t) in practice. Therefore, under different approximate assumptions, many methods for determining the inverse filter factor a(t) have been studied, which can be basically divided into two categories: one is to obtain the seismic wavelet b(t) first, and then obtain A (t) according to b(t); Another category is to find a(t) directly from seismic records. There are many different methods in each category (only the number of inverse filtering methods shows the difficulty of inverse filtering processing). In the following, several representative inverse filtering methods are discussed.

3.3.3.1 formation inverse filtering

formation inverse filtering belongs to the method of finding wavelet b(t) first, and then finding a(t). This method requires logging data and good seismic records beside wells. Firstly, the formation reflection coefficient sequence r(t) matched with the seismic trace x(t) near the well is converted from acoustic logging data, and the frequency domain equation of r(t) and x(t) is

X (ω) = R (ω) B(ω) (3.3-37)

, that is, there is

. Where A(ω) is the frequency spectrum of the inverse filter factor. Written as Z transformation, A(z)=, it can be seen that A(z) is a rational fraction. To make A (z) stable, the root of denominator polynomial B(z) must be outside the unit circle, that is, wavelet b(t) is required to be the minimum phase.

The wavelet and inverse filtering factor can be obtained by using logging and seismic traces near the well, and then the inverse filtering factor can be used to perform inverse filtering on other traces of the survey line.

3.3.3.2 Least square inverse filtering

Least square inverse filtering is the application of least square filtering (or Wiener filtering and optimal filtering) in the field of inverse filtering.

the basic idea of least square inverse filtering is to design a filter operator, which can be used to convert a known input signal into an output that is best close to a given expected output signal in the sense of minimum square error.

let the input signal be x(t), which is convolved with the filter factor h(t) to get the actual output y(t), that is, y(t)=x(t)? h(t)。 For various reasons, the actual output y(t) can't be exactly the same as the expected output (t) given in advance, and only the two can be required to be optimally close. There are many criteria to judge whether it is the best approach, and the least square error criterion is one of them, that is, when the sum of the squares of their errors is the smallest, it means that they are the best approach. In this sense, the filtering by finding the filter factor h(t) is the least square filtering.

if the filter factor to be obtained is the inverse filter factor a(t), the expected output after inverse filtering of the input wavelet b(t) is d(t), and the actual output is y(t), according to the least square principle, the inverse filter factor obtained by minimizing the sum of the squares of their errors is called the least square inverse filter factor, and the inverse filtering of the seismic record x(t) with it is the least square inverse filter.

3.3.3.2.1 Basic equation of least square inverse filtering

Let the input discrete signal be seismic wavelet b(n)={b(), b(1), ..., b(m)}, and the inverse filtering factor to be solved is a(n)={a(m), a. The convolution of b(n) and a(n) is the actual output y(n), that is,

principles, methods and explanations of seismic exploration

principles, methods and explanations of seismic exploration

In order to minimize Q, it is mathematically a problem of seeking the extreme value of Q, that is, seeking to satisfy the principles and methods of seismic exploration.

principles, methods and explanations of seismic exploration

principles, methods and explanations of seismic exploration

are autocorrelation functions of seismic wavelet, while

principles, methods and explanations of seismic exploration

are cross-correlation functions of seismic wavelet and expected output, so equation (3.3-41) can be written as

principles, methods and explanations of seismic exploration

. In this equation, the coefficient matrix is a special positive definite matrix, which is called the general Tobritz matrix. This matrix equation can be quickly solved by levinson recursive algorithm.

equation (3.3-45) is the basic equation of least square inverse filtering. The equation adapts wavelet b(n) to minimum phase, maximum phase and mixed phase. In the formula, the starting time m of the inverse filter factor a(n) is related to the phase of the wavelet, and its value rule is determined by the Z transformation of the wavelet and the inverse filter factor.

because m+1 seismic wavelets b(n)={b(), b(1), ..., b(m)} are physically realizable signals, whose z-transform is B(z), the inverse signal of wavelet b(n) is a(n), and the z-transform of a(n) is A(z).

1) when b(n) is the minimum phase, all the roots αi satisfy | α I | > 1, and each term in the formula (3.3-46) is expanded into a power series of z, and the lowest power of z in each factor is , and the lowest power of z is still after the multiplication of m factors, so this formula can be written as (taking m+1 term)

a.

2) when b(n) is the maximum phase, all αi satisfy | α I | < 1, and the power series of expanding A(z) into z is

a (z) = …+a-m-3z-m-3+a-m-2z-m-2+a.

3) when b(n) is mixed phase, the root | α I | exists both inside and outside the unit circle, and | α I | ≠ 1. If the roots in the unit circle are treated as the maximum phase and the roots outside the unit circle are treated as the minimum phase, the inverse filtering factor of the mixed-phase wavelet has values on both sides of the time coordinate, and the number of values on both sides depends on the distribution number of roots | α i | inside and outside the unit circle.

It can be seen that the inverse filtering factors of seismic wavelets are quite different with different phases. The corresponding relationship between different phase wavelet and inverse filtering factor is shown in Figure 3-18.

in addition, when solving equation (3.4-16), it is necessary to select the function form that you want to output d(t). Generally, d(t) can be selected as

Principles, methods and explanations of seismic exploration

When d(t)=δ(t) in the equation, the output is a pulse, and a(n) is called a pulse. If the second item in the formula is selected, the output is a zero-phase waveform with the main frequency of? , which can play the role of wavelet shaping or phase conversion, then a(n) is called the inverse filter factor of wavelet shaping. The above processes are all based on the assumption that the wavelet is known, so they are collectively called wavelet processing or wavelet inverse filtering. In practical application, if Equation (3.3-44) is written as

Principles, methods and explanations of seismic exploration

, the inverse filter factor to be solved is bilateral inverse filter factor.

3.3.3.2.2 Least Square Inverse Filtering with Unknown Wavelet

Generally, seismic wavelet is unknown.

in order to find the inverse filter factor without knowing the wavelet, it is necessary to add certain restrictions or assumptions to the seismic wavelet and reflection coefficient sequence, which include

1) assuming that the reflection coefficient sequence r(t) is a random white noise sequence, that is, its autocorrelation is

Principles, methods and explanations of seismic exploration

2) assuming that the seismic wavelet is the smallest phase.

according to the first hypothesis, the autocorrelation rbb(τ) of seismic wavelet can be replaced by the autocorrelation of recorded x(t), because

seismic exploration.