APPLICATION OF CORRELATION ANALYSIS METHODS FOR AUTOMATIC PROCESSING OF SEISMIC DATA

V. P. Ivanchenkov, A. M. Kuvshinov, V. A. Filinova

(Presented by the seminar of the LVTA NIIYaF)

One of the most important tasks in the interpretation of seismic data is the tracing of oscillation records caused by one and the same wave from one recording point to another. This process is nothing other than correlation processing, and here a person acts as the correlation analyzer. In this case, a double correlation of seismic records is performed. This is, first, correlation in the selection of a certain group of waves on one seismic line and, second, selection from the axes of in-phase alignment over many seismograms of those that correspond to definite stratigraphic horizons. Wave correlation is one of the most labor-intensive operations; therefore the development of algorithms for automatic wave correlation is of great practical interest [1]. The large volume of information contained in seismic records, even when powerful electronic computers are used, makes it necessary to seek simplified techniques and more advanced technical methods for carrying out correlation analysis.

In the proposed article we consider an algorithm developed by us that makes it possible to carry out automatic correlation of seismic waves and, from the data of correlation processing, to construct a compressed seismic section.

The discussion is carried out for data obtained by the reflected-wave method (RWM), although the algorithm developed can also be used for processing data obtained by other methods, for example, COT.

A separate trace of a seismogram at the output of the ARU may be represented as the sum of two independent processes: a “signal” process \(Z(t)\), which is a superposition of a large number of pulses of identical form, and noise \(N(t)\) [2].

\[ X(t)=Z(t)+N(t), \tag{1} \]

where

\[ Z(t)=\sum_k a_k S(t-t_k); \]

\(S(t)\) — the shape of an individual pulse;
\(a_k\) — the amplitude of the pulse;
\(t_k\) — the arrival time.

In this case, for the set of records obtained at receiver points \(x_1=0; x_2,\ldots,x_N\), \(Z(t)\) at the point \(x\) is written as

\[ Z(t,x)=\sum_k a_k S[t-t_k-\Delta t_k(x)] , \tag{2} \]

where \(t_k\) is the arrival time of the \(k\)-th wave at the receiver point \(x_1=0\);
\(\Delta t_k(x)\) is the increment of the hodograph times of the reflected waves.
The times \(t_k\), according to experimental data [2], form a Poisson process. Let us divide the \(N\) traces of one seismogram into several sufficiently extended segments \(m\), with a time interval \(T\), within which the traces \(x_1(t), x_2(t), \ldots, x_N(t)\) may be regarded as a stationary random process.

Then an estimate of the mutual correlation function of adjacent traces can be written as

\[ B_{x_1,x_2}(\tau)=\frac{1}{2T}\int_{-T}^{T} x_1(t)\cdot x_2(t+\tau)\,dt . \tag{3} \]

The essence of the algorithm adopted by us is as follows. In the course of processing, pairwise computation is performed of estimates of the mutual correlation functions \(B_{x_1,x_2}(\tau)\), \(B_{x_2,x_3}(\tau), \ldots, B_{x_{N-1},x_N}(\tau)\), with appropriate filtering of the traces for the first segment of the seismogram. Subsequently, the positions of the principal maxima of the mutual correlation functions of the traces \(B_{x_i,x_j}(\tau)\) are sought, and the time shifts \(\tau_m^{ij}\) at which \(B_{x_i,x_j}(\tau)\) attain their extremal value are determined. The computed values \(\tau_m^{ij}\) are an estimate of the mean shift \(\Delta t_{k,\mathrm{av.}}^{ij}\) between traces for the selected segment \(m\):

\[ \Delta t_{k,\mathrm{av.}}^{ij}=\tau_m^{ij}+\xi, \tag{4} \]

\(\xi\) is the error in computing the mean shift \(\Delta t_{k,\mathrm{av.}}\).

As shown in [2], finding the mean shift between traces from the maximum of the mutual correlation function is an asymptotically optimal procedure for estimating the indicated shift. To increase the reliability of computing \(\tau_m^{ij}\), when determining the mutual correlation functions a restriction is introduced on the change in the lag \(\tau\), determined by the maximum shift of waves on the seismogram:

\[ |\tau|\leq \rho . \tag{5} \]

For the first segment of the seismogram, directed summation of the traces is carried out in accordance with the computed value of \(\tau\), i.e., the summed signal is determined for definite time shifts \(\tau_m^{ij}\):

\[ F(t)=\sum_{i=1}^{N} x_i(t+\tau_m^{ij}) . \tag{6} \]

Then a search is made for the maximum excursions of the summed signal at some prescribed threshold \(n\), and their time positions are determined. The indicated procedure is repeated for subsequent segments of the seismogram, and a tying of the obtained resul-

... in accordance with the number of the segment \(m\). When tying individual segments, summation of traces for each of them is carried out not over the entire time interval \(T\), but over an interval \(T' < T\). In this case \(T'\) may be chosen, for example, from the relation

\[ T = T' + t_b, \tag{7} \]

where \(t_b\) is the visible period of oscillation of the seismic wave.

The indicated criterion excludes loss of information when tracing waves at the ends of a segment.

It should be noted that the mean shift between traces, determined in our case from the cross-correlation functions, may depend on several unknown parameters. Thus, at small angles of dip of the horizon it will determine the static and the averaged kinematic correction. The latter makes it possible to process a seismic record without preliminary introduction of static and kinematic corrections.

Fig. 1. Segment of the seismogram being processed

In order to determine the effectiveness of the proposed algorithm, it was studied on a BESM-4 electronic computer using materials from the Tomsk Geophysical Trust.

In testing the algorithm, the task was set of estimating the accuracy of determining the mean shift \(\tau_m^{ij}\) from the principal maximum (minimum) of the cross-correlation function of seismic traces, and of determining the possibility of further using the results of correlation analysis to identify regular waves by means of directed summation (multiplication) of seismic records.

For processing, segments of seismograms of profile \(72 \cdot 3 \cdot 14\) were selected, including 8 traces on which a group of oscillations belonging to one wave was traced (Fig. 1). Digitization of the seismic traces was carried out with a sampling interval of 1.67 msec. Table 1 gives

1. Call IS-2

2. Conversion of the initial data from \(10 \to 2\)

3. Formation of cells for computation

4. Reading seismic information from magnetic tape

5. Transfer of information for computing \(\Sigma, \Pi\)

6. Computation

\[ B(z)=\frac{1}{T_0}\int x_i(t)\,y_j(t+z)\,dt, \qquad i=1,2,\ldots,j,\quad N|z|<\tilde c \]

1. Determination of \(\tilde{\tau}^{ij}_m\)

2. Computation

\[ K(\tau)=\frac{1}{T_0}\max \int [x_i(t)+y_j(t+\tau)]\,dt \]

1. Formation of \(\tilde{\tau}^{ij}_m\)

2. End of seismic record?

no — return to reading seismic information from magnetic tape.

yes — proceed further.

1. Reading information for computing \(\Sigma, \Pi\)

2. Computation

\[ \Sigma=\sum_i f_i\!\left(t+\sum_{j=1}^{n}\tau^{ij}_m\right) \]

1. Computation

\[ \Pi=\prod_{1}^{n} f_k\!\left(t+\sum_{i=j}\tau^{kj}_m\right) \]

1. Printing plots of \(\Sigma\) and \(\Pi\)

2. Program stop

Fig. 2. Program flowchart

statistical characteristics for eight seismic traces, computed on the BESM-4.

To study the proposed algorithm, a program was developed, the block diagram of which is shown in Fig. 2. Figure 3 presents the cross-correlation functions of eight seismic traces, computed with an averaging interval \(T\) equal to 0.83 sec. The same figure

Table 1

\(m\) — mathematical expectation; \(\sigma^2\) — variance;
\(\sigma\) — root-mean-square deviation.

shows the plot of the autocorrelation function \(B_{x_1,x_1}(\tau)\), the position of the principal maximum of which determines the start of the count of the mean shifts \(\tau_m^{ij}\). The computed values of \(\tau_m^{ij}\) for the seismogram segment are given in Table 2. The summed signal \(F(t)\), obtained for the wave shifts \(\tau_m^{ij}\) given in Table 2, is shown in Fig. 4. Comparison of the values of \(\tau_m^{ij}\) and \(t_0\), obtained on the computer and by measurement under a microscope, showed their good agreement. In the course of processing the algorithm, the resolution of waves under directed multiplication of seismic records was investigated. In this case \(t_0\) was determined from the position of the maximum values of the signal

\[ G(t)=\prod_{i=1}^{N} x_i(t+\tau_m^{ij}). \tag{8} \]

Fig. 3. Plot of correlation functions.

The plot of the function \(G(t)\) is shown in Fig. 5. From a comparison of Figs. 4 and 5 it is evident that the shape of the regular seismic wave under directed multiplication of trace records is substantially distorted; however, a considerably better suppression of noise waves is observed than in summation. On the computer, a check was also carried out of the possibility of applying approximate estimates of the cross-correlation functions of seismic traces for carrying out phase tracking. In this case the computation of approximate estimates of the cross-correlation functions was performed according to the following expressions:

\[ R_{x_i,x_j}(\tau)=\frac{1}{T}\int_{0}^{T}\max\left[x_i(t),\bigl(x_j(t+\tau)\bigr)\right]dt; \tag{9} \]

\[ R_{x_i,x_j}(\tau)=\frac{1}{T}\int_{5}^{T}\left|x_i(t)-x_j(t+\tau)\right|\,dt, \tag{10} \]

The computation of functions (9), (10) can also be implemented quite simply by optical correlation for records of seismic signals on transparencies by the variable-area method and the “wide trace” method. Such methods of recording are widely used in analog processing of seismic information. The values of \(\tau_m^{i,j}\) computed by expression (9) for one segment of the seismogram are given in Table 3.

Fig. 4. Plot of the summed signal \(F(t)\)

Table 2

A comparison of the results obtained in computing the “exact” and approximate estimates of the cross-correlation functions showed the possibility of using them in the automatic processing of seismic-

seismic data. The use of functions (9), (10) in the correlation processing of data on an electronic digital computer makes it possible to substantially reduce the expenditure of machine time. The performance of approximate correlation also proves preferable in optical processing. Thus,

Fig. 5. Graph of the function \(G(t)\)

the studies carried out on the materials of seismic profile 72 3 14 showed the good resolving power of the proposed algorithm and the possibility of using it in the mass processing of seismic information. It should be noted that, in the processing—

Table 3

in seismic materials with substantial curvature of the phase hodograph, according to this algorithm it is advisable to carry out a preliminary introduction into the records of calculated kinematic and static corrections. In this case, the introduction of calculated corrections makes it possible to increase the reliability of determining the mean shift from the position of the maximum of the cross-correlation functions of seismic traces.

References

1. I. I. Gurvich. Seismic Exploration. Nedra, 1970.

2. A. K. Yanovsky. Use of cross-correlation functions of seismic traces for tying seismic observations. — In: Discrete Correlation of Seismic Waves. Novosibirsk: Nauka, 1971.