This invention relates generally to the exploration and production of fossil fuels from subterranean reservoirs, and more particularly to improved methods for detecting the presence of elemental concentrations using gamma ray logging.
Over the past years, those involved in the exploration and production of fossil fuels have developed complex methods and tools for evaluating the presence of underground resources. The concentrations of radioactive isotopes of elements such as potassium, uranium and thorium in subsurface earth formations provide valuable geophysical and petrophysical information. The determination of the concentrations of these isotopes is made with radioactive well logging techniques.
In many cases, a logging tool is used to measure the amount of naturally occurring radiation from the formation. Shales often emit more gamma rays than other sedimentary rocks because shales include radioactive potassium, uranium and thorium. Using a multichannel detector, the naturally occurring radiation can be evaluated at a number of different energies and then compared against known spectra corresponding to constituent components that are expected in the formation.
A conventional approach for determining the concentrations from gamma ray spectra is based on minimizing the square of residuals using the weighted least square method expressed by the following equation:
F=½Σi=1NW(Ax−p)2 (1)
where “p” is a vector representing the measured spectrum, “A” is a matrix of the standards for each element, “W” is a weight matrix, N is number of channels in the spectra, with each channel corresponding to a range of energy, and x is a matrix that represents the elemental concentrations. To minimize equation 1, the derivative of F is taken with respect to x. Setting the result equal to zero, x can be estimated as follows:
x=(ATWA)−1ATWp (2)
This equation can be referred to as a direct solver and represents the conventional method for determining the measured elemental concentrations. Importantly, however, the minimization of the cost function (F) without any conditions can give negative concentrations with sensitive spectra. The presence of these non-physical results affects the rest of the analysis. The negative results of any of the elemental concentrations can lead to the wrong value for other elements, so that the conventional approach of limiting the negative concentrations to zero does not solve the problem.
Accordingly, there remains a need for an improved method for determining elemental concentrations using spectral natural gamma ray logging. It is to this and other deficiencies in the prior art that the present invention is directed.
In an embodiment, the present invention includes an improved method for determining elemental concentrations based on measurements from spectral natural gamma ray logging. In an embodiment, the method includes the step of measuring radiation levels in the underground formation with a detector that includes (N) channels. The measured spectral radiation from the underground formation is stored in a measured matrix (p) that includes the measured radiation level at each channel (N). The method continues with the step of providing a standards matrix (A) that includes the established radiation levels for each of the plurality of elements at each energy level corresponding to each channel (N). The process continues with the step of calculating a proportion matrix (x) that provides the concentrations of each of the plurality of elements in the underground formation. The step of calculating the proportion matrix (x) includes using a direct solver equation to determine initial values for the proportion matrix (x), applying the initial values of the proportion matrix (x) to the standards matrix (A) to create an initial proportion model (Ax). Next, the difference between the initial proportion model (Ax) and the measured matrix (p) is determined. A weight matrix (W) is provided and then applied to the difference between the proportion model (Ax) and the measured matrix (p) to normalize the established radiation levels for each of the plurality of elements. The weighted difference between the proportion model (Ax) and the measured matrix (p) across all (N) channels provides the basis for a cost function (F). The cost function (F) is then minimized over a series of iterations to determine a best solution for the proportion matrix (x). During each iteration, the derivatives of the cost function (F) are calculated with respect to (x) to obtain the gradient (dF/dx). A gradient factor is then determined by multiplying the gradient (dF/dx) by a step size factor (α). The values of the proportion matrix (x) are updated by subtracting the gradient factor and iteratively repeating the determination of the proportion matrix (x). The process further includes the step of displaying the calculated values for the proportion matrix (x) on a computer.
Embodiments of the present invention include an improved method for determining elemental concentrations using spectral natural gamma ray logging. Referring to
In an embodiment, the sensor 104 is configured and positioned to detect the emission of naturally occurring gamma ray radiation from the constituent components within the formation 112. The sensor 104 is configured to output signals to the computer 110 representative of the measured radiation. The logging system 100 may optionally include an emitter that irradiates the formation 112 to produce the release of characteristic radiation from the formation.
In an embodiment, the sensor 104 includes about 256 channels that are each configured to measure quantity (counts) of radiation at different energies on the spectrum.
The determination of the proportions of potassium, uranium and thorium in the formation is made using an iterative method 200 shown in
Generally, the method 200 provides an inversion algorithm that seeks to minimize the cost function F (equation 1). Giving an initial guess to the elemental concentrations, the gradient of the cost function is used to update the solution at each iteration. If any of the elemental concentrations are negative, the corresponding derivative is set to zero and the process continues. The iteration is terminated if the difference between the model and measurement becomes suitably small (e.g., 1e-6) or if a suitably large number (e.g., 1000) of iterations have been taken.
Thus, the method 200 begins with a measured matrix (p) that includes the output from the sensor 104 and represents the measured radiation spectrum, a standards matrix (A) that represents that previously established model or signature spectra for each of the elements under evaluation and a weight matrix (W) that is used to normalize the radiation levels within the various spectra. In an embodiment, the measured matrix (p) constitutes a single row matrix with 256 columns that correspond to each channel of the sensor 104, the standards matrix (A) includes three rows (each corresponding to a unique element) with 256 columns that correspond to the established spectra across the 256 channels, and weight matrix (W) is a diagonal matrix that is applied to the difference between the proportion model (Ax) and the measured matrix (p) to normalize the established radiation levels for each of the plurality of elements.
At step 202, the initial values (x0) for a proportion matrix (x) are calculated using the direct solver (equation 2). In the embodiment, proportion matrix (x) is a single row, three-column matrix with each entry corresponding to the proportion of a different one of the three elements under evaluation. The iterative process begins at step 204 with the first iteration (g) with the initial values for the proportion matrix (x) defined as (xold).
At step 206, the cost function (F) (equation 1) is evaluated with the initial values (xold) using a weighted least squares method for each channel (N). Notably, the difference between the initial proportion model (Ax) and the measured matrix (p) is determined. The weight matrix (W) is provided and then applied to the difference between the proportion model (Ax) and the measured matrix (p). The weighted difference between the proportion model (Ax) and the measured matrix (p) across all (N) channels provides the basis for the cost function (F).
At step 208, the derivatives of the cost function (F) are calculated with respect to (x) to obtain the gradient (dF/dx). A provisional new solution for the proportion matrix (xnew) is then calculated by subtracting a gradient factor from the current solution for the proportion matrix (xold). The gradient factor is defined as the product of the gradient (dF/dx) and a step size factor (α). The step size factor (α) is in an embodiment small so that the incremental change between iterations is well controlled. In an embodiment, the step size factor (α) is set at 0.01. The value of the step size factor (α) can be adjusted to change the rate of convergence around a solution.
Next, the method 200 moves to a decision step 212 that queries whether the provisional solution (xnew) includes a negative entry. If an entry (i) within (xnew) is negative, the derivative of (dF/dx) for that element is set to 0, and the value for that entry is returned to the former value (i.e., (xnew(i)=xold(i)). At step 216, (xold) is then set for the subsequent iteration as equal to the determined value of (xnew). If, on the other hand, no entries (i) within (xnew) are negative, the method moves directly to step 216 without the intervening step 214 and (xold) is updated for the subsequent iteration as the value of (xnew).
Next, the method 200 moves to two decision steps 218, 220. At decision step 218, the method queries whether the value of the cost function (F) determined at step 206 during the current iteration is sufficiently small. In an embodiment, the decision step 218 queries whether the cost function returned a result less than 1×10−6. If so, the method 200 moves to step 222 and the results of the proportion matrix (x) are displayed and the method 200 ends. If not, the method 200 progresses to step 220, which queries whether a predefined number of iterations have occurred. In an embodiment, the maximum number of iterations is set at 1,000. If less than 1,000 iterations have occurred, the method moves to step 224 and the iteration count (g) is incremented by 1 before returning to step 206. If the predefined number of iterations has occurred (e.g., g=1000), the method moves to step 222 and the results of the proportion matrix (x) are displayed and the method 200 ends. It will be appreciated that the results of the method 200 may be displayed, printed, recorded or automatically ported as inputs into additional calculations.
Thus, the method 200 provides an iterative process for solving the cost function (F) that eliminates the possibility of physically impossible negative elemental proportions that jeopardize the determination of the proportions of the remaining elements. A comparison of the method 200 is compared against the conventional “direct solver” approach in
It is to be understood that even though numerous characteristics and advantages of various embodiments of the present invention have been set forth in the foregoing description, together with details of the structure and functions of various embodiments of the invention, this disclosure is illustrative only, and changes may be made in detail, especially in matters of structure and arrangement of parts within the principles of the present invention to the full extent indicated by the broad general meaning of the terms in which the appended claims are expressed. It will be appreciated by those skilled in the art that the teachings of the present invention can be applied to other systems without departing from the scope and spirit of the present invention.
This written description uses examples to disclose the invention, including the preferred embodiments, and also to enable any person skilled in the art to practice the invention, including making and using any devices or systems and performing any incorporated methods. The patentable scope of the invention is defined by the claims, and may include other examples that occur to those skilled in the art. Such other examples are intended to be within the scope of the claims if they have structural elements that do not differ from the literal language of the claims, or if they include equivalent structural elements with insubstantial differences from the literal languages of the claims.
Filing Document | Filing Date | Country | Kind |
---|---|---|---|
PCT/EP2015/059024 | 4/27/2015 | WO | 00 |