3.3.3. 1 seismic wavelet
Seismic wavelet is a signal with definite starting time, limited energy and limited duration, and it is the basic unit of seismic wave in seismic records. It is generally believed that the seismic wave generated when the source is excited is only a sharp pulse with a very short duration. When the sharp pulse propagates in viscoelastic medium, the high frequency component of the sharp pulse decays rapidly, and the waveform increases accordingly, becoming a seismic wavelet with limited frequency band width and certain duration. Seismic waves propagate underground in the form of seismic wavelets.
3.3.3. 1. 1 seismic wavelet mathematical model
In practice, seismic wavelet is a very complicated problem, because seismic wavelet is related to formation lithology, which itself is a complex. 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 and called Lake wavelet. The mathematical model of 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 = | m 1/m2 | is the ratio of the maximum peak m 1 to the maximum trough m2. The wavelet shape calculated by formula (3.3-25) is shown in Figure 3- 16.
Fig. 3- 16 schematic diagram of lake wavelet
Phase problem of 3.3.3.2 wavelet
If the spectrum of seismic wavelet b(t) is found to be B(ω) by Fourier transform, there are
B(ω)=| B(ω)| ejφ(ω) (3.3-26)
Where | b (ω) | is the amplitude spectrum of wavelet, and φ (ω) is the phase spectrum of wavelet. 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 speaking, seismic wavelet has a relatively stable amplitude spectrum, but a relatively large phase spectrum. If you take it with you
ψ(ω)=-φ(ω) (3.3-27)
ψ (ω) is called phase delay spectrum. Wavelet 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): wavelet B 1 (n) = {1, 0.5}, and wavelet B2 (n) = {0.5, 1}, which can be written as the general formula b(n)={a0, a/kloc-0.
B(z)=a0+a 1z (3.3-28)
The amplitude spectrum and phase spectrum can be obtained by z = e-j ω δ.
Principles, methods and explanations of seismic exploration
Substituting a0= 1, a 1=0.5 and a0=0.5, a 1= 1, the result is
Principles, methods and explanations of seismic exploration
And phase delay spectrum.
Principles, methods and explanations of seismic exploration
Where δ is the sampling rate. It can be seen that the amplitude spectra of the two wavelets b 1(n) and b2(n) are the same, but the phase delay spectra are different. Because the phase delay of φ1(Ω) is less than φ 2 (Ω), it is called minimum delay phase, which is called minimum phase for short. The phase delay of φ2(ω) is greater than φ 1(ω), which is called the maximum delay phase.
In the Z-domain (Z-plane), the minimum and maximum phase wavelets can be described to find the zero α of the wavelet Z-transform equation (a0+a 1z)=0, that is, there are
Principles, methods and explanations of seismic exploration
When | α| > 1, wavelet b(n) is called minimum phase wavelet.
When | α| < 1, wavelet b(n) is called maximum phase wavelet.
When | α| = 1, wavelet b(n) is called mixed phase wavelet.
The phase classification of binomial wavelet can be extended to N-term wavelet classification, that is, let b(n)={b0, b 1, …, bn}, and its z transform is
b(z)= B0+b 1z+b2z 2+…+bnzn(3.3-33)
Find all zeros (i.e. roots) of B(z) polynomial αi, i= 1, 2, …, n. If all zeros are outside the unit circle (| α i | > 1, i= 1, 2, …, n), then b(n) If all zeros are in the unit circle (| α I | < 1, i= 1, 2, …, n), then b(n) is the maximum phase wavelet. If there are zeros inside and outside the unit circle, then b(n) is a mixed-phase wavelet. It can be seen that the phase of wavelet can be judged completely according to the position of the zero point of wavelet Z transform polynomial on the Z plane. So the Z plane can be divided into two areas, with the unit circle as the external minimum phase area and the internal maximum phase area, 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.
Fig. 3-Phase delay characteristics of zero position indicating wavelet in Z plane 17
Principle and method of 3.3.3.3 inverse filtering
From the above, it can be seen that the seismic record is the convolution of the strata reflection coefficient sequence r(t) and the 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. The inverse filter factor is a(t), and a(t) and b(t) are required to satisfy the following relationship.
a(t)？ b(t)=δ(t)(3.3-35)
Figure 3- Signal Classification and Counter-signal Characteristics of Project 18m+ 1
Filtering seismic record x(t) with a(t)
x(t)？ a(t)=r(t)？ b(t)？ a(t)=r(t)？ δ(t)=r(t)(3.3-36)
The result is a sequence of reflection coefficients. The above is the basic principle of inverse filtering.
When realizing inverse filtering, the core is to determine the inverse filtering factor a(t). Because of the uncertainty of seismic wavelet and the existence of noise interference in seismic records, it is difficult or even impossible to determine the accurate a(t) in practice. Therefore, under different approximate assumptions, many methods to determine the inverse filter factor a(t) have been studied, which can be basically divided into two categories: one is to first find the seismic wavelet b(t), and then find A (t) according to b(t); The other is to find a(t) directly from earthquake records. There are many different methods for each category (the difficulty of inverse filtering can be seen only from the number of inverse filtering methods). Several typical inverse filtering methods are discussed below.
3.3.3.3.1formation inverse filtering
Formation inverse filtering belongs to the method of finding wavelet b(t) first and then finding a(t). This method needs logging data and good seismic records beside wells. Firstly, the formation reflection coefficient sequence r(t) matched with seismic trace x(t) near the well is converted from acoustic logging data, and the frequency domain equations of r(t) and x(t) are obtained as follows.
X(ω)=R(ω) B(ω) (3.3-37)
Have it at once
Principles, methods and explanations of seismic exploration
Where B(ω) is the frequency spectrum of wavelet b(t), the relationship between wavelet and inverse filter factor is
Principles, methods and explanations of seismic exploration
A(t) is obtained by inverse Fourier transform. Where A(ω) is the frequency spectrum of the inverse filter coefficients. Written as z transformation, A(z)=, so 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.
Wavelet and inverse filtering factor can be obtained by using logging and seismic traces near the well, and the inverse filtering factor can be used to filter other traces of the survey line.
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 to convert the known input signal into the output closest to the given expected output signal in the sense of minimum square error.
Let the input signal be x(t) and convolve with the filter factor h(t) to get the actual output y(t), that is, y(t)=x(t)? Hamilton, due to various reasons, the actual output y(t) can't be exactly the same as the expected output (t) given in advance, and only the optimal approximation of the two can be required. There are many criteria to judge whether they are the best approximation, and the least square error criterion is one of them, that is, when the sum of their errors is the smallest, it means that they are the best approximation. 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 their error squares is called the least square inverse filter factor, and it is the least square inverse filter to filter the seismic record x(t).
3.3.3.3.2. 1 Least Square Inverse Filtering Basic Equation
Let the input discrete signals be seismic wavelets b(n)={b(0), b( 1), ..., b(m)}, and the inverse filter factor of the solution is required to be a(n)={a(m0), a(m0+ 1), a (m0.
Principles, methods and explanations of seismic exploration
The sum of the squares of the errors between the actual output and the expected output is
Principles, methods and explanations of seismic exploration
Minimizing q is mathematically a problem of finding the extreme value of q, that is, finding satisfaction.
Principles, methods and explanations of seismic exploration
The filter factor a(t).
Principles, methods and explanations of seismic exploration
Principles, methods and explanations of seismic exploration
Is the autocorrelation function of the seismic wavelet, and
Principles, methods and explanations of seismic exploration
Is the cross-correlation function between seismic wavelet and expected output, so equation (3.3-4 1) can be written as follows.
Principles, methods and explanations of seismic exploration
This is a system of linear equations, written in matrix form as follows
Principles, methods and explanations of seismic exploration
The symmetry of autocorrelation function is used in the formula. In this equation, the coefficient matrix is a special positive definite matrix called generalized Tobriz matrix. This matrix equation can be solved quickly by levinson recursive algorithm.
Equation (3.3-45) is the basic equation of least square inverse filtering. The equation applies wavelet b(n) to minimum phase, maximum phase and mixed phase. In the formula, the starting time m0 of the inverse filter factor a(n) is related to the phase of the wavelet, and its value rule is determined by the Z transform of the wavelet and the inverse filter factor.
Because m+ 1 seismic wavelet b(n)={b(0), b( 1), ..., b(m)} are physically realizable signals, their z-transform is B(z), and the inverse signals of wavelet b(n) are a(n), a (.
Principles, methods and explanations of seismic exploration
If the inverse Z transform of A(z) is found, it is the inverse filter factor a(n), and the position of a(n) is determined by the position of the root ai of wavelet Z transform B(z) in the Z plane.
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, the lowest power of z in each factor is 0, and the lowest power of z is still 0 after the multiplication of m factors, so this formula can be written as (take
a(z)= a0z 0+a 1z 1+…+anzn+…+amzm(3.3-47)
This shows that the starting time of the inverse filter factor a(n) is m0=0, and when n < m0, a(n)=0.
2) when b(n) is the maximum phase, all αi satisfy | α i | < 1, and the power series of A(z) expanded into z is
a(z)=…+a-m-3z-m-3+a-m-2z-m-2+a-m- 1z-m- 1+a-mz-m
If m+ 1 is also taken, the corresponding inverse filtering factors a(n)={a(-m-m), a(-m+1), ..., a(-m-3), a(-m-2), a (.
3) When b(n) is a mixed phase, the root | α i | exists both inside and outside the unit circle, and | α i | ≠ 1. If the root in the unit circle is taken as the maximum phase and the root outside the unit circle as the minimum phase, the inverse filter 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 inside and outside the unit circle | α i |.
It can be seen that the inverse filtering factors of seismic wavelets in different phases are quite different. The corresponding relationship between different phase wavelets and inverse filtering factors 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 speaking, you can choose d(t) as.
Principles, methods and explanations of seismic exploration
When d(t)=δ(t) in the formula, the output is a pulse, and a(n) is called the pulse inverse filter factor. If the second term in the formula is selected, the output is the main frequency? 0, 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 wavelet is known, so they are collectively called wavelet processing or wavelet inverse filtering. If Formula (3.3-44) is written in practical application,
Principles, methods and explanations of seismic exploration
The inverse filter factor solved is a bilateral inverse filter factor.
3.3.3.3.2.2 Least Square Inverse Filtering of Unknown Wavelet
Usually, the seismic wavelet is unknown.
In order to find the inverse filter factor without knowing the wavelet, it is necessary to add some restrictions or assumptions to the seismic wavelet and reflection coefficient sequence, including
1) suppose 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) Assume that the seismic wavelet is the minimum phase.
According to the first hypothesis, the autocorrelation rbb(τ) of seismic wavelet can be replaced by the autocorrelation of record x(t), because
Principles, methods and explanations of seismic exploration
Principles, methods and explanations of seismic exploration
According to the second hypothesis, we can know that the zeros of Z-transform B(z) of seismic wavelet b(t) are all outside the unit circle, that is, the poles of Z-transform A(z)= 1/B(z) of inverse wave factor a(t) are all outside the unit circle, so a(t) is stable and physically realizable. Therefore, when t < 0, m0=0, and the free term becomes b(0), b(- 1), …, b (-m) T. Because b (-m) t) must be physically realizable, b(- 1)=0, b (-. Let a'(t)=a(t)/b(0), then the basic equation becomes
Principles, methods and explanations of seismic exploration
This is the basic equation for obtaining the inverse filter factor when the wavelet is unknown, and the elements in the coefficient matrix can be obtained directly from the seismic records. The difference between the obtained inverse filtering factor A ′ (t) and a(t) is only several times, which does not affect the inverse filtering effect of compressing wavelet and improving resolution. It is also commonly called sharp pulse inverse filtering.
3.3.3.3.2.3 Statistical method for wave calculation
Although wavelets are generally unknown, seismic records contain wavelets, so wavelets can be obtained from seismic records. At present, there are many methods to find seismic wavelet. The following is the method of finding wavelet by multi-channel statistics.
If wavelet is regarded as a general signal, it can also be expressed by s(t). It has been proved in the previous section that if the reflection coefficient is a random white noise sequence and the autocorrelation x(t) of the seismic record and the autocorrelation s(t) of the wavelet are equal, then the recorded amplitude spectrum | x (ω) | and the wavelet amplitude spectrum | s (ω) | are equal.
| S(ω)|=| X(ω)| (3.3-53)
The logarithm of the sum spectrum is also equal.
ln| S(ω)|=ln| X(ω)|
Theoretically, it is proved that when wavelet is the minimum phase, its logarithmic spectrum sequence (or re-matched spectrum) S (n) is a causal sequence.
Principles, methods and explanations of seismic exploration
Because ln | s (ω) | is a real even function and (n) is a real causal sequence.
Any real sequence can be written as the sum of odd and even sequences, so s (n) can be written as:
Principles, methods and explanations of seismic exploration
That is, the odd part (n) and the even part (n) of the wavelet logarithmic spectrum sequence (n) have the following two properties:
First of all, due to the causality of (n), the odd and even parts of (n) have the following relationship.
Principles, methods and explanations of seismic exploration
In the formula,
Principles, methods and explanations of seismic exploration
Secondly, the Fourier transform of even and odd parts of S (n) is the real part and imaginary part of its Fourier transform. Let the Fourier transform of (n) be (? ), (? )= (? )，S^i(？ )= (? ), (? ) is the logarithmic spectrum of wavelet, then
Principles, methods and explanations of seismic exploration
The properties obtained by Fourier transform are
Principles, methods and explanations of seismic exploration
So there is
Principles, methods and explanations of seismic exploration
that is
Principles, methods and explanations of seismic exploration
So the method of finding wavelet can be summarized as follows:
1) The real part of the reliable wavelet log spectrum is obtained by multi-channel statistical method. Wanderer spectrum
S(ω)=| S(ω)| ejφ(ω) (3.3-62)
Then there is
Principles, methods and explanations of seismic exploration
The logarithm ln | s (ω) | of wavelet amplitude spectrum is determined by the geometric average of several amplitude spectra (or the average of correlation functions of multi-channel records).
2) Find the wavelet phase spectrum φ (ω) from the logarithm of the wavelet amplitude spectrum. The calculation formula is
Principles, methods and explanations of seismic exploration
3) calculate wavelet S(t). Through | s (ω) | and φ (ω)
Principles, methods and explanations of seismic exploration
Because of the influence of interference and the incomplete correlation of reflection coefficient sequence, it is necessary to shape the amplitude spectrum and phase spectrum of wavelet. In addition, this method is only suitable for the minimum phase in theory. In order to adapt to mixed phase recording, the phase of seismic record profile can be minimized by exponential attenuation method, and then the obtained wavelet is weighted inversely exponentially.
Predictive inverse filtering
The prediction problem is to estimate the future value of a physical quantity, and get its estimated value (predicted value) at a certain time in the future by using the known past value and present value of the physical quantity. This is a very important scientific and technical issue. Weather forecast, earthquake forecast and anti-missile automatic tracking all belong to this kind of problem. Prediction is essentially a kind of filtering, which is called predictive filtering.
3.3.3.3.1principle of predictive inverse filtering
According to the prediction theory, if the seismic record x(t) is regarded as a stationary time series, the seismic wavelet b(t) is the physically realizable minimum phase signal, the reflection coefficient r(t) is uncorrelated white noise, and the seismic record x (t+α) at (t+α) is
Principles, methods and explanations of seismic exploration
Let j=s-α.
Principles, methods and explanations of seismic exploration
The first term of the analytical formula (3.3-66)
Principles, methods and explanations of seismic exploration
It can be seen that this term is determined by the future value of the reflection coefficient r(t). If the second item is
Principles, methods and explanations of seismic exploration
(t+α) is determined by t and the value of r(t) at the moment before t, that is to say, (t+α) can be predicted from the present and past data, and (t+α) is called the predicted value. Find the difference between x(t+α) and (t+α) as follows.
Principles, methods and explanations of seismic exploration
ε(t+α) is called prediction error or new record. Comparing equations (3.3-66) and (3.3-67), when the predicted value is known, the new record ε (t+α) formed by subtracting the predicted value from the original record x(t+α) is smaller than that of the original record, and the interference to the waveform after wavelet convolution is light, so the waveform is easy to distinguish, that is, the resolution is improved. In the above formula, α is called prediction distance or prediction step size. When α= 1,
Principles, methods and explanations of seismic exploration
Have it at once
Principles, methods and explanations of seismic exploration
At this time (t+ 1), there is only a constant b(0) difference between the prediction error and the reflection coefficient. Therefore, the prediction distance α= 1, and the prediction error is the reflection coefficient to achieve the purpose of inverse filtering, which is called predictive inverse filtering at this time. When α > 1, the prediction error is the result of prediction filtering. Predictive filtering is mainly used to eliminate multiples, especially sea ringing.
3.3.3.3.3.2 Calculation method of predicted value (t+α)
In predictive filtering and predictive inverse filtering, the key is to calculate the predictive value (t+α), and the method is as follows.
Inverse filtering equation
Principles, methods and explanations of seismic exploration
Expression for substituting predicted value (t+α)
Principles, methods and explanations of seismic exploration
Principles, methods and explanations of seismic exploration
Where: let τ=s-j, c(s)= b(j+α)a(s-j) be called a predictor; A(t) is the inverse filter factor; The predicted value (t+α) is the convolution of the prediction factor c(s) and the seismic record.
Now it is necessary to design an optimal prediction factor c(s) so that the predicted value (t+α) is closest to x(t+α), even if the sum of squares of prediction errors (error energy).
Principles, methods and explanations of seismic exploration
Minimum. According to the principle of least square method, that is, finding
Principles, methods and explanations of seismic exploration
Available linear equations
Principles, methods and explanations of seismic exploration
j=0， 1，2，…，m
Where: Rxx(τ) is the autocorrelation function of seismic records.
Principles, methods and explanations of seismic exploration
T is the correlation window length, and m+ 1 is the predicted value length. Write (3.3-73) in matrix form as follows
Principles, methods and explanations of seismic exploration
By solving this set of equations, the prediction filter factor c(t) can be obtained, and by convolving the seismic record x(t), the best prediction value (t+α) in the future time (t+α) can be obtained.
3.3.3.3.3 Thoughts on Inverse Filtering
Improving vertical resolution is an important task in seismic exploration. The ideal result is that the seismic wavelet is compressed into a sharp pulse and the seismic record becomes a reflection coefficient sequence. If this result can be obtained, it is equivalent to completing the inversion work. However, although many filtering methods have been developed, the practical application results are different and the universality is poor.
An important reason is that all kinds of inverse filtering methods must have some assumptions. Seismic record x(t) is the convolution result of seismic wavelet b(t) and reflection coefficient sequence r(t). Usually, in seismic exploration filtering, only the seismic record x(t) is known, but the wavelet b(t) is not known. At this time, it is impossible to find a unique r(t) from x(t). Therefore, it is necessary to make some restrictions on b(t) or r(t), that is, to find the unique optimal solution under these assumptions. The effect of inverse filtering is closely related to whether the actual situation conforms to these assumptions. For example, the least square inverse filtering and predictive inverse filtering require that the wavelet is the minimum phase and the reflection coefficient sequence is white noise. In practical work, seismic wavelets are often mixed, and the reflection coefficient sequence is not completely white noise, so it is certainly impossible to get ideal inverse filtering results. Therefore, the effort to study inverse filtering is to develop and apply inverse filtering methods whose assumptions are as close to reality as possible.
The second reason is the convolution model problem of reflecting seismic records. The seismic wavelet in convolution model is the impulse response of the earth filter. However, the geodetic filtering is very complicated, so the seismic wavelet also changes. It is generally assumed that there is a stable wavelet in inverse filtering, which has a certain gap with reality.
The third reason is the influence of noise interference. The general inverse filtering equation does not contain interference factors, so it is not a big problem for noiseless seismic records in theory. However, there is noise interference in actual seismic records. Because noise interference is random interference, its spectral characteristics are similar to reflection coefficient, and the signal-to-noise ratio of seismic records will be reduced after general filtering. Therefore, the resolution and signal-to-noise ratio should be considered in the process of inverse filtering, and the signal-to-noise ratio should not be reduced as much as possible while improving the resolution.
The fourth reason is the quality of the original seismic data. The ability to improve the resolution through processing means is limited. Using the same inverse filtering method, not all kinds of data can get the same result, but the processing effect is directly related to the original quality of seismic records. The richer the information collected on site, the better the processing effect. It can be seen from the field collection that the real high-resolution records should be collected and processed according to the requirements of high resolution.
Finally, it should be pointed out that the seismic data of each region has an optimal resolution. When improving the resolution processing, don't blindly pursue how high the resolution is, but find the best point of the resolution, and the effect is the best. If you don't handle it well, you will have an illusion.