Amplitude processing and processing methods to improve signal-to-noise ratio and resolution

In seismic data processing, maintaining the true amplitude characteristics of seismic waves and improving the signal-to-noise ratio and resolution of seismic records as much as possible are called "three high" processing. This has always been the pursuit of seismic data processing personnel. goal. Because the quality of the "three high" processing directly affects the accuracy and effect of lithological parameter extraction and seismic exploration.

10.3.1 True amplitude recovery

Maintaining the true amplitude characteristics of seismic waves (referred to as amplitude-preserving processing) should include two major aspects in a broad sense: that is, true amplitude recovery (or called amplitude preservation processing) Amplitude compensation) and other amplitude maintenance issues in various processes. This section mainly discusses the methods of true amplitude recovery. For other processing methods that affect the amplitude characteristics, corresponding measures must be taken to keep the relative relationship of the amplitudes unchanged as much as possible.

After gain recovery processing, the amplitude characteristics of seismic records are consistent with the amplitude characteristics of seismic waves received by surface geophones. This amplitude is still not called true amplitude. What we call true amplitude refers to the amplitude of the reflected wave generated by the difference in formation wave resistance, that is, the amplitude that can reflect the changes in formation lithology. In addition to the changing factors of formation wave impedance, the amplitude received at the surface also has spherical diffusion factors and inelastic attenuation factors. Therefore, it is necessary to eliminate the effects of spherical diffusion and inelastic attenuation and restore the true amplitude characteristics of seismic waves.

Spherical dispersion is the attenuation of amplitude due to wavefront expansion as the wave propagates away from the source. Such amplitude attenuation (A) is inversely proportional to the propagation distance r

Exploration Technology Engineering

where v is the average velocity of the overlying medium at the interface; t is the recording time of the reflection. Correction for spherical diffusion requires multiplying the data by the time-varying function vt.

Inelastic attenuation is the result of elastic wave energy being dissipated into heat and absorbed by the formation due to internal friction when it propagates in the rock. The principle section has explained that this attenuation is an exponential function of frequency and propagation distance

Exploration Technology Engineering

where α is the inelastic attenuation coefficient (absorption coefficient)

Exploration Technology Engineering

Therefore, inelastic attenuation can be corrected by multiplying the data by the function eαvt. At this point, the true amplitude recovery process is completed.

The coefficient α can be determined from the amplitude-time function after gain recovery and spherical diffusion correction. In order to obtain a good statistical estimate of α, a set of seismic traces are used to measure the energy and obtain the attenuation curve.

There is another method of true amplitude recovery, in which velocity information is not required. After gain recovery, the amplitude decay is assumed to be an exponential function. Therefore, by fitting the gain-corrected records with an exponential function according to the least squares method, the true amplitude correction function (that is, including both spherical diffusion and inelastic attenuation corrections) is obtained.

As mentioned above, the wavefront divergence factor K is

Exploration Technology Engineering

where r and t are the propagation distance and propagation time of the wave respectively, C and a are constants related to formation velocity.

The absorption attenuation factor is

Exploration Technology Engineering

where α is the absorption coefficient; b is a constant to be determined. The total impact of wavefront divergence and absorption attenuation is

Exploration Technology Engineering

The method for obtaining a and b is as follows.

Read the amplitude extreme value (peak or trough) of the reflected wave from the seismic record, and use (10.3-5) as the regression equation, we get

Exploration Technology Engineering

In the formula: ut=lnAi-lnti; Ai, ti are the amplitude extreme values ??and their corresponding times; N is the number of points of the amplitude extreme values. The correction function is a-1tebt.

In order to obtain representative true amplitude recovery parameters, the selected seismic traces should have no multiple sweeps and a high signal-to-noise ratio. For areas with stable geological conditions, a set of parameters can represent the entire area. When the geological conditions in the work area change significantly, these parameters need to be recalculated.

10.3.2 Digital filtering to improve signal-to-noise ratio

In seismic exploration, the seismic waves used to solve geological tasks are called effective waves, while other waves are collectively called interference waves. Suppressing interference and improving signal-to-noise ratio is a task that runs through the entire process of seismic exploration. In addition to taking corresponding measures to suppress interference in field data collection, digital filtering is also a very important measure to improve the signal-to-noise ratio in digital processing of seismic data.

The digital filtering method uses the difference in frequency and apparent speed between the effective wave and the interference to suppress the interference, which are called frequency filtering and apparent speed filtering respectively. And because frequency filtering only needs to operate on a single channel of data, it is called one-dimensional frequency filtering. Implementing apparent velocity filtering requires processing multiple channels of data at the same time, so it is called two-dimensional apparent velocity filtering. This section mainly introduces these two filtering methods.

10.3.2.1 One-dimensional frequency filtering

The so-called one-dimensional digital filtering refers to the use of computers to filter single-variable signals. The single variable can be time or frequency, or it can be Space or wave number. Taking time or frequency as an example to discuss one-dimensional digital filtering, other principles are the same.

(1) Principle of one-dimensional digital filtering

Suppose the seismic record x(t) is composed of effective wave s(t) and interference wave n(t), that is,

Exploration Technology Engineering

The spectrum is

Exploration Technology Engineering

Where: X(f) is the spectrum of x(t) ;S(f) and N(f) are the spectra of s(t) and n(t) respectively. If the amplitude spectrum of X(f) |X(f)| can be represented by Figure 10-6. It shows that the amplitude spectrum of the effective wave |S(f)| is in the low frequency band, while the amplitude spectrum of the interference wave is in the high frequency band.

Figure 10-6 Schematic diagram of spectrum distribution of effective waves and interference waves

If the amplitude spectrum of a frequency domain function H(f) is designed to be |H(f)|,

Exploration Technology Engineering

The graph is shown in Figure 10-7(a).

Order

Exploration technology engineering

and

Exploration technology engineering

In the time domain there are ( Using the convolution theorem of Fourier transform)

Exploration Technology Engineering

H(f) is called the one-dimensional filter frequency response, and (10.3-9) is the frequency domain The filter equation, h(t) is the time domain function of H(f), which is called the one-dimensional filter filter factor (Figure 10-7(b)). (10.3-11) is the time domain filtering equation, y(t) and Y(f) are respectively the seismic records and spectra with only effective waves after filtering, φx(f), φy(f), φh(f) are respectively Seismic records before and after filtering and the phase spectrum of the filter. The above filtering mainly uses the frequency difference between the effective wave and the interference wave to eliminate the interference wave, so it is also called frequency filtering.

Figure 10-7 Filter frequency response and filter factor

The above-mentioned filter is called ideal low-pass filter. According to the different frequency band distribution of effective wave and interference wave, it can also be Filters are divided into ideal band-pass filters, ideal high-pass filters, etc. The so-called ideal means that the frequency response of the filter is a rectangular gate. The effective waves inside the gate pass without distortion, which is called the passband, and the interference waves outside the gate are all eliminated. This is actually not possible in digital filtering. Because the filter factors that can be processed during digital filtering can only be finitely long, while the filter factors of an ideal filter composed of discontinuous functions are infinitely long. In practical applications, it can only be truncated to a finite length. After truncation, a truncation effect will occur, that is, the frequency response corresponding to the truncated filter factor is no longer an ideal rectangular gate, but a line close to a rectangular gate, but with amplitude fluctuations. Curve, this phenomenon is called the Gipps phenomenon.

Since the frequency response curve is a fluctuating curve in the passband, the effective wave will definitely be distorted after filtering. In addition, the curve outside the passband is also fluctuating, which must not effectively suppress interference. To avoid the Gipps phenomenon, several methods can be used, one of which is the flanging method.

10.3.2.2 Two-dimensional apparent velocity filtering

(1) Proposal of two-dimensional apparent velocity filtering

In seismic exploration, sometimes effective waves and interference waves The spectrum components are very close or even overlap. In this case, frequency filtering cannot be used to suppress the interference. It is necessary to use other differences between the effective wave and the interference wave for filtering. If there is a difference in apparent velocity distribution between the effective wave and the interference wave, apparent velocity filtering can be performed. This kind of filtering needs to calculate several channels at the same time to get the output, so it is a two-dimensional filtering.

The seismic wave received by the surface is actually a two-dimensional function g(t, x) of time and space, that is, a combination of vibration diagram and wave profile, between the two through

Exploration technology and engineering

are intrinsically linked. In the formula, k is the spatial wave number, which represents the number of wavelengths per unit length, f is the frequency, which describes the number of vibrations per unit time, and v is the wave speed.

In actual seismic exploration, observations are always made along the ground survey line. The above wave numbers and velocities should be substituted with the wave number component kx and the apparent velocity v*. Then there are

Exploration Technology Engineering

Since the seismic fluctuation is a two-dimensional function of the space variable x and the time variable t, and there is a close relationship between space and time, no matter which one is performed separately Filtering in one dimension will cause changes in the characteristics of another dimension (for example, frequency filtering alone will change the shape of the wave profile, and wavenumber filtering alone will affect the vibration pattern and produce frequency distortion), producing undesirable effects. Then only by forming a time-space domain (or frequency wavenumber domain) filter based on the intrinsic relationship between the two, can the purpose of suppressing interference and highlighting effective waves be achieved. Therefore, two-dimensional filtering should be performed.

(2) Principle of two-dimensional apparent velocity filtering

The principle of two-dimensional filtering is based on the two-dimensional Fourier transform. The seismic wave g(t,x) observed along the direct measurement line on the ground is a wave that changes with time and space. Its frequency wave number spectrum G(ω,kx) and space-time function are obtained through two-dimensional forward and inverse Fourier transform. .

Exploration Technology Engineering

The above formula shows that g (t, x) is composed of countless plane simple harmonics with angular frequency ω = 2πf and wave number kx. They propagate along the survey line with apparent speed v*.

If the plane simple harmonic components of the effective wave and the interference wave are different, and the plane harmonic component of the effective wave propagates at a different apparent speed from the plane harmonic component of the interference wave, as shown in Figure 10-8, then They can be separated by two-dimensional apparent velocity filtering to suppress interference and improve the signal-to-noise ratio.

(3) Calculation of two-dimensional filtering

Figure 10-8 Propagation of effective waves and interference waves with plane simple harmonics of different components

Two-dimensional linear The properties of the filter are determined by its space-time characteristics h(t,x) or frequency-wavenumber characteristics H(ω,kx). Like same-dimensional filtering, in the time-space domain, two-dimensional filtering is implemented by the two-dimensional convolution operation of the input signal g (t, x) and the filter operator h (t, x). In the frequency -In the wavenumber domain, it is completed by multiplying the spectrum G(ω,kx) of the input signal and the frequency wavenumber characteristic H(ω,kx) of the filter.

Exploration Technology Engineering

Due to the discrete nature of seismic observations and the limited arrangement length, the summation of a limited number (N) of record traces must be used to replace the spatial coordinates integral.

Exploration Technology Engineering

In the formula, n is the original track number; m is the result track number.

It can be seen from equation (10.3-15) that two-dimensional convolution can be reduced to the summation of the results of one-dimensional convolution. Therefore, the result of two-dimensional filtering at any point on the survey line can be obtained by adding the one-dimensional filtering results of N seismic traces. At this time, each channel is processed by its own filter, and its time characteristic hm-n(t) depends on the distance between the channel and the output channel. Calculating along the survey line in sequence can obtain the two-dimensional filtering results on the entire survey line (Figure 10-9).

Same as ideal one-dimensional filtering, ideal two-dimensional filtering also requires that the amplitude spectrum of the frequency-wavenumber response within the pass band is 1, and is 0 outside the pass band, and the phase spectrum is also 0, that is, Zero phase filtering. Therefore, the frequency-wavenumber response of a two-dimensional ideal filter is a positive real symmetric function (two-dimensional symmetry, that is, symmetrical to both parameters), and the space-time factor must be a real symmetric function. Two-dimensional filtering also has pseudo-gate phenomena and Gipps phenomena, which can also be solved by the flanging method and the multiplication factor method. Because it is a two-dimensional function, the situation is much more complicated. Usually, only the method of reducing the sampling interval (including the time sampling interval Δt and the frequency sampling interval Δf) and increasing the number of calculation points (including the number of points M and N in both the time and space directions) is adopted. method.

Figure 10-9 Schematic diagram of two-dimensional filtering calculation (N=5)

(4) Fan filter

The most commonly used two-dimensional filter is fan filter. It filters out low-speed and high-frequency interference. Its frequency wave number response is

Exploration Technology Engineering

Figure 10-10 Frequency wave number response of the sector filter

The pass band is on the f-kx plane A sector-shaped area starting from the coordinate origin and symmetrical with the f-axis and kx-axis is formed (Figure 10-10). Therefore this filter is called a sector filter.

Using the inverse Fourier transform, the factor can be obtained as

Exploration Technology Engineering

When the operation is implemented on a computer, discretization is required. For time sampling: t=nΔ, n=0, ±1, ±2,..., Δ is the time sampling interval, Δ=1/2fc. The spatial sampling interval is the track spacing Δx of the input track.

A standard sector filter can be used to construct a pie-cut filter that suppresses both high apparent speed interference and low apparent speed interference. Furthermore, a band-shaped filter that suppresses high and low frequency interference at the same time can also be constructed. Pass fan filter and bandpass pie filter.

Apply sector filtering before superposition. The suppression target can be surface waves, scattered waves, refracted waves or waves generated by cable vibration. As for applications after superposition, multiple reflections or side waves generated from inclined interfaces can be suppressed.

10.3.3 Inverse filtering process to improve longitudinal resolution

According to the propagation theory of seismic waves, seismic waves propagate underground in the form of seismic wavelets in the medium. The reflected wave seismic record received by the ground is the convolution of the formation reflection coefficient and the seismic wavelet. Therefore, the formation acts as a filter, turning the reflection coefficient sequence into a seismic record composed of wavelets, reducing the longitudinal 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 formation filtering and improve the longitudinal resolution of the seismic records.

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

Exploration technology engineering

Due to the wavelet problem, 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. Let the inverse filter factor be a(t), and require a(t) and b(t) to satisfy the following relationship

Exploration Technology Engineering

Use a(t) to inversely filter the seismic record x(t)

Exploration Technology Engineering

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

In the specific implementation of inverse filtering, the core is to determine the inverse filtering factor a(t). Due to the uncertainty of seismic wavelets and the presence of noise interference in seismic records, it is very difficult or even impossible to determine the exact a(t) in practice.

To this end, many methods for determining the inverse filter factor a(t) have been studied under different approximate assumptions. These methods can basically be divided into two categories: one is to first obtain the seismic wavelet b(t), Then find a(t) based on b(t); another major category is to find a(t) directly from earthquake records. There are many different methods in each category (there are only so many inverse filtering methods, which illustrates the difficulty of inverse filtering processing). Several representative inverse filters among the inverse filtering methods are discussed below.

10.3.3.1 Formation inverse filtering

The formation inverse filtering is a method of first finding the wavelet b(t) and then finding a(t). This method requires well logging data and good wellside seismic records. First, the acoustic logging data is converted into a formation reflection coefficient sequence r(t) that matches the wellside seismic record x(t). The frequency domain equation can be obtained by calculating the spectrum of r(t) and x(t): p>

Exploration Technology Engineering

That is

Exploration Technology Engineering

where B (ω) is the wavelet b (t) Spectrum, and then the relationship between wavelet and inverse filter factor is:

Exploration Technology Engineering

A(t) is obtained by inverse Fourier transform. In the formula, A (ω) is the spectrum of the inverse filter factor. Written as z transformation, it is A(z)=. It can be seen that A(z) is a rational fraction. To make A(z) stable, the root of the denominator polynomial B(z) must be outside the unit circle, that is, a wavelet is required. b(t) is the minimum phase.

Use well logging and side-well seismic traces to obtain the wavelet and inverse filtering factors, that is, use the inverse filtering factors to inversely filter other traces of the survey line.

10.3.3.2 Least squares inverse filtering

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

The basic idea of ??least squares inverse filtering is to design a filter operator and use it to convert the known input signal into a given expected output signal that is optimally close in the sense of least square error output.

Suppose the input signal is x(t), which is convolved with the filter factor h(t) to be found to obtain the actual output y(t), that is, y(t)=x(t)*h (t). Due to various reasons, the actual output y(t) cannot be exactly the same as the pre-given expected output (t), and the two can only be required to be optimally close. There are many standards for judging whether the two are optimally close, and the least square error criterion is one of them. That is, when the sum of the squares of the errors between the two is the minimum, it means that the two are optimally close. In this sense, the filtering performed by finding the filter factor h(t) is the least square filter.

If the filter factor to be found 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 squares principle, the inverse filtering factor obtained when the sum of the squares of the errors is minimized is called the least squares inverse filtering factor. The inverse filtering used for seismic records x(t) is least squares inverse filtering.

Suppose the input discrete signal is seismic wavelet b(n)={b(0), b(1),...,b(m)}, and the inverse filtering factor to be found is a(n)= {a(m0), a(m1), a(m2),..., a(mm)}, m0 is the starting time of a(t), (m+1) is a( The continuation length of t), the convolution of b(n) and a(n) is the actual output y(n), that is,

Exploration Technology Engineering

is the seismic wavelet and The cross-correlation function of the desired output.

According to the least squares principle, the basic equation of least squares inverse filtering can be obtained through derivation:

Exploration technology engineering

Where,

Exploration technology engineering

is the autocorrelation function of the seismic wavelet,

Exploration technology engineering

is the interaction between the seismic wavelet and the desired output Related functions.

Equation (10.3-24) is a linear system of equations, written in matrix form:

Exploration Technology Engineering

The symmetry of the autocorrelation function is used in the equation sex. In this equation, the coefficient matrix is ??a special positive definite matrix, called the general Tobritz matrix. This matrix equation can be quickly solved by Levinson's recursive algorithm.

Equation (10.3-27) adapts wavelet b(n) to the minimum phase, maximum phase and mixed phase. In the formula, the starting time m0 of the inverse filtering 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 filtering factor.

10.3.3.3 Predictive inverse filtering

The prediction problem is to estimate the future value of a certain physical quantity, using the known past value and present value of the physical quantity to obtain its value at a certain future date. The problem of estimated value (predicted value) at a moment. It is a very important issue in science and technology. Weather forecasting, earthquake forecasting, automatic tracking of anti-missiles, etc. all belong to this type of problem. Prediction is essentially a kind of filtering, called predictive filtering.

(1) Principle of predictive inverse filtering

According to the prediction theory, if the earthquake record x(t) is regarded as a stationary time series, the seismic wavelet b(t) is a physical The minimum achievable phase signal, the reflection coefficient r(t) is uncorrelated white noise, and based on the convolution model of the seismic record, the seismic record x(t+α) at ??(t+α) is

Exploration Technology Engineering

Analyze the first term of (10.3-28)

Exploration Technology Engineering

It can be seen that this term is caused by reflection Determined by the future value of coefficient r(t). If the second term is

Exploration Technology Engineering

^x (t+α) is determined by the value of r (t) at the time before t and t, that is to say ( t+α) can be predicted from current and past data, and (t+α) is called the predicted value. Find the difference between x (t+α) and (t+α) as

Exploration Technology Engineering

ε (t+α) is called the prediction error, or new Record. Comparing the two formulas (10.3-28) and (10.3-29), when the predicted value is known, the new record ε(t+ α) involves less reflection coefficients than in the original record, the degree of interference with the waveform after convolution is light, and the waveform is easy to distinguish, that is, the resolution is improved.

In the above formula, α is called the prediction distance or prediction step size. When α=1,

Exploration technology and engineering

There is

Exploration technology and engineering

At this time (t+1 ) time, there is only a constant b(0) difference between the prediction error and the reflection coefficient.

Therefore, the prediction distance α=1 is selected, and the prediction error is the reflection coefficient, which achieves the purpose of inverse filtering. This is called prediction inverse filtering.

When α>1, the prediction error is the result of prediction filtering. Prediction filtering is mainly used to eliminate multiple waves, especially to eliminate offshore tremors.

(2) Method of calculating predicted value (t+α)

In predictive filtering and predictive inverse filtering, the key is to calculate the predicted value (t+α). The method is as follows .

From the inverse filtering equation

Exploration Technology Engineering

Put it into the expression of the predicted value (t+α)

Exploration Technology Engineering Learn

Let τ=s-j, c(s)=b(j+α)a(s-j) be called the predictor. a(t) is the inverse filter factor. The predicted value (t+α) is the convolution of the predictor c(s) with the seismic record.

Now we need to design an optimal predictor c (s) so that the obtained predicted value (t+α) is closest to x (t+α), even if the sum of squares of the prediction errors (error energy )

Exploration technology engineering

is the smallest. According to the least squares principle, a system of linear equations can be obtained

Exploration Technology Engineering

where Rxx (τ) is the autocorrelation function of seismic records

Exploration Technology Engineering

T is the length of the relevant time window, and m+1 is the length of the predictor. Write (10.3-34) in matrix form as

Exploration Technology Engineering

Solving this set of equations, we can obtain the prediction filter factor c(t), which can be used to record the earthquake x (t) Convolution can find the best predicted value (t+α) at ??the future time (t+α).