Reliability measures for statistical prediction of geophysical and geological parameters in geophysical prospecting

Abstract

A method for assessing reliability of a prediction model constructed from N training data attribute vectors and N associated observed values of a specified parameter. Each training data attribute vector includes seismic attributes obtained from seismic data traces located at or near a well and each associated observed value is obtained from well log or core data from the well. A predicted value of the specified parameter is determined for each of the N training data attribute vectors, from the training data attribute vectors and the prediction model. A residual is determined for each of the N training data attribute vectors, as the difference between the associated observed value and the predicted value of the specified parameter for the training data attribute vector. An attribute vector for the designated location is determined. The predicted value of the specified parameter at the designated location is determined, from the attribute vector for the designated location and the prediction model. N basic probability distributions are determined from the N training data attribute vectors, the N associated observed values, the N residuals, and the predicted value. N basic probability assignments are determined for each of three hypotheses that the predicted value is reliable, unreliable, and unpredictable, respectively, from the N basic probability distributions. A reliability value, an unreliability value, and an unpredictability value for the predicted value are determined as combinations of the N basic probability assignments for each of the three hypotheses.

Description

FIELD OF THE INVENTION

[0002] This invention relates generally to the field of statistical prediction models used in geophysical prospecting. In particular, the invention relates to a method for assessing reliability of predicted values for geophysical or geological parameters from statistical predictive models used in geophysical prospecting.

BACKGROUND OF THE INVENTION

[0003] In the ongoing search for subsurface hydrocarbons, geophysicists have developed methods for evaluating and interpreting the structure and characteristics of the subsurface formations of the earth. One such method is the analysis of seismic data traces. Of particular importance to geophysicists is the ascertainment of formation structure and characteristics consistent with the presence of hydrocarbon deposits.

[0004] The seismic data traces that are analyzed using the method of the present invention are obtained by any conventional means. One of ordinary skill in the art would understand that seismic data traces are usually obtained by the reflection of acoustic waves from geologic layers of differing acoustic impedance. Acoustic impedance is dictated by the physical properties of the material through which the acoustic wave travels. Such properties include lithology, porosity, and fluid content.

[0005] Seismic data are generally obtained by imparting seismic energy into the earth, then receiving and recording the energy reflected by the subsurface formations. This received energy is then processed to produce seismic signals or traces, which depict reflection information as a function of the time lapse from signal generation to reception and the embedded seismic pulse. As is known to those of ordinary skill in the art, processing of seismic data may vary, but typically includes stacking, migration, and deconvolution.

[0006] Originally, seismic data traces were used simply for ascertaining formation structure. However, exploration geophysicists have developed a number of methods to obtain a variety of characteristics that describe the seismic data traces. Such characteristics are termed attributes. These attributes provide quantitative measures of the shape of the seismic data traces. Attributes are said to be instantaneous when values for the attributes are obtained for each data point (i.e., each time sample) or within a small time window of data points. Examples include amplitude, phase, frequency, dip, and power. Attributes are said to be interval, when values for the attributes are taken from a specified seismic interval within a seismic trace. Examples include averages, maximums, and cycle widths of measured characteristics of the seismic signal over the seismic interval.

[0007] The objectives of seismic attribute analysis include identifying the boundaries of geologic intervals of differing acoustic impedance and assigning characteristic values, which may exhibit lateral variations, for the physical rock properties within each of these geologic intervals. There are two key steps used in the art to assign these characteristic values for the physical rock properties. The first step is a seismictosynthetic well tie, which compares the seismic signal shapes and attributes identified on at least one seismic data trace at or near a well location with those identified in a synthetic seismogram which is developed for that well. The well being used is termed a calibration well. This synthetic seismogram is generated using well data such as well log data or core data, coupled with standard techniques familiar to those skilled in the art. The second step is termed seismic attribute calibration, which involves statistically relating the attributes obtained from seismic data traces that are presumed to represent the formation properties at a well for any seismic interval, with the measured rock properties in that well over that same interval. Thus, a seismictosynthetic well tie relates a real seismic data trace to a synthetic seismogram, while calibration relates a real seismic data trace to actual rock properties as determined by the well data, such as well log or core data.

[0008] In seismic attribute calibration, a statistical prediction model is used to investigate a relationship between observable seismic attribute variables and a variable representing a geological or geophysical rock property. This involves establishing a predictive relationship between the independent variables observed at densely distributed seismic locations and the dependent variable of interest observed at typically more sparsely distributed well locations. The statistical prediction model may be linear regression or any other statistical technique. This relationship may then be used to predict the variable of interest at locations away from the wells by using the observed seismic attributes.

[0009] Suppose the relationship between the rock property variable and attribute variables can be described by a statistical prediction model. In general, let z be the dependent variable of interest and let x represent the independent variables as an attribute vector. It is assumed that z and x are related by a functional relationship, z=f(x). A statistical prediction model represents the relationship as {circumflex over (z)}={circumflex over (f)}(x), where {circumflex over (f)} is an estimating function for the prediction model and {circumflex over (z)} is the predicted value for a given attribute vector x. This relationship will be illustrated using a multiple linear regression model, although the present invention applies to any statistical prediction technique. In a multiple linear regression model, the dependent variable z is represented as a linear combination of the independent attribute vector x:

z=a T x+b+e.   (1)

[0010] Here a is a vector of linear coefficients, b is an intercept, and e is a model or observation error. The superscript T represents the transpose of a vector.

[0011] For a set of N training data points, (xi, zi), i=1, . . . , N, of independent attribute vectors xi and dependent variable zi, the classical regression approach uses a least squares criteria to estimate the linear coefficients a and b. This yields a prediction model:

{circumflex over (z)}=â T x+{circumflex over (b)},   (2)

[0012] where {circumflex over (z)} is the predicted value for attribute x, and â and {circumflex over (b)} are prediction coefficients.

[0013] Reliability measures of statistical prediction models such as confidence intervals can be calculated only when certain statistical assumptions are imposed on the training data, which are rarely met in real situations. If the expected value of the error term e of equation (1) is zero, then the prediction coefficients â and {circumflex over (b)} of prediction equation (2) are unbiased estimates of the linear coefficients a and b of equation (1). This assumption also implies that the underlying relationship between the variables x and z in equation (1) is linear. In this case, the predicted value {circumflex over (z)} is an unbiased estimate of z. However, if this linear assumption does not hold, then the predicted value {circumflex over (z)} becomes a biased estimate of z. There are cases where the physical relationship between the dependent variable and the independent variables supports a linear regression model. In these cases, it would be useful to have statistical confidence intervals for the predicted value {circumflex over (z)} in prediction equation (2). A confidence interval gives upper and lower bounds between which there is a given probability, say 95%, of finding the variable. However, obtaining these confidence intervals requires that the errors that populate the error term e in equation (1) be normally distributed, have constant variance, and be statistically independent of one another. These conditions are rarely met in practical data analysis situations, and the confidence intervals computed from a multiple linear regression model will usually give false information. The confidence intervals associated with multiple linear regression analysis are valid measures for assessing the reliability of predicted values only when linearity and the above statistical assumptions are met.

[0014] Similarly, for a predictive neural network, such as a back-propagation neural network, no standard method for determining reliability exists. Thus, there exists a general need for quantifying the reliability of predicted values from statistical prediction models without the strict assumptions needed in classical confidence interval calculations.

SUMMARY OF THE INVENTION

[0015] The present invention is a method for assessing the reliability associated with a statistical prediction of a specified geological or geophysical parameter at a designated location. The statistical prediction is obtained from a prediction model constructed from a set of N training data points. Each training data point comprises a training data attribute vector and an associated observed value of the specified geological or geophysical parameter. Further, each training data attribute vector includes one or more seismic attributes obtained from seismic data traces located at or near a well and each associated observed value of the specified parameter is obtained from well log or core data from that well. First, a residual is determined for each of the training data points. The residual is the difference between the associated observed value of the specified parameter for the training data point and a predicted value of the specified parameter for the training data point obtained from the prediction model. Next, an attribute vector is determined for the designated location. Next, a predicted value of the specified parameter at the designated location is obtained using the attribute vector for the designated location and the prediction model. Next, N basic probability distributions are determined from the N training data attribute vectors, the N associated observed values, the N residuals, and the predicted value of the specified parameter at the designated location. Next, N basic probability assignments are determined for each of three hypotheses that the predicted value of the specified parameter at the designated location is reliable, unreliable, and unpredictable, respectably, from the N basic probability distributions. Finally, reliability, unreliability, and unpredictability values for the predicted value of the specified parameter at the designated location are determined as combinations of the N basic probability assignments for each of the hypotheses that the predicted value of the specified parameter at the designated location is reliable, unreliable, and unpredictable, respectively.

BRIEF DESCRIPTION OF THE DRAWINGS

[0016] The present invention and its advantages may be more easily understood by reference to the following detailed description and the attached drawings in which:

[0017] FIG. 1 shows the three reliability measures (reliability, unreliability and unpredictability) for a linear prediction from ten hypothetical training data points, in which the attribute vector consists of only one element;

[0018] FIG. 2 is a flowchart illustrating the processing steps of one embodiment of the method of the present invention for determining reliability measures for prediction values for a specified geophysical or geological parameter from a statistical prediction model;

[0019] FIG. 3 is a flowchart illustrating a preferred embodiment of the method of the present invention for optimizing the prediction error level γ, standard deviation σ, and maximum prediction distance dmax parameters;

[0020] FIG. 4 is a flowchart illustrating a preferred embodiment of the method of the present invention for calculating the reliability measures from FIG. 2 using the optimized values of the prediction error level γ, the standard deviation σ, and the maximum prediction distance dmax from FIG. 3; and

[0021] FIG. 5 illustrates the application of the method of the present invention to a synthetic data set for quantifying the reliability of predicted sand percentages at seismic locations.

[0022] While the invention will be described in connection with its preferred embodiments, it will be understood that the invention is not limited to these embodiments. On the contrary, it is intended to cover all alternatives, modifications and equivalents that may be included within the scope of the invention, as defined by the appended claims.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0023] The present invention relates to a method for assessing reliability of predicted values of geologic or geophysical parameters from statistical predictive models used in geophysical prospecting. These predictive models could include multiple linear regression and neural networks, but are not restricted to them. The reliability analysis method of the present invention is model free, that is, independent of the predictive model used. The method of the present invention derives reliability measures directly from the topology of the training data points and the predicted values, instead of the residual statistics at the training data points and the predicted values. Therefore, the method is free from assumptions on the form of the residual statistics or the linearity of the underlying predictive model. Although the method of the present invention applies to any statistical prediction model, it will be illustrated using a multiple regression model.

[0024] Consider the case of N training data points (xi, zi), i=1, . . . , N, with a multiple regression function {circumflex over (z)}=f(x), in which {circumflex over (z)} is a predicted value of z for an attribute vector x. A residual ri at the ith training point (xi, zi) is defined by:

r 1 =z i −{circumflex over (z)} i ,   (3)

[0025] for i=1, . . . , N. In the method of the present invention, a linear relationship between x and z is not assumed. Instead, the reliability of a predicted value {circumflex over (z)} is assessed solely based on the N training data points and the residual values observed at the training data points. For a predicted value {circumflex over (z)}, if the attribute vector x is dissimilar from all the training data attribute vectors, xi, i=1, . . . , N, then there may not be enough information available to judge if {circumflex over (z)} is correct or incorrect. This situation is labeled as a prediction {circumflex over (z)} with high unpredictability. On the other hand, reliability of a predicted value {circumflex over (z)} will be high if its attribute vector x is close to some of the training data attribute vectors and the corresponding residuals ri at the neighboring xi are small. This is the case of a prediction {circumflex over (z)} with high reliability. If an attribute vector x is close to some of the training data attribute vectors, but the residuals at these neighboring attribute vectors are large, it is highly probable that predicted value {circumflex over (z)} will be incorrect. This is the case of a prediction {circumflex over (z)} with high unreliability. The above described reliability measures are useful in assessing reliability associated with a predicted value, but these measures are nonetheless qualitative. In the method of the present invention, a mathematical model that quantifies these reliability measures is described.

[0026] To test the validity of a predicted value {circumflex over (z)} relative to the magnitude of prediction error tolerable, as given by a parameter γ, an approach in classical probability theory may start with two hypotheses:

[0027] H 0 : |r|=|z−{circumflex over (z)}|≦γ,   (4a)

[0028] and

H 1 : |r|=|z−{circumflex over (z)}|>γ.   (4b)

[0029] Here, the residual r is defined as in residual equation (3) by r=z−{circumflex over (z)}. Reliability equation (4a) represents the hypothesis H0 that the prediction residual is less than or equal to the prediction error limit γ, while unreliability equation (4b) represents the hypothesis H1 that the prediction residual is greater than γ. Hypothesis H0 is related to a notion of reliability of the predicted value {circumflex over (z)}, while hypothesis H1 is related to a notion of unreliability of the predicted value {circumflex over (z)}. In equations (4a) and (4b), the parameter γ must be appropriately determined, as will be further described later.

[0030] In conventional probability theory, the probability of each hypothesis is calculated based on the residual statistics. In turn, the residual statistics are computed from the residuals observed at the training data points. In this prior art approach, all residuals participate equally in computing residual statistics.

[0031] Conventional residual statistics computation is based on the values of the residuals ri, but ignores their locations xi in attribute space. The method of the present invention calculates reliability measures that are sensitive not only to the residual values but also to the locations of the training data point attribute vectors in attribute space.

[0032] In the method of the present invention, the training point triplets, that is, the training data points and their residuals, (xi, zi, ri), i=1, . . . , N, are considered as N pieces of evidence from which to evaluate the reliability in a predicted value {circumflex over (z)}. The Dempster-Shafer theory of evidence is used to aggregate uncertainties on the predicted value {circumflex over (z)} from each training point triplet (xi, zi, ri), i=1, . . . , N. The Dempster-Shafer theory of evidence is described in Dempster, A. P., “Upper and Lower Probabilities Induced by Multivalued Mapping”, Annals Math. Statistics, Vol. 38, No. 2, p. 325-339, 1987 and Shafer, G., A Mathematical Theory of Evidence, Princeton University Press, 1976.

[0033] In classical probability theory, the sample space of hypotheses consists of all the standard hypotheses: H0 and H1 as given in equations (4a) and (4b). In Dempster-Shafer theory, the sample space of standard hypotheses is enlarged to include all possible subsets of the set of standard hypotheses. This enlarged set of hypotheses is called a frame of discernment θ={Hk}. In our case, one more hypothesis is needed:

Hθ=H0-or- H1.   (4c)

[0034] Here the notation H0 -or- H1 means that H0 and H1 cannot be distinguished. Hypothesis Hθ is related to a notion of unpredictability of the predicted value {circumflex over (z)}. Thus, in our case, the frame of discernment θ={H0, H1, Hθ}. A frame of discernment in Dempster-Shafer theory is a generalization of a sample space of standard hypotheses in conventional probability theory.

[0035] The hypothesis Hθ represents noncommitment to either H0 or H1. In this expanded hypothesis set, the impact of a piece of evidence on each hypothesis is represented by a function m(Hk) called a basic probability assignment. For each piece of evidence represented by a training data triplet (xi, zi, ri), the sum of the basic probability assignments assigned to all three hypotheses, H0, H1, and Hθ, becomes unity. Thus, for each i=1, . . . , N,

m i (H0)+mi(H1)+mi(Hθ)=1   (5)

[0036] This framework provides a means to represent lack of relevant information associated with each piece of evidence by the basic probability assignment mi(Hθ). The smaller the value of mi(Hθ), the more definite the information on H0 or H1.

[0037] The present invention quantifies reliability associated with a predicted value in terms of the values of three reliability measures. These measures are termed reliability, unreliability, and unpredictability, and are associated with the hypotheses given in equations (4a), (4b), and (4c), respectively. Further, equation (5) dictates that the values of the measures add up to one for each i=1, . . . , N. The complement of unpredictability, (1-unpredictability), measures the degree with which the predicted value is supported by the training data points. This concept is not available in the standard confidence interval approach.

[0038] In our case, each basic probability distribution mi(z) is the likelihood of a dependent variable z being evidenced by an observation of a training data triplet (xi, zi, ri). This ith observed training data triplet provides valuable information in that the predicted value {circumflex over (z)}i at xi misses the true value zi by its residual ri. It is assumed that the ith training data triplet (xi, zi, ri) will predict that the probability distribution of the true value z at x will have a mean value of ({circumflex over (z)}+ri). It is also assumed that the prediction power of the ith training data triplet as a piece of evidence concerning z will diminish as the distance between the two attribute vectors x and xi increases. Thus, a normalized distance di, for i=1, . . . , N, between attribute vectors x and xi is defined by: 1$\begin{array}{cc}{d}_1=\sqrt{\sum _{j=1}^J{\left(x_j^n-x_{i,j}^n\right)}^2/J}& \left(6\right)\end{array}$

[0039] where 2$x_j^n=100\left(x_j/\left(x_j^{\max }-x_j^{\min }\right)\right),\mathrm{and}$$x_{i,j}^n=100\left(x_{i,j}/\left(x_j^{\max }-x_j^{\min }\right)\right).$

[0040] Here xj and xi,j are the jth element of the attribute vectors x and xi, respectively, xjmax and xjmin are the maximum and minimum values, respectively, of the jth attribute element for the whole training data set, xjn and xi,jn are normalized attribute elements, and J is the dimension of the attribute vector x. The basic probability assignment on Hθ, mi(Hθ), will increase as di increases.

[0041] It is further assumed that the basic probability distribution of the dependent variable z takes the form of a normal distribution. Based on the notation described above, a basic probability distribution mi(z), as evidenced by the ith training data triplet (xi, zi, ri), is defined as:

m i (z)=θc(di)(({circumflex over (z)}+ri),σ)   (7a)

[0042] for i=1, . . . , N, where 3$\begin{array}{cc}\begin{array}{cc}{\theta }^c\left(d_i\right)=1-\left(d_i/d_{\max }\right)& \mathrm{for}{d}_i<d_{\max }\\ =0& \mathrm{for}{d}_i\ge d_{\max }.\end{array}& \text{(7b)}\end{array}$

[0043] Here θc(di) is the prediction power of the ith training data triplet (xi, zi, ri), dmax is the maximum prediction distance, and (({circumflex over (z)}+ri),σ) represents a normal distribution over z with mean ({circumflex over (z)}+ri) and a standard deviation σ. In equations (7a) and (7b), the parameters σ and dmax are not yet known and must be appropriately determined, as will be described below.

[0044] With the formulations given above, the basic probability assignments of the hypotheses from the evidence given by the ith training data triplet (xi, zi, rj), i=1, . . . , N, are computed as follows: 4$\begin{array}{cc}{m}_i\left(H_0\right)=m_i\left(\&LeftBracketingBar;z-\stackrel{⋒}{z}\&RightBracketingBar;\le \gamma \right)={\int }_{\stackrel{⋒}{z}-\gamma }^{\stackrel{⋒}{z}+\gamma }{m}_i\left(z\right)dz& \text{(8a)}\\ m_i\left(H_1\right)=m_i\left(\&LeftBracketingBar;z-\stackrel{⋒}{z}\&RightBracketingBar;>\gamma \right)={\int }_{-\infty }^{\stackrel{⋒}{z}-\gamma }{m}_i\left(z\right)dz+{\int }_{\stackrel{⋒}{z}+\gamma }^{+\infty }{m}_i\left(z\right)dz& \text{(8b)}\\ m_i\left(H_{\theta }\right)=1-{\theta }^c\left(d_i\right)=1-m_i\left(H_0\right)-m_i\left(H_1\right)& \text{(8c)}\end{array}$

[0045] Here mi(Hθ) represents the amount of belief not committed to either hypothesis H0 or H1. The above three quantities defined in basic probability assignment equations (8a), (8b), and (8c) are termed γ-level reliability, γ-level unreliability, and unpredictability, respectively.

[0046] Dempster-Shafer theory has a formal process, called the rule of combination, to combine individual basic probability assignments induced by independent evidence. If m1 and m2 are two basic probability assignments from two pieces of evidence, and Ai and Bj are hypotheses induced by the two evidences, then a combined basic probability assignment m1⊕m2 (C) for a hypothesis C is computed as: 5$\begin{array}{cc}{m}_1\oplus m_2\left(C\right)=\left(1/K\right)\sum _{A_i\bigcap B_j=C}{m}_1\left(A_i\right)m_2\left(B_j\right)& \left(9\right)\end{array}$

[0047] where K is a normalization factor defined by: 6$K=1-\sum _{A_i\bigcap B_j=\phi }{m}_1\left(A_i\right)m_2\left(B_j\right)$

[0048] and Φ is the null set. The combination is designated by the orthogonal sum operation ⊕.

[0049] For the specific case of hypotheses H0, H1, and Hθ, we have

H0=H0∩H0=H0∩Hθ,

H1=H1∩H1=H1∩Hθ,

Hθ=Hθ∩Hθ, and

Φ=H0∩H1.

[0050] Then, the rule of combination in equation (9) becomes:

m 1 ⊕m 2 (H0)=[m1(H0)m2(H0)+m1(H0)m2(Hθ)+m2(H0)m1(Hθ)]/K,

m 1 ⊕m 2 (H1)=[m1(H1)m2(H1)+m1(H1)m2(Hθ)+m2(H1)m1(Hθ)]/K,

[0051] and

m1⊕m2(Hθ)=m1(Hθ)m2(Hθ)/K

[0052] with normalization factor K defined by:

K= 1−m1(H0)m2(H1)−m1(H1)m2(H0)

[0053] The commutative property of multiplication ensures that this rule yields the same value regardless of the order in which basic probability assignments are combined. In the present invention, the total belief on a hypothesis is computed by combining the N basic probability assignments from the N training data triplets. In the present case, this gives:

Bel(H0)=m1⊕m2⊕ . . . ⊕mN(H0),   (10a)

Bel(H1)=m1⊕m2⊕ . . . ⊕mN(H1), and   (10b)

Bel(Hθ)=m1⊕m2β . . . ⊕mN(Hθ).   (10c)

[0054] Belief Bel (Hk), for k=0, 1, and θ, may be interpreted as the total support that is attributed to the hypothesis Hk by the evidence from all N training data triplets.

[0055] Belief equations (10a)-(10c) define the three reliability measures of the present invention. Reliability is defined by Bel (H0) in equation (10a), unreliability is defined by Bel (H1) in equation (10b), and unpredictability is defined by Bel (Hθ) in equation (10c). By their construction, these three reliability measures have the following characteristics for a prediction point x. Reliability, Bel (H0), will be large if the attribute vector x is near the vicinity of some of the training data point attribute vectors xi and the residuals ri at these nearby training data points are small. Unreliability, Bel (H1), will be large if the attribute vector x is near the vicinity of some of the training data point attribute vectors xi and the residuals ri at the nearby training data points are large. Unpredictability, Bel (Hθ), will be large if the attribute vector x has only a few or no training data point attribute vectors xi in its vicinity.

[0056] The lower part of FIG. 1 shows 10 training data points 100 and the linear prediction line 102, in which the horizontal axis 104 represents the attribute vector x of one dimension and the vertical axis 106 represents the predicted variable z. The upper part of FIG. 1 displays the three reliability measures computed for the prediction line. These measures are reliability 108, unreliability 110 and unpredictability 112. One can observe the three characteristics described above in this plot. A comparison of three reliability measures identifies regions of high reliability 114 or low reliability 116 on the predicted values. It also indicates regions 118 where poorly constrained extrapolation or interpolation might occur, which is signified by a high unpredictability value. In this example, the regions of high reliability 114 or high unpredictability 118 are relatively easy to identify visually without the help of reliability measures, because the attribute vector has only one element or dimension. However, when the dimension of an attribute vector becomes greater than two, this type of visual analysis becomes very difficult.

[0057] In the present invention, three parameters affect the quantification of the three reliability measures. With reference to equations (7) and (8), the choices of parameter values for prediction error level γ, standard deviation σ, and maximum prediction distance dmax will result in different values for the reliability measures. The following discusses effects of these parameters and establishes a preferred procedure for determining an optimal set of values for them. However, other methods for determining the parameter values are within the scope of the present invention.

[0058] First, an appropriate value for prediction error level γ is determined. The parameter γ defines the amount of prediction error one is willing to allow in the definition of γ-level reliability or γ-level unreliability in basic probability assignment equations (8a) and (8b). If γ is unreasonably large, then all γ-level reliabilities will tend to have high values. Conversely, if γ is too small, then all γ-level reliabilities will tend to have low values. Therefore, improper choice of γ will hamper the resolution of reliability measures, even though the relative magnitudes will still have meaning.

[0059] A sample variance σt2 of the residuals ri from the set of N training data triplets is defined as: 7$\begin{array}{cc}{\sigma }_t^2=\sum _{i=1}^N{r}_i^2/\left(N-1\right)& \left(11\right)\end{array}$

[0060] If one would like to assess the prediction reliability as compared to one sample standard deviation of the residuals of the training data triplets, the preferred embodiment would be to set the prediction error level γ to

γ=σ5   (12)

[0061] Alternatively, if a limit on acceptable error tolerance exists, such as 5 porosity units for predicted porosities, then an alternative embodiment would be to set the prediction error level γ to the error tolerance limit, such as 5 in our example.

[0062] Next, a range of acceptable values for the standard deviation σ is determined. In basic probability distribution equations (7a) and (7b), σ2 is the variance of the quantity z−({circumflex over (z)}+ri), which can be expanded into:

z −({circumflex over (z)}+ri)=(z−{circumflex over (z)})−(zi−{circumflex over (z)}i)   (13)

[0063] Without assuming certain statistical properties for the residuals, such as constant variance or correlation structure of residuals, the variance of the quantity in equation (13) cannot be estimated. Therefore, a constant variance of residuals is assumed to find an approximate range of values for σ. If the statistical prediction model {circumflex over (z)}={circumflex over (f)}(x) describes the underlying relationship between attribute vector x and dependent variable z fairly well everywhere, then variance Var(z−{circumflex over (z)}) will be the same as variance Var(zi−{circumflex over (z)}i). If it can further be assumed that the residuals (z−{circumflex over (z)}) and (zi−{circumflex over (z)}i) are uncorrelated, then 8$\begin{array}{cc}\begin{array}{c}{\sigma }^2=\mathrm{Var}\left(z-\left(\stackrel{⋒}{z}+r_i\right)\right)=\mathrm{Var}\left(\left(z-\stackrel{⋒}{z}\right)-\left(z_i-{\stackrel{⋒}{z}}_i\right)\right)\\ =\mathrm{Var}\left(z-\stackrel{⋒}{z}\right)+\mathrm{Var}\left(z_i-{\stackrel{⋒}{z}}_i\right)\approx 2{\sigma }_t^2\end{array}& \left(14\right)\end{array}$

[0064] On the other hand, if it is assumed that the residuals (z−{circumflex over (z)}) and (zi−{circumflex over (z)}i) are positively correlated, that is, the values of the residuals are not random, then 9$\begin{array}{cc}{\sigma }^2<2{\sigma }_t^2& \left(15\right)\end{array}$

[0065] The assumption in a conventional statistical analysis is that the residuals are completely uncorrelated, implying that one could set 10${\sigma }^2=2{\sigma }_t^2.$

[0066] In general, however, σ should be set to a value less than {square root}{square root over (2)}σt if correlations exist among the residuals. This determines a range of values for the standard deviation σ.

[0067] If an unreasonably high value is chosen for the standard deviation σ, then the normal distribution in basic probability distribution equations (7a) and (7b) flattens, resulting in loss of distinction between reliability and unreliability measures. Conversely, if σ is too small, then the contrast between reliability and unreliability measures will be exaggerated. In either case, results are less desirable, although still meaningful.

[0068] Next, a range of acceptable values for the maximum prediction distance dmax is determined. In basic probability distribution equations (7a) and (7b), the prediction power, θc(di), of the ith piece of evidence on the true value z is a function of dmax. The prediction power of θc(di) becomes 1 when the distance di is zero and it decreases linearly to zero when di=dmax. Therefore, a large value of dmax implies that prediction power of evidences reaches far in attribute space. If one believes that the residuals (z−{circumflex over (z)}) and (zi−{circumflex over (z)}i) are correlated at large distances in attribute space, then dmax should be assigned a large value, as a percentage of the maximum of the normalized attribute distances di defined in equation (6) above. For example, a value for dmax of 30 implies that the correlation between residuals is assumed to be zero at and beyond the 30 percent of the maximum normalized attribute distance. In the preferred embodiment, dmax is set to a value less than 100. This determines a range of values for the maximum prediction distance dmax.

[0069] Finally, optimal values for σ and dmax are determined within the ranges of values for σ and dmax determined above. So far, it has been noted that the two parameters σ and dmax should be a function of the correlation between the two residuals r and ri, which may not be constant and is difficult to estimate with a limited number of training data points. For a given value of γ, the parameters σ and dmax will be determined so as to maximize the average correlation between probabilities based on an observed residual at a training data point and γ-level reliabilities based on the rest of the observed residuals. In the present invention, random variables zit for i=1, . . . , N are used. These random variables zit are defined to have a normal distribution with a mean zi and a standard deviation σ. This assumption allows one to compute the probability that the random variable Zit is between {circumflex over (z)}i+γ and {circumflex over (z)}i−γ, given the observation of zi, for i=1, . . . , N. This is expressed as

P(|zit−{circumflex over (z)}i|<γ|zi)   (16)

[0070] The γ-level reliability based on the restricted set of training data triplets (xj, zj, rj) for j=1, . . . , N and j≠i, can be computed as the degree of belief that the residual between the random variable Zit and the predicted value {circumflex over (z)}i is less than γ. This modified belief is then given by

[0071] Beli′(H0)=mi⊕ . . . ⊕mi-1⊕mi+1⊕ . . . ⊕mN(H0)   (17)

[0072] Note that the probability in equation (16) is based on the single observed parameter value zi, while the belief in equation (17) is based on the remaining observed parameter values Zj for j=1, . . . , N and j≠i. Further, the belief in equation (17) is a function of σ and dmax, while the probability in equation (16) is independent of σ and dmax.

[0073] Now, for fixed values of σ and dmax, the average of the correlation coefficients between the probabilities from equation (16) and the γ-level reliabilities from equation (17) for i=1, . . . , N, is computed as: 11$\begin{array}{cc}\rho \left(\sigma ,d_{\max }\right)=\frac{1}{N}\sum _{i=1}^{N}P\left(\&LeftBracketingBar;z_i^t-{\stackrel{⋒}{z}}_i\&RightBracketingBar;<\gamma |z_i\right){\mathrm{Bel}}_i^{\prime }\left(H_0\right).& \left(18\right)\end{array}$

[0074] Here, the average correlation ρ(σ, dmax) is a function of the two parameters σ and dmax.

[0075] By knowing an approximate range of values for σ and dmax, 12${\sigma }^2<2{\sigma }_t^2$

[0076] and dmax <100, the average correlation ρ in correlation equation (18) can be computed for various values of σ and dmax with the training data points. Because the computation of the average correlation is very fast, a parameter set of σ and dmax which maximizes the average correlation is easily found by an exhaustive search. With an appropriate value for the parameter γ and optimum values for the parameters σ and dmax, equations (7) through (10) may be used to calculate the desired reliability measures.

[0077] FIG. 2 is a flowchart which illustrates a preferred embodiment of the method of the present invention for determining reliability measures for prediction values for a specified geophysical or geological parameter from a statistical prediction model. First, at step 200, N locations for training data are selected. Preferably, the training point locations are at or near N well locations. Second, at step 202, N attribute vectors xi for i=1, . . . , N are selected at the N training data locations from step 200. Preferably, these attribute vectors comprise one or more seismic attributes obtained from seismic data traces collected at or near N wells. Next, at step 204, N associated observed values zi for i=1, . . . , N are selected at the N training data locations from step 200 and corresponding to the N attribute vectors from step 202. Preferably, these associated observed values are taken from well data at the N wells. Next, at step 206, these attribute vectors from step 202 and associated observed values from step 204 are combined to form N training data points (xi, zi) for i=1, . . . , N. Next, at step 208, a statistical prediction model is constructed from the N training data points from step 206. Preferably, this prediction model is constructed by multiple regression on the N training data points from step 206. However, the method of the present invention is not limited to any particular statistical prediction model. Next, at step 210, predicted values {circumflex over (z)}i for i=1, . . . , N, are determined for each of the N attribute vectors xi from step 202, using the statistical prediction model from step 208. If the prediction model is defined by a function f, then {circumflex over (z)}i=f(xi) for i=1, . . . , N. Next, at step 212, residuals ri for i=1, . . . , N are determined for each of the N training data points from step 206. Preferably, each of these residuals ri is defined by the difference between an associated observed value zi from step 204 for the specified parameter for the training data point and the corresponding predicted value {circumflex over (z)}i from step 210 of the specified parameter for the training data point from the statistical prediction model. Thus ri=zi−{circumflex over (z)}i for i=1, . . . , N. Next, at step 214, the training data points (xi, zi) for i=1, . . . , N from step 206 and the residuals ri for i=1, . . . , N from step 212 are combined to form N training data triplets (xi, zi, ri) for i =1, ..., N.

[0078] Next, at step 216, a designated location is selected. Preferably, this designated location is in an area in which the same one or more seismic attributes collected at the N training data locations in step 202 can be obtained from seismic data traces. Next, at step 218, an attribute vector x is determined at the designated location from step 216. Preferably, this attribute vector comprises the one or more seismic attributes obtained from seismic data traces collected at or near the designated location. Next, at step 220, a predicted value {circumflex over (z)} of the specified parameter is determined at the designated location from step 216. Preferably, the predicted value {circumflex over (z)} of the specified parameter is obtained by applying the statistical prediction model from step 208 to the attribute vector x from step 218. If the prediction model is defined by a function f, then {circumflex over (z)}=f(x). Next, at step 222, N basic probability distributions mi(z) for i=1, . . . , N are defined from the N training data triplets (xi, zi, ri) from step 214 and the predicted value {circumflex over (z)} from step 220. The preferred construction of the N basic probability distributions as in equations (7a) and (7b) will be covered in more detail in the discussion of FIG. 4 below.

[0079] Lastly, the reliability measures are calculated. First, at step 224, a hypothesis H0 that the predicted value {circumflex over (z)} from step 220 is reliable is defined. Preferably, the hypothesis H0 is defined by equation (4a), above. Similarly, at steps 226 and 228, respectively, hypotheses H1 and Hθ that the predicted value {circumflex over (z)} from step 220 is unreliable and unpredictable, respectively, are defined. Preferably, the hypotheses H1 and Hθ are defined by equations (4b) and (4c), respectively, above. Next, at step 230, respectively, basic probability assignments mi(H0) for i=1, . . . , N are determined for the predicted value {circumflex over (z)} from step 220 and for the hypotheses H0 from step 224, using the N training data triplets (xi, zi, ri) from step 214. Preferably, basic probability assignments mi(H0) are defined by equation (8a), above. Similarly, at steps 232 and 234, respectively, basic probability assignments mi(H1) and mi(Hθ) for i=1, . . . , N are determined for the predicted value {circumflex over (z)} from step 220 and for the two hypotheses H1 and Hθfrom steps 226 and 228, respectively, using the N training data triplets (xi, zi, ri) from step 214. Preferably, basic probability assignments mi(H1) and mi(Hθ) are defined by equations (8b) and (8c), respectively, above. Finally, at step 236, a reliability measure Bel (H0), an unreliability measure Bel (H1), and an unpredictability measure Bel (Hθ) for the predicted value {circumflex over (z)} from step 220 are determined. The reliability measure Bel (H0), unreliability measure Bel (H1), and unpredictability measure Bel (Hθ) are determined as combinations of the N basic probability assignments mi(Hk) for i=1, . . . , N and k=0, 1, and θ, from steps 230, 232, and 234. Preferably, reliability measure Bel (H0), unreliability measure Bel (H1), and unpredictability measure Bel (Hθ) are calculated using equations (10a), (10b), and (10c), respectively.

[0080] FIG. 3 is a flowchart that illustrates a preferred embodiment of the method of the present invention for fine-tuning the prediction error level γ, standard deviation σ, and maximum prediction distance dmax parameters. First, at step 300, a sample variance σt2 is calculated for the N residuals ri from step 212 of FIG. 2. Preferably, the sample variance σt2 is defined as in equation (11) above. Second, at step 302, a prediction error level γ is set. Preferably, the prediction error level γ is set equal to the sample standard deviation σt from step 300. Alternatively, if a limit on error tolerance exists, then the value of γ is set equal to the error tolerance limit.

[0081] Next, at step 304, a modified γ-level reliability Beli′(H0) is calculated using the N−1 training data points (xj, zj, rj) for j=1, . . . , N and j≠i, from step 206 of FIG. 2 and the prediction error level γ from step 302. Preferably, the modified γ-level reliability Beli′(H0) is calculated using equation (17), above. Next, at step 306, random variables Zit are defined for i=1, . . . , N. Preferably, each random variable zit is defined having a normal distribution with a mean equal to the associated observed value zi from step 204. The normal distribution of the random variables zit will each have a standard deviation σ, which will be determined below. Next, at step 308, a residual rit is determined for i=1, . . . , N. Preferably, each residual rit is defined by the difference between the random variable zit from step 304 and the predicted value {circumflex over (z)}i from step 210 of FIG. 2. Next, at step 310, probabilities P (|zit−{circumflex over (z)}i|<γ|zi) are calculated for the residuals rit from step 308, for i=1, . . . , N and the prediction error level γ from step 302. Next, at step 312, an average correlation ρ between the γ-level reliabilities from step 304 and the probabilities from step 310 are calculated. Preferably, the average correlation ρ is defined as in equation (18) above.

[0082] Next, at step 314, the ranges of the variables σ and dmax are limited. Preferably, if the residuals are correlated, the limit on σ is given by variance inequality (15), using the sample variance σt2 from step 300. Alternatively, if the residuals are uncorrelated, σ is approximated by variance approximation (14), again using the sample variance σt2 from step 300. Preferably, dmax is limited to less than 100. Next, at step 316, the values of σ and dmax which maximize the average correlation ρ from step 312 is determined by an optimization algorithm. Preferably, the optimization algorithm is an exhaustive search over the range of limited values for σ and dmax from step 314.

[0083] FIG. 4 is a flowchart which illustrates a preferred embodiment of the method of the present invention for calculating the reliability, unreliability, and unpredictability measures in step 236 of FIG. 2, using the selected values of the prediction error level γ from step 302 of FIG. 3 and the standard deviation σ and maximum prediction distance dmax from step 316 of FIG. 3. First, at step 400, a normalized distance di for i=1, . . . , N is defined for each training data point from step 206 of FIG. 2. Preferably, the normalized distance di is defined as in normalized distance equation (6) above, using values for the attribute vectors xi and x from steps 202 and 218, respectively, of FIG. 2. Next, at step 402, prediction powers θc(di) for i=1, . . . , N are defined for each training data point. Preferably, the prediction powers θc(di) are defined as in equation (7b) above, using values for the normalized distance di from step 400 and the maximum prediction distance dmax from step 316 of FIG. 3. Next, at step 404, normal distributions (({circumflex over (z)}+ri),σ) for i=1, . . . , N are defined, based on each training data point. Preferably, these normal distributions (({circumflex over (z)}+ri),σ) are defined with a mean of {circumflex over (z)}+ri, using values for the residuals ri and the predicted value {circumflex over (z)} from steps 212 and 220, respectively, of FIG. 2, and the standard deviation σ from step 316 of FIG. 3. Then, at step 406, basic probability distributions mi(z) for i=1, . . . , N are defined for each training data point. These are the basic probability distributions referred to in step 222 of FIG. 2. Preferably, the basic probability distributions mi(z) are defined by multiplying the prediction powers θc(di) from step 402 with the normal distributions (({circumflex over (z)}+ri),σ) from step 404, as in probability equation (7a) above. Next, at step 408, basic probability assignments mi(Hk) are calculated from basic probability distributions mi(z) from step 406 for each of the hypotheses H0, H1, and Hθ from steps 230, 232, and 234, respectively, of FIG. 2. Preferably, the basic probability assignments mi(H0), mi(H1), and mi(Hθ) are defined as in equations (8a), (8b), and (8c), respectively, above, using values for the predicted value {circumflex over (z)} from step 220 of FIG. 2 and for the prediction error level γ from step 302 of FIG. 3. These basic probability assignments from step 408 define γ-level reliability, γ-level unreliability, and unpredictability, respectively, as evidenced by a single training data triplet (xi, zi, ri). Finally, at step 410, total belief is calculated for each of the hypotheses H0, H1, and Hθ, as a combination of the basic probability assignments from step 408. Preferably, the total beliefs Bel (H0), Bel (H1), and Bel (Hθ) are calculated by using the Dempster-Shafer theory rule of combination, as in equations (10a), (10b), and (10c) above. These total beliefs from step 410 define reliability, unreliability, and unpredictability measures, respectively, for a given prediction error level γ, as evidenced by the N training data triplets (xi, zi, ri) for i=1, . . . , N.

EXAMPLES

[0084] The method of the present invention is further described with the following examples. With a 3-D seismic survey, seismic data and the accompanying measured attributes are available at close intervals in the space domain. Even though well-log data in the same area may provide important geophysical or geological information, typically wells are relatively sparsely located. Geoscientists often use linear regression to establish relationships between seismic attributes and geophysical or geological properties observed at wells. This regression model, along with seismic attributes, is then used to predict the unknown geophysical or geological properties away from the wells. The first example demonstrates the value of the present invention in this type of seismic attribute calibration. Also shown in this example is a demonstration of how the parameter optimizing methods are employed.

[0085] A 100×100 synthetic seismic attribute data set, which simulates a fluvial environment, is provided, along with sand percentages at 30 well locations. In this example, four seismic attributes, named χ1, χ2, χ3, and χ4, were used to predict the sand percent:

[0086] x=(χ1, χ2, χ3, χ4)

[0087] z=sand-percent

[0088] (xi, zi), i=1, . . . , 30 available at well locations.

[0089] The linear regression model, calculated using standard least squares techniques applied to the 30 training data points, was:

{circumflex over (z)}=â T x+{circumflex over (b)}   (19)

[0090] where

â T =(−1.512, −1.4242, −0.7811, 1.2334)

{circumflex over (b)}T=171.77

[0091] σt=7.51 sand percentage.

[0092] The latter was calculated from residuals, using equation (11). The prediction model in model equation (19) was used to predict sand percentages at all seismic locations. The actual sand percent 500, predicted sand percent 502, and prediction error 504 are plotted in FIG. 5. In a real situation, only the predicted sand percentage map would be available because the actual sand percentages would not be known. The task of the present invention is to evaluate the reliability of the predicted sand percentages.

[0093] By following the procedure for optimal parameter selection described in conjunction with FIG. 3, γ=8 sand percentage was chosen, which is close to the sample standard deviation of the residuals of the training data points. In this case, the γ-level reliability computes the degree of belief that the difference between a true sand percentage and a predicted sand percentage is less than 8 sand percent. For a fixed value of γ=8, average correlations ρ between the probabilities and the γ-level reliabilities for various values of σ and dmax are calculated from equation (18) and listed in Table 1.

	TABLE 1
	σ
	5	8	10
	dmax	3	0.420	0.404	0.392
		6	0.421	0.406	0.394
		9	0.432	0.408	0.390
		12	0.423	0.390	0.360
		15	0.354	0.331	0.284

[0094] In Table 1, the parameter set σ=5 and dmax =9 maximizes the average correlation.

[0095] Next, three reliability measures were computed for the rest of the seismic attribute data. The three lower maps in FIG. 5 display the degrees of reliability 506, unreliability 508, and unpredictability 510 associated with the predicted sand percentage map 502 in the middle of the upper row. In this synthetic data case, one can visually inspect how well the method determined the reliability of the predicted sand percentages by comparing the error map 504 with the three reliability maps 506, 508, and 510. It is apparent that there tend to be large errors in areas, such as 512, where unreliability measures are higher than reliability measures. Conversely, error tends to be small in areas, such as 514, where the reliability is high and the unreliability is low. Areas with high unpredictability, such as 516, indicate places where lack of training data points prohibits reliable predictions. These high unpredictability regions correspond to an area with no nearby training data points in the attribute space. This example demonstrates that the technique of the present invention can identify areas of reliable and unreliable predictions. It also identifies areas of unpredictability, which correspond to the area with no attribute training data.

[0096] In this synthetic data case, true sand percentages are known at every attribute location. To quantify how well the belief model worked, we computed the average correlation from correlation equation (18) for the whole data set. The results are shown in Table 2.

	TABLE 2
	σ
	5	8	10
	dmax	6	0.60	0.56	0.52
		9	0.59	0.55	0.51

[0097] Table 2 shows that the average correlation between the probabilities from equation (16) and the reliabilities from equation (17) is 0.59 for γ=8, σ=5, and dmax=9. Because the reliability measure is most meaningful when the unpredictability measure is not high, the correlation was also calculated for only the data points at which unpredictability is less than 0.9.

	TABLE 3
	σ
	5	8	10
	dmax	6	0.84	0.77	0.73
		9	0.83	0.77	0.73

[0098] Table 3 shows that the correlation is then 0.83 for γ=8, σ=5, and dmax=9. This relatively high correlation coefficient indicates that the method of the present invention is able to determine the reliability of predictions with high accuracy.

[0099] The technique of the present invention does not assume a linear relationship between the independent variables, the attributes, and the dependent variable. The present invention is model independent. Therefore, the present invention can be applied to evaluate reliabilities of nonlinear prediction results. A certain type of artificial neural network, such as a back-propagation, feedforward neural network (BPN), can be considered as a nonlinear prediction model. Even though BPN is a powerful tool for establishing a nonlinear relationship between input (attribute) and output (dependent) variables, the particular nonlinear relationship created is not apparent to the user. In a neural network application, only the elemental functional form (sigmoid) is known, not the relationship created with multiple layers of its elemental (sigmoidal) functions. Thus, in BPN applications, an additional test data set must be set aside to verify the model's performance after training is completed. Even then, the reliability of the model's prediction is unknown. Thus, the technique of the present invention will be valuable in evaluating predictions made from a BPN.

[0100] The three reliability measures derived in the present invention will provide more meaningful reliability measures than standard confidence intervals when:

[0101] (1) the training data points are clustered into groups separated by data gaps or have insufficient coverage of the attribute space,

[0102] (2) the residuals are correlated,

[0103] (3) the variance of the residuals is not constant,

[0104] (4) the dimension of the attribute space is greater than two, which makes it difficult to visualize gaps in training data points or a certain trend in residual variances, and

[0105] (5) the degree of freedom in a prediction model is greater than 8, so that the average correlation in equation (18) between the probabilities in equation

[0106] (17) and the modified γ-level reliabilities in equation (17) is meaningful.

[0107] The foregoing description has been directed to particular embodiments of the invention for the purpose of illustrating the invention. It will be apparent to persons skilled in the art, however, that many modifications and variations to the embodiments described herein are possible. For example, in any of the foregoing embodiments, any statistical prediction model can be used. All such modifications and variations are intended to be within the scope of the present invention, as defined by the appended claims.

Claims

1. A method for assessing reliability of a predicted value of a specified geological or geophysical parameter at a designated location, the predicted value obtained from a prediction model constructed from N training data attribute vectors and N associated observed values of the specified parameter; each training data attribute vector including one or more seismic attributes obtained from seismic data traces located at or near a well, each associated observed value of the specified parameter obtained from well log or core data from the well, the method comprising the steps of: determining a predicted value of the specified parameter for each of the N training data attribute vectors, from the training data attribute vectors and the prediction model; determining a residual for each of the N training data attribute vectors, the residual being the difference between the associated observed value of the specified parameter for the training data attribute vector and the predicted value of the specified parameter for the training data attribute vector ; determining an attribute vector for the designated location; determining the predicted value of the specified parameter at the designated location, from the attribute vector for the designated location and the prediction model; determining N basic probability distributions from the N training data attribute vectors, the N associated observed values, the N residuals, and the predicted value of the specified parameter at the designated location; determining N basic probability assignments for each of three hypotheses that the predicted value of the specified parameter at the designated location is reliable, unreliable, and unpredictable, respectably, from the N basic probability distributions; and determining a reliability value, an unreliability value, and an unpredictability value for the predicted value of the specified parameter at the designated location as combinations of the N basic probability assignments for each of the hypotheses that the predicted value of the specified parameter at the designated location is reliable, unreliable, and unpredictable, respectively.

2. The method of claim 1, wherein the basic probability distributions mi(z) over the specified parameter z are defined for i=1, . . . , N, by

3. The method of claim 2, wherein the normalized distance di is defined for i=1, . . . N, by

4. The method of claim 3, wherein the normalized jth elements of the attribute vector x, xjn and xijn, are defined by

5. The method of claim 1, in which the basic probability assignments mi(Hk) are defined for i=1, . . . , N and k=0, 1, and θ, by

6. The method of claim 5, wherein the reliability value Bel (H0), the unreliability value Bel (H1), and the unpredictability value Bel (Hθ) are calculated by

7. The method of claim 6, wherein

8. The method of claim 5, further comprising the steps of: calculating a sample standard deviation σt for the N residuals; and setting the prediction error level γ equal to the sample standard deviation σt.

9. The method of claim 5, further comprising the steps of: selecting an error tolerance limit; and setting the prediction error level γ equal to the error tolerance limit.

10. The method of claim 2, further comprising the steps of: defining a random variable zit having a normal distribution with a mean zi and the standard deviation σ, for i=1, . . . , N; calculating a modified γ-level reliability Beli′(H0), for i=1, . . . N , based on a restricted set of the training data triplets (xj, zj, rj) with j=1, . . . , N and j≠i; calculating a probability P (|zit−{circumflex over (z)}i|<γ|zi), for i=1, . . . , N; calculating an average correlation between the modified γ-level reliabilities Beli′(H0) and the probabilities P (|zit−{circumflex over (z)}i|<γ|zi); and determining the values of the standard deviation C and the maximum prediction distance dmax, which maximize the average correlation by an optimization algorithm.

11. The method of claim 10, wherein the modified γ-level reliability Beli′(H0) is defined by

12. The method of claim 10, wherein the average correlation ρ is defined by

13. The method of claim 10, wherein the optimization algorithm is an exhaustive search.