Spectral Data Processing Method for Detection of Hydrocarbons

This invention relates to a method for processing spectral data for the earth's surface such that the results of such data processing can be used to predict the presence of hydrocarbons.

It is known to collect and process spectral data for the earth's surface for the purpose of locating hydrocarbons. According to one prior method, intensities for different wavebands of electromagnetic radiation reflected from the earth's surface are detected by an appropriate means, such as the LANDSAT satellite. This intensity data is then processed to give alteration indexes which can be used to predict the presence of hydrocarbons in certain areas.

Although the above-described prior method provides excellent results, improvement and further development would be desirable.

It is, therefore, an object of the invention to provide an improved spectral data processing method for identifying areas which are potentially productive with respect to hydrocarbons.

The object is realized by a method of processing spectral data for a preselected set of pixels which correspond to adjacent areas of the earth's surface, wherein the method comprises: selecting n waveband(s) as represented by the positive integers 1, . . . , n; obtaining intensity values for each waveband with respect to electromagnetic radiation reflected from the areas of the earth's surface corresponding to each pixel so as to yield an intensity data set {x.sub.1, . . . , x.sub.n } for each pixel, where x.sub.1 . . . , x.sub.n corresponds to the intensity values for respective wavebands 1, . . . , n; dividing said intensity data sets into a predetermined number of disjoint data set clusters, wherein any one intensity data set is contained in only one cluster; determining the value of each of monotonic functions f.sub.1 (x), . . . , f.sub.n (x), where x represents intensity, for each respective waveband for each pixel so as to yield a function data set {f.sub.1 (x.sub.1), . . . , f.sub.n (x.sub.n)} for each pixel, wherein f.sub.1 (x.sub.1), . . . , f.sub.n (x.sub.2) are hereafter denoted as y.sub.1, . . . , y.sub.n ; selecting certain productive pixels and dry pixels, associated with each of the clusters, previously determined to correspond to hydrocarbon-bearing and nonbearing areas of the earth' s surface respectively; determining a discriminant probability function, being a function of y.sub.1, . . . , y.sub.n, for each cluster by utilizing the function data sets corresponding to the productive pixels and dry pixels associated with each respective cluster; and determining for each cluster the value of its corresponding discriminant probability function with respect to each function data set, the value of the discriminant probability function being indicative of the probability of the presence of hydrocarbons. Pixels can be determined to be either potentially productive or dry with respect to hydrocarbons based on their corresponding discriminant probability function values. Additionally, it is within the scope of certain aspects of the invention to plot a map based on probability values for the pixels. Such a map enables easy and fast identification of areas with high production potential.

The previously determined productive and dry pixels can be, by way of example, pixels containing productive and dry wells respectively. By using spectral data corresponding to these pixels, the discriminant probability function so obtained can be employed as an effective exploration tool in identifying potentially productive areas, as is further explained in the following description and as is demonstrated in a subsequent example.

FIG. 1 is a map of an area for which spectral data was processed in accordance with the invention. Letters A-U indicate areas which have high discriminant probability function values associated therewith and are therefore considered to be good production prospects. The map also shows stream drainage with solid dark lines and producing trends with dashed lines.

FIG. 2 is an enlarged depiction of one area (township) shown in the map of FIG. 1. This Figure is a contour map representative of discriminant probability function values as will be discussed further below.

FIG. 3 is a graph which correlates the number of dry wells and successful wells with various probability values obtained according to the invention.

FIG. 4 is a graph which illustrates the chi-squared test with respect to the invention which, as is discussed further below, shows the statistical significance of results obtained in accordance with the invention.

FIG. 5 is a flow chart illustrating the various steps performed in accordance with a preferred embodiment of the invention.

A preferred embodiment of the invention will now be described. First however, it is appropriate at this point to define several terms which are used repeatedly throughout the following description.

A "pixel" as used herein is an element of an image that corresponds to an area of the earth's surface, typically square in shape.

A "waveband" is range of wavelengths in the electromagnetic spectrum.

A "cumulative probability density function (CPDF)" is a function such that for an arbitrary parameter Z associated with a number of observations (corresponding to a number of pixels, for examples) CPDF (Z) for a particular value of Z is the fraction of the total number of observations having parameter values less than the particular value of Z.

A "data set" is groups of data elements which are listed hereafter by enclosing them in braces.

A "mean data set" is that data set derived by averaging the data elements of the various data sets in a group of data sets.

A "vector" is a means by which a data set can be expressed. For example, the data set {x.sub.1, . . . , x.sub.n } is expressible as vector x which has components x.sub.1, . . . , x.sub.n. Any vector subsequently discussed is assumed to be a column vector, unless stated otherwise, which can be defined as an n.times.b 1 matrix having n rows and one column.

The "transpose" of a vector is denoted by a prime symbol ('), which indicates interchanging of rows and columns. The transpose of a column vector (n.times.1 matrix) is a row vector (1.times.n matrix).

The various steps of an embodiment of the invention are as follows. Some of the steps may be performed in a different order if desired.

1. Select Waveband(s)

It is within the scope of the invention to select any number, including only one, and/or combination of wavebands since the presence of hydrocarbons will have at least some effect on the spectral reflective properties of the earth's surface with respect to virtually any waveband. The number of wavebands selected is hereafter denoted as n such that the wavebands are represented by the positive integers 1, . . . , n. However, it is preferable to use at least two different wavebands, and it is most preferable to use at least three different wavebands. It is also desirable that at least one of the wavebands is a green waveband in the wavelength range of about 500 nm to about 600 nm since the reflective properties of the earth's surface have been found to be highly affected for this waveband. Employing the green waveband in combination with a red waveband in the wavelength range of about 600 nm to about 700 nm is particularly advantageous in giving results which closely correlate with the presence of hydrocarbons. Such a combination of wavebands along with a third infrared waveband (i.e. wavelength greater than 700 nm) is demonstrated to be particularly effective in a subsequently described example.

2. Obtain Intensity Values for Each Waveband

Spectral data, as hereafter described in detail, is collected with respect to a preselected section of the earth's surface which can be divided for the purpose of spectral data collection into a plurality of pixels. This section of the earth's surface can contain from, for example, about 500 to about 100,000 pixels, and more preferably from about 30,000 to 50,000 pixels.

Most conveniently, spectral data for each pixel is collected by means of one of the LANDSAT satellites. One type of LANDSAT satellite uses a multispectral scanner which has the capability of detecting radiation in several different wavebands, including a green waveband in the wavelength range of about 500 nm to about 600 nm, a red waveband in the wavelength range of about 600 nm to about 700 nm, and infrared wavebands in the wavelength ranges of 700-800 nm and 800-1,100 nm.

In operation, the scanner aboard LANDSAT scans the earth's surface so as to receive radiation being reflected from each of the separate pixels. For each pixel, the scanner produces an electrical signal proportional to the energy reflected from the pixel for each waveband. These signals are then transmitted in digital form to a ground station for recording. The signal produced and transmitted as described is generally proportional to the intensity of radiation received by the multispectral scanner for each pixel. However, the recorded signal is divided into 64 unitless increments, such that a signal for a particular pixel is assigned a digit from 0 to 64. The parameter thus detected and measured by the LANDSAT scanner is usually called "brightness", which is a nonphysical parameter, but is generally proportionally related to intensity. Therefore, throughout this description, the term intensity will be regarded as synonymous with the term brightness.

LANDSAT data can be obtained in digital or image form and can be purchased from the EOSAT Corp. at the EROS Data Center, Sioux Falls, S.D. If data is purchased in image form, a photographic representation of the terrain is supplied which is typically converted to digital form for processing. Most conveniently for the purposes of this invention however, data is obtained in digital form on computer-compatible tapes. Typically, as purchased, a LANDSAT tape will contain data for what is called a "scene". Such a scene is a relatively large area of the earth's terrain which includes about 7 millions pixels. Thus, only a portion of section of the terrain is usually selected from the scene whose associated data is processed according to the present method.

Although the preferred embodiment is described as employing LANDSAT data, the above data could alternatively be collected by means of a scanner aboard an aircraft similar to that described in U.S. Pat. No. 4,421,981 of Hough.

By using the LANDSAT data for the section of interest, an intensity value is obtained for each selected waveband with respect to each pixel of the section so as to yield an intensity data set {x.sub.1, . . . , x.sub.n } for each pixel where x.sub.1, . . . , x.sub.n corresponds to the intensity values for respective wavebands 1, . . . , n. It should be noted that these integers as they are arbitrarily assigned to particular wavebands are solely for convenience of notation to differentiate between the wavebands. An integer 1, . . . , or n can be arbitrarily assigned to each of the n wavebands. For example, for two selected red and green wavebands, the integer 1 could be assigned to the green waveband and the integer 2 could be assigned to the red waveband, or vice versa.

3. Divide Intensity Data Sets Into Disjoint Clusters

The various intensity data sets are now divided into a preselected number, preferably 3 to 5, of disjoint clusters such that any one intensity data set is in only one cluster. It is most advantageous to group data sets corresponding to pixels with similar spectral properties into a particular cluster. The final clusters generally correspond to pixels having similar terrain. For example, one cluster may correspond to areas of scrub oak, another cluster may correspond to pasture, and yet another may correspond to heavy vegetation along river or creek banks. Therefore, the effects hydrocarbons have on vegetation and/or soil in any particular cluster is relatively uniform among the various pixels so as to enable derivation of statistically significant discriminant probability functions for each cluster.

A preferred procedure (FASTCLUS) for clustering data sets is described below. It should be understood, however, that any other clustering technique is within the scope of the invention.

The FASTCLUS procedure is discussed in detail on pages 377-400 of the SAS.RTM. User's Guide: Statistics, Version 5 Edition, published by the SAS Institute Inc. of Cary, N.C., copyrighted in 1985. Software employing this procedure is available from the SAS Institute. The FASTCLUS procedure, as carried out in accordance with a preferred embodiment of the invention, performs a disjoint cluster analysis on the basis of Euclidean distances. By way of example, Euclidean distance between two data sets and their associated vectors .alpha. (having components .alpha..sub.1, . . . .alpha..sub.n) and .beta. (having components .beta..sub.1, . . . .beta..sub.n) can be defined as .sqroot.(.alpha..sub.1 -.beta..sub.1).sup.2 + . . . (.alpha..sub.n -.beta..sub.n).sup.2. For the purpose of the following description of FASTCLUS, the term "distance" is understood to mean Euclidean distance.

By way of example and for the sake of illustration of FASTCLUS, assume that four clusters are desired. Four intensity data sets are randomly selected. These randomly selected data sets are called "seeds". The distance between each remaining intensity data set and each respective seed is then determined and each data set is assigned to the nearest seed. Assignment of a data set to a particular seed means that data set is a member of the cluster associated with that seed. Each data set is assigned to only one cluster. The various intensity values associated with the data sets in each cluster are now averaged to derive a means data set which replaces the seed. Distance from each remaining data set to the mean data sets are determined to decide which data sets are closest to which means such that the data sets are reassigned to the nearest mean. This procedure is repeated as many times as is necessary until the means either do not move anymore, or as is more typical, move less than a predetermined distance. At this final state of FASTCLUS, the final four clusters according to this illustration are formed by assigning each data set to the nearest mean such that the data set is a member of the cluster associated with that nearest mean.

4. Determine the Value of Functions for Each Pixel

The next step of the method is to determine the value of each of functions f.sub.1 (x), . . . , f.sub.n (x) for each respective waveband for each pixel in each cluster so as to yield a function data set {f.sub.1 (x), . . . , f.sub.n (x.sub.n)} for each pixel in each cluster. The function f(x) should be monotonic. That is, the function should either increase or decrease in value as the intensity value x increases. The functions f.sub.1 (x), . . . , f.sub.n (x) may be the same or different, but are most conveniently the same, and are most preferably cumulative probability density functions (CPDF's) for each cluster. Using this particular type of function is advantageous primarily because it changes in value rapidly around the mean intensity data set of each cluster. Thus with respect to a particular cluster, and a particular intensity data set in that cluster, CPDF (x.sub.1) is the fraction of the total number of pixels in that cluster which have intensity values for the first waveband less than a particular intensity value x.sub.1. The function is applied similarly to intensities for the other wavebands and for other pixels in the cluster, and function data sets are obtained for the remaining clusters.

For ease of notation f.sub.1 (x.sub.1), . . . , f.sub.n (x.sub.n) will hereafter be denoted y.sub.1, . . . ,y.sub.n.

5. Select Productive and Dry Pixels

Certain pixels are now selected from each of the clusters which have previously determined to be hydrocarbon-bearing and nonbearing respectively. It is preferable to select at least one such productive or dry pixel for each approximately six to ten square mile area, and most preferably for each three square mile if possible.

Such pixels can be determined to be either dry or productive by employing any reliable technique. Of course, if the information is available it is most desirable to assign certain pixels as being productive or dry based on the presence of productive or dry wellsites located with the respective pixels. However if such information is not available, other means of selecting such pixels can be used, such as geochemical techniques. One particularly suitable geochemical technique is known as MOST (Microbial Oil Survey Technique) which is described in U.S. Pat. No. 2,880,142 of Hitzman. Certain pixels could be assigned to be "productive", or hydrocarbon-bearing, or "dry", nonbearing, based on geochemical survey results for such pixels.

5. Determine Discriminant Probability Function for Each Cluster

A discriminant probability function (being a function of y.sub.1, . . . , y.sub.n which represents f(x.sub.1), . . . , f(x.sub.n) as noted above) is now determined for each cluster by utilizing the function data sets which correspond to the productive and dry pixels associated with each respective cluster. The value of such a function for a pixel should be indicative of the probability of the presence of hydrocarbons in that pixel. Although any probability function having this property is within the scope of the invention, a function obtained by the DISCRIM procedure is presently preferred.

The DISCRIM procedure is described on pages 317-333 of the SAS.RTM. User's Guide, supra, which is hereby incorporated by reference. This procedure as applied to a particular cluster develops a discriminant probability function corresponding to that cluster by making use of a measure known as generalized square distance. The notation below is used to described generalized squared distance:

Covariance matrices are well known in the field of statistics (i.e. multivariate analysis) and are discussed in, for example, the text Probability and Statistics by Morris H. DeGroot, Addison-Wesley Publishing Co., 1975, pp. 530-532, such pages being hereby incorporated by reference. A covariant matrix of y (having variable components y.sub.1, . . . , y.sub.n) is defined to be the n.times.n matrix such that for i=1, . . . , n and j=1, . . . , n, the element in the ith row and jth column is the covariance of y.sub.i with respect to y.sub.j, or Cov (y.sub.i,y.sub.j). The concept of covariance is discussed in the Degroot text on pages 172-173. A "within group" covariant matrix is the covariant matrix constructed from data sets in a particular group, whereas the pooled covariant matrix is the covariant matrix derived by averaging the elements of the p and d within group matrices.

The generalized square distance from a vector y to, for example, the mean vector m.sub.p is

where

if the within-group covariance matrix is used, or

if the pooled covariance matrix is used; and

if the prior probabilities are not all equal, or

if the prior probabilities are equal. As discussed above, the prior probabilities are assumed to be equal, so that g.sub.2 (p)=0 and

The option of using within-group or pooled covariance matrices can be built into a suitable computer program (later set forth) for performing the various data processing steps. Of course, the generalized squared distance for group d is determined similarly.

The discriminant probability function for a particular cluster and for the productive group can be expressed in terms of generalized squared distance as follows: ##EQU1## The value P of this function for a particular vector y corresponding to a particular function data set {y.sub.1, . . . , y.sub.n } is the probability that the pixel corresponding to this data set bears hydrocarbons.

6. Determine Value of Discriminant Probability Function With Respect to Each Pixel

For each cluster, the value P of the corresponding discriminant probability function is determined with respect to each function data set and its corresponding pixel. Thus, a value P is assigned to each pixel.

7. Determine the Value of CPDF(P) for Each Pixel

It is preferable to determine the cumulative probability density function CPDF(P) for each pixel in each cluster rather than relying solely on the P values themselves. By determining CPDF(P) for each pixel, the result is a fractional value corresponding to each pixel which essentially ranks any particular pixel with respect to the other pixels. For example, a CPDF(P) value of 0.90 for a particular pixel means that 90% of the total number of pixels have P values less than the P value for this particular pixel.

8. Map CPDF(P) Values

A particularly effective but optional technique for visually illustrating the CPDF(P) values for pixels in a section involves mapping the CPDF(P) values by plotting on a map, of the areas corresponding to the pixels, locations having associated therewith pixels with CPDF(P) values above a predetermined value. Most preferably, different ranges of CPDF(P) values can be plotted so as to be visually distinguishable using various black and white notations or different colors. Such mapping is discussed and illustrated in a subsequent example.

9. Select Pixels as Potentially Productive Based on CPDF(P) Values

Certain pixels are now selected as potentially productive if their corresponding CPDF(P) values are above a predetermined value, which is preferably greater than 0.05. These selected pixels are now targeted as drilling prospects and can be further investigated with other hydrocarbon locating techniques if desired.

This example is provided to further illustrate the invention and should not be construed to limit the invention in any manner.

LANDSAT spectral data was obtained for a section in Kansas which is represented by the map in FIG. 1. The section contains 540,000 pixels and is divided by solid lines into 6 mile.times.6 mile square areas ("townships"), such as the one shown at 10. Intensity values were obtained for each pixel with respect to a green waveband (500 nm-600 nm), a red waveband (600 nm-700 nm) and an infrared waveband (800-1,100 nm), thus yielding three intensity values x.sub.1, x.sub.2 and x.sub.3 for each pixel. Pixels in LANDSAT imagery correspond to areas of about one acre. To increase the pixel size to approximately twenty acres, the image was divided into 5.times.5 pixel blocks, and the intensity values x.sub.1, x.sub.2, x.sub.3 for all the pixels in each block were averaged together to produce average intensities which characterized the larger twenty acre area as one single pixel.

This intensity data was processed in accordance with a preferred embodiment of the invention, described previously, by employing four clusters (derived with FASTCLUS) and CPDF functions for f.sub.1 (x.sub.1), f.sub.2 (x.sub.2) and f.sub.3 (x.sub.3) to derive the function data sets. In addition 258 dry pixels and 223 productive pixels were selected based on the presence of either dry or productive wellsites in the respective pixels. A discriminant probability function of the form of equation (7) was determined for each cluster using the DISCRIM procedure, and probability P values were computed for each pixel by utilizing these functions and within-group covariant matrices. Finally, CPDF(P) values were determined for each pixel. Of course, all of the data processing was performed by a computer using a computer program which is presented following this example.

Referring again to FIG. 1, "productive trends" are indicated by dashed lines as shown at 12. These trend lines are drawn so as to lie along actual producing fields. Drainage streams (mostly small creeks) are represented by solid lines as at for example 14. Letters A-U indicate areas having associated therewith CPDF(P) values of at least 0.60.

A and C-J correspond to producing fields and are used for control. It can be seen that these lie on or near producing trends. K-U are sites which have not been drilled, but are excellent prospects based on their high CPDF(P) values. K-R can be seen to lie on producing trends. S and T do not lie on producing trends, but are on circular drainage anamolies. That is, these are surrounded by drainage streams. Such an area is usually associated with a topographic high and corresponding subterranean structure (i.e. anticline) favorable to the presence of hydrocarbon deposits. U is not on a producing trend or drainage anamoly but corresponds to a very large area with high (above 0.06) CPDF(P) values, thus making it an attractive prospect. B corresponds to a producing field which does not lie on a trend or drainage anamoly. The fact that this area is indicated to be a good prospect according to the invention particularly demonstrates the effectiveness and usefulness of the invention. B is indicated by the invention to be likely to bear hydrocarbons, which in fact it does; and moreover the area is not indicated by conventional analysis (considered here) to be a good prospect.

Referring to FIG. 2, there is shown a contour plot of area 10 in which areas with CPDF(P) values within certain ranges are plotted. As shown, areas with CPDF(P) values of 0-0.60 are illustrated as white background, areas with values of 0.60-0.90 are cross-hatched, and areas with values of 0.90-1.00 are stippled (dotted). Reference character 16 denotes a particularly large area with high CPDF(P) values which has at its center smaller areas with the highest CPDF(P) values of 0.90-1.00. This area is thus considered a very good prospect and corresponds to P in FIG. 1.

Referring to FIG. 3, there is shown a graph which clearly illustrates the effectiveness of the invention in predicting the presence of hydrocarbons. The horizontal axis corresponds to CPDF(P) values, and the vertical axis corresponds to the percentage of the total number (491) of wells in section 10. For example, for a CPDF(P) value of 0.0125, about 2% of the total number of wells which are located in pixels having this CPDF(P) value are dry holes, whereas only about 0.2% are productive. As the CPDF(P) values increase, it can be seen that the productive percentage generally increases. Thus, the CPDF(P) value produced by the present invention is clearly correlated with the presence of hydrocarbons.

A table is set forth below which gives the prediction results obtained for the wellsite pixels (those containing wells) if an arbitrary threshold level of 0.50 is selected for CPDF(P). That is, if one of the wellsite pixels has a CPDF(P) value of above 0.50, it is classified as productive.

Thus, the method of the invention predicts that of the 491 total wellsites, a total of 248 will be productive and 243 will be dry. Of the 248 predicted to be productive, 146 are actually productive, but 102 are dry. Dividing 146 into 248 gives a success rate of about 58% with respect to predicting whether or not certain sites will be productive. This compares favorably with the actual success rate of only 45%, derived by dividing the total number of actual productive wells (223) into the total number of wells drilled (491). Thus, success rates could be improved up to 13% simply by screening drilling targets in accordance with the present invention.

The statistical significance of the data in the above Table can be shown using the chi-squared test. The chi-squared test is a well known test in the field of statistics which measures the extent to which a contingency table, like the table above, departs from random chance. The larger the value of chi-squared, the larger the departure from randomness. A discussion of chi-squared can be found in, for example, the text Statistical Methods, sixth edition, by Snedecor and Cochran, Iowa University Press, 1967, on pages 215-220. The chi-squared test was carried out by performing a series of 51 separate runs in which the 491 wellsites were assigned to be either productive or dry on a purely random basis.

The results of these random runs as compared with the results obtained according to the invention are shown in FIG. 4. The points as shown at 18 indicate the chi-squared values obtained with the random runs, and the rectangle at 20 indicates the chi-squared value obtained by application of the invention to real field data. These values are plotted versus CPDF(.chi..sup.2), where .chi. is the symbol chi, and where the CPDF(.chi..sup.2) values for a particular run indicates the fraction of runs having smaller chi-squared values. The graph shows that the results of the invention as set forth in the above table give a chi-squared value of about 26.2 as compared with much lower chi-squared values for the random runs. This shows the statistical significance of the results obtained according to the invention, and proves that these results were not obtained due to random chance.

Referring now to FIG. 5, there is shown a flow chart which outlines the various steps of a preferred embodiment of the invention. The invention, as has been demonstrated above, provides a technique of processing spectral data which effectively predicts the presence of hydrocarbons in selected areas.

Obviously many modifications and variations of the present invention are possible in light of the above teachings. It is therefore to be understood that within the scope of the appended claims the invention may be practiced otherwise than as specifically described.

There immediately follows a preferred computer program for carrying out the various data processing steps of the invention. This program is written for an IBM 3090 computer and is self-explanatory to one skilled in the use of the IBM 3090 computer. The input required into the computer program is the imagery data over an area (section) of interest, photogrammatic and geographic coordinates of specific locations on the earth's surface; geographic coordinates of a set of points on the earth's surface known to be productive; similar coordinates for a set of points known to be nonbearing; certain special mapping data such as the geographic coordinates of specific section, township and range boundaries. ##SPC1## PG,48