Method for Stratigraphic Correlation and Reflection Character Analysis of Setsmic Signals

1. Field of the Invention

The invention relates generally to the field of digital seismic signal processing and, more particularly, but not by way of limitation, it relates to improved methods for pattern recognition as applied to the classification and analysis of seismic signals, particularly seismic section data.

2. Description of the Prior Art

There are various prior procedures for examination of seismic data in order to establish meaningful data indications. All of the prior approaches to seismic energy analysis, whether effected manually, on specialized processing machines, or on digital or analog computers under control of any of various software dictates, all utilize either visual or linear correlation of data groups or indications in order to establish meaningful data findings. In general, the prior methods each embody an inadequacy due to the fact that the spacings of stratigraphic reflectors will change non-linearly from region to region on a seismic section. Thus, the prior methods attempted to operate with a built-in or inherent defect in their operation. Applicant is not aware of any prior art that attempts to characterize seismic reflection characteristics using a stochastic modeling procedure.

The present invention relates to a method for examining seismic data within a statistical decision framework in order to distinguish selected areas of seismic data on the basis of reflection character. Thus, the procedure uses stochastic models to characterize seismic reflection characteristics within a given training area of interest and, thereafter, this reflection characteristic data can be compared with like data from another seismic area of interest in order to establish a degree of similarity or other distinguishing criteria. The technique may establish a stochastic model of either noncausal or causal type depending upon how the seismic signal at a sample point is predicted based on the information from its neighborhood.

Therefore, it is an object of the present invention to develop stochastic seismic signal data models for subsequent use in a probablistic framework for such as performing stratigraphic correlations and delineation of anomalous structure on seismic section data.

It is also an object of the invention to provide a statistical analysis method for the detection and identification of reflection continuity in the presence of data loss on the section.

Also, it is an object of the invention to provide a method for segmentation or classification of different groups of reflection characteristics on a seismic section.

Finally, it is an object of the invention to classify seismic reflection data on a statistical basis without regard to visual or linear correlation of data.

Other objects and advantages of the invention will be evident from the following detailed description when read in conjunction with the accompanying drawings which illustrate the invention.

FIG. 1 is a flow diagram of the general steps of the process of the present invention;

FIG. 2 is a schematic illustration of a seismic section region with its neighboring periods;

FIG. 3 is a schematic drawing of a sample point of a seismic section region showing its two-dimensional neighborhood in surround;

FIG. 4 is a flow diagram illustrating computation of the noncausal model parameter estimation procedure;

FIG. 5 is a flow diagram illustrating the sequence of steps for scanning the rows of a seismic reflection data region;

FIG. 6 illustrates a region of a seismic section as it is utilized in characterization of the causal model;

FIG. 7 is a black and white reproduction of a seismic section illustrating correlation of reflection packages across a salt dome;

FIG. 8 is a black and white reproduction of a photograph illustrating the delineation of the salt dome with reflection character identification;

FIG. 9 is a black and white reproduction of a photograph illustrating segmentation of seismic reflection data into plural packages of different reflection characteristics; and

FIG. 10 is a black and white production of the seismic section illustrating reflection character identification by image segmentation.

FIG. 1 illustrates the method of the present invention for stratigraphic correlation and reflection character analysis. A problem common to the objectives is to classify a test region 10 into one of the given n types of training region 12, 14 through 16. The region characterization scheme of FIG. 1 illustrates the basic procedures wherein it is desirable to classify the character of test region 10 relative to plural selected training regions 12, 14 and 16.

The reflection signal characteristics of test region 10 as well as training regions 12-16 are stochastically modeled and a set of parameters are estimated for each region. The test region 10 is classified into one of the training types, i.e. regions 12-16, by using a decision scheme, either (i) that based upon maximum likelihood concepts or (ii) that based upon distance between the distributions of the data in the test and training regions. Both noncausal and causal models may be used for region characterization.

Within a seismic section there may be any number of different anomalies of interest and these can generally be broken down into plural training regions such as regions 12-16 which exhibit particular characteristics visually. Such characteristics which it may be desirable to classify or relate more broadly either over the entire seismic section or to individual test regions 10. The training regions 12-16 are then stochastically modeled on a probability basis to estimate parameters for each individual training region 12-16 as at estimate stages 18, 20 and 22. In like manner, the data of test region 10 would also undergo estimation of parameters, as at stage 24, and the test parameters would then be compared to the groups of parameters from the training regions 12-16 to determine most similarity. Thus, the stage 26 would determine similarities on either a maximum likelihood or minimum distance decision scheme, as will be further described. The decision will then output data identifying the particular region type as at output stage 28.

Region parameter estimation may be effected using a non-causal model such as that of FIG. 2. FIG. 2 illustrates a region 30 having neighboring periods as shown. Let us allow that the region 30 contains N.sub.1 columns and N.sub.2 rows and that it is periodic in nature.

The signal value at a sample point of the region 30 is modeled as a weighted combination of signal values in a two-dimensional window around the sample point and an additive noise term. Let .theta. contain the weights and the noise is taken as Gaussian with zero mean and unit variance. The noise is independent and identically distributed from sample to sample. Let Y be the row stacked version of the elements of the region 30 and let u be the row stacked version of noise elements. Using the periodic property of the region we can write

where .beta. is a constant that specifies the variance of the noise. The matrix A(.theta.) is block circulant and of size N.sub. N.sub.2 .times.N.sub.1 N.sub.2. That is ##EQU1##

The matrix elements A.sub.j of A(.theta.) are circulant and of size N.sub.1 .times.N.sub.1. Thus, we can write ##EQU2##

Referring to FIG. 3, a sample point of the region at the position (0,0) and its rectangular neighborhood of size K.sub.1 .times.K.sub.2 are shown. This amounts to a sample point 32 of the region 30 and its two-dimensional neighborhood. Let K.sub.1 and K.sub.2 be odd, and let K.sub.12 =K.sub.1.vertline.2 and K.sub.22 =K.sub.2.vertline.2 be the integers obtained by truncation. The first rows of circulant matrix elements A.sub.O, A.sub.1, . . . , A.sub.N.sbsb.2.sub.-1 of the block circulant matrix A(.theta.) are shown in FIG. 4. The indices of the parameter .theta. are given in FIG. 3. The elements of the index vector IND are used to form a weighted combination of sample values in the window. The elements of the index vector IND are identified in FIG. 4 and would be equated to specific matrix elements as shown in Table I.

With u being Gaussian, the likelihood function of Y given the parameters .theta. and .beta. can be written as ##EQU3## where .lambda..sub.k are the eigenvalues of the block circulant matrix A(.theta.). .lambda..sub.k are the functions of the parameter vector .theta.. .OMEGA. denotes the region being considered, and s specifies the central point of a two-dimensional neighborhood, and x(s) is a vector of sample values formed by using the index vector IND. The eigenvalues .lambda..sub.k are computed by exploiting the block circulant properties of A(.theta.). That is, they are the fourier coefficients of the elements of the first row of A(.theta.). The maximum likelihood estimates of the parameters and are obtained by iteratively maximizing Equation (4). The initial values of .theta. and .beta. can be determined using least squares approximation. They are ##EQU4##

The computations consists of scanning the rows of the region by setting up a two-dimensional array SD whose row dimesion is equal to the row dimension of the two-dimensional window. The first K.sub.1 rows of the regions are transferred into SD after taking the periodicity of the data into account. Thereafter, the computations are repeated by adding a row to and deleting a row from SD until all of the rows of the region have been scanned. The computations for model parameter estimation using noncausal models are summarized in FIGS. 4 and 5.

Region characterization may also be estimated with causal models derived from a region of seismic section in trace format as shown in FIG. 6. Thus, a region 36 is defined by selecting a time portion of a plurality of successive traces of a seismic section. Letting x(i,t) be the tth sample of an ith trace, and letting x(i,t), be an estimate of x(i,t). Then, the estimate x(i,t) may be written as ##EQU5## where f.sub.j is an lx1 vector of weights, l is the number of channels and n is the number of samples of each channel. The data vector X(i,t) of l channels at time t is given by ##EQU6## The error e(i,t) in the estimation of x(i,t) can be written as

Let e(i,t) be Gaussian with mean zero and variance s. Let it also be independent from sample to sample. The the variation in the data essentially arises because of the variation in e(i,t). The likelihood function of the data, given the parameters f.sub.j, n, and s, can be written as ##EQU7## Taking the logarithm of Equation (10), we get ##EQU8## Differentiation of Equation (11) with respect to s and equating the resulting expression to zero yields ##EQU9## Using Equation (12) in Equation (11), it can easily be shown that the maximization of Equation (11) is equivalent to the minimization of .epsilon..sup.2 given by the following: ##EQU10## Differentiation of Equation (13) with respect to the coefficient vectors f.sub.j and equating the resulting expressions to zero yields ##EQU11## Equations (14) can be efficiently solved using a multichannel generalization of the Levinson algorithm. Substitution of (14) into (13) yields an expression for the minimum value of .epsilon..sup.2 as ##EQU12## From Equations (15) and (16), we get

From Equation (12), s can be written as ##EQU13##

Region recognition may be effected as at stage 26 (FIG. 1), and the test region can be classified into one of the training region types either on the basis of the maximum likelihood decisions or the distance between distributions of data. Thus, based on the maximum likelihood decision, the test region may be identified. The likelihood that the data from a selected test region comes from a distribution of the kth training region with the parameters .theta..sub.k and .beta..sub.k can be written as ##EQU14## Where .lambda..sub.k l are eigenvalues of A(.theta..sub.k) and .OMEGA. is the kth training region. Taking the logarithm of Equation (19), the kth discriminate function DF.sub.k becomes ##EQU15## where ##EQU16## The quantity DF.sub.k serves as a measure of similarity between the kth training region and the test region. The test region is recognized as being of the type of Kth training region if

A test region can also be classified as similar to a specific training region using causal modeling based on distance between data distributions in the regions. Thus, the distance J.sub.k can be computed between the data distribution of Kth training region and the test region fom the expression ##EQU17## where (.theta..sub.k, .beta..sub.k) are the estimated parameters of the kth training region and (.theta..sub.t, .beta..sub.t) are the estimated parameters of the test region. Also, d.sub.k .sbsb.l are the eigenvalues of A(.theta..sub.k) and d.sub.t.sbsb.l are the eigenvalues of A(.theta..sub.t). The test region is recognized as the type of training region K for which ##EQU18##

In the above, the divergence information measured between the distribution of the data in each of the training regions and the test region is used as a measure of similarity or distance between the distribution in recognizing the test region type as one of the training regions. Treating the distribution of the data in the training and test regions as Gaussian with zero mean, we can then write

and

where N stands for normal distribution and

and

For zero mean Gaussian distribution, the divergence information between the kth training region and the test region can be written as

Let W be a matrix that diagonalizes the block circulant matrices A(.theta..sub.k) and A(.theta..sub.t) of the training and test regions.

Where D.sub.k is the diagonal matrix having the eigenvalues of A(.theta..sub.k) as its diagonal elements. We also have ##EQU19## From Equations (28) and (29), we can write ##EQU20##

Similar relations exist for block circulant matrices A(.theta..sub.t). From Equations (26) and (27)-(30), we obtain ##EQU21## where d.sub.k.sbsb.l and d.sub.t.sbsb.l are the diagonal elements of D.sub.k and D.sub.t, respectively. Substitution of Equation (31) into Equation (27) then yields Equation (23).

FIGS. 7-10 illustrate some test results that were achieved using the pattern recognition of the present invention. The testing was performed on compressional wave seismic data profiles selected from a region having anomalies that particularly benefit from the evaluation scheme of the present invention. Thus, the basic data of FIGS. 7 and 8 show a subsurface region having a distinct salt dome that divides the subsurface region. The seismic data of FIGS. 9 and 10 illustrates data from a different overall section that is similar in slope and character to subterrain adjacent a dome structure.

The first application of the present reflection character analysis technique is illustrated in FIG. 7 wherein it is the aim to detect continuation of specific horizons from one side to the other of the salt dome in the presence of loss of data due to the salt dome. Thus, a training region 40, as outlined in crayon was, selected as a region of particular interest and the procedure was directed toward establishing that same reflection character in a test region over on the left hand side of the salt dome 42. Thus, three test regions 44, 46 and 48 were selected as possible candidates for the continuity of reflection horizons in the training or template region 40. A non-causal stochastic model was first fit to the template or training region 40. Thirty trace by thirty sample area test regions 44, 46 and 48 were then statistically modeled. The divergence distance measure was used as a measure of similarity between each test subregion 44, 46 and 48 and the template region 40, and the location of maximum correlation was identified as indicated by a cross in test region 44. After initial identification of the general region, more careful examination verified the conclusion that the continuity of reflection horizons was correctly identified.

Another application of the present method is illustrated by the photo reproduction of FIG. 8 which shows the delineation of salt dome 42 and reflection character identification in color overlay of one of four classes of reflection character. Four training regions were selected for their four types of distinctly different criteria. One such region is selected to represent the salt dome, while three remaining different character training regions were identified. The fcolor differentiation of the photograph is in four shades of green or gray tone as shown. Non-causal stochastic models were fit to the training regions. The seismic section was then divided into subregions of size 25 traces by 25 samples and the statistical models for each of the four training regions were compared to each of the subregions to establish best fit by computation of the measure of similarity between each subregion and the most similar training region. It may readily be seen that the salt dome 42 is fairly accurately delineated as to reflection character and other areas of specific reflection character may also be classified as at darkest section 50 on the left of the dome and an intermediate green section 52 on the right, each of these being predominately of one shade of indicator.

Referring to FIG. 9, the technique of the present invention is applied to carry out segmentation of the seismic section 54 into individual reflection packages of certain characteristics. Four training regions were initially selected and the seismic section was divided into subregions of size 25 traces by 25 samples. Stochastic models were then determined for the training and test regions, and test subregions were identified into one of the training region types having maximum similarity. Thus, FIG. 9 illustrates a seismic section segmented into packages of different reflection characteristics, four different classes in all.

The seismic image of FIG. 10 illustrates the same seismic section as classified into three different types of reflection characteristics. Here again, we have a segmentation operation wherein similarity of characteristics was determined over a different training area. In FIG. 10 the test and training regions were selected in a basic size of 75 traces by 25 samples, this contributing to the elongated, rectangular general pattern.

The foregoing discloses a novel method of seismic data interpretation using machine processing which employs basic pattern recognition techniques in the classification and correlation of seismic data. Use of the present method brings still another tool to the geophysicist for use in interpreting subsurface data representations, and the application of pattern recognition techniques enables identification of new parameters and type comparisons between various forms of data.

Changes may be made in combination and arrangement of elements as heretofore set forth in the specification and shown in the drawings; it being understood that changes may be made in the embodiments disclosed without departing from the spirit and scope of the invention as defined in the following claims.