GAMMA RAY LOGGING

FIELD OF THE INVENTION

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.

BACKGROUND

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 = - ₂ ∑ ₌₁ W(Ax - p) ² (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= A ^T WA ¹ A ^T Wp (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.

SUMMARY OF THE INVENTION

In a preferred embodiment, the present invention includes an improved method for determining elemental concentrations based on measurements from spectral natural gamma ray logging. In the preferred 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 (a). The values of the proportion matrix (x) are updated by subtracting the gradient factor and iterative ly 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.

BRIEF DESCRIPTION OF THE DRAWINGS FIG. 1 depicts a downhole logging system constructed in accordance with preferred embodiments.

FIG. 2 is a graph depicting standard spectra for thorium.

FIG. 3 is a graph depicting sample standard spectra for uranium.

FIG. 4 is a graph depicting sample standard spectra for potassium. FIG. 5 is a flowchart depicting a preferred embodiment of the iterative method for determining elemental concentrations from measured spectra in spectral natural gamma ray logging.

FIG. 6 is a graph depicting a measured spectra and a comparison of a fit using a direct solver method and a fit using the iterative method of the preferred embodiments.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

The preferred embodiments of the present invention include an improved method for determining elemental concentrations using spectral natural gamma ray logging. Referring to FIG. 1 , shown therein is a logging system 100 configured to carry out the analytic methods of the preferred embodiments. The logging system 100 preferably includes a wireline 102 and a multichannel sensor 104 that are disposed in a wellbore 106. The wireline 102 and sensor 104 are connected to surface facilities 108 and a computer 1 10. Although the computer 1 10 is depicted at the surface and in close proximity to the wellbore 106, it will be appreciated that the computer 1 10 may be positioned at a remote location and connected to the sensor 104 via a networked connection. Alternatively, the computer 1 10 may by embodied by a processor or computer located within the downhole portion of the logging system 100. In yet another embodiment, the sensor 104 is contained within a downhole pumping system or drilling system. In a preferred 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 1 10 representative of the measured radiation. The logging system 100 may optionally include an emitter that irradiates the formation 1 12 to produce the release of characteristic radiation from the formation.

In a particularly preferred embodiment, the sensor 104 includes about 256 channels that are each configured to measure quantity (counts) of radiation at different energies on the spectrum. FIGS. 2-4 depict the radiation spectra of thorium (Th), potassium (K) and uranium (U) for the sensor 104 and for one unit of concentration. The 256 channels of the sensor 104 are preferably selected to measure energies within the signature radiation spectra of these radioactive elements. The determination of the proportions of potassium, uranium and thorium in the formation is made using an iterative method 200 shown in FIG. 5. As explained below, the iterative method 200 solves for the elemental concentrations of potassium, uranium and thorium using a gradient method with non-negative constraints. Notably, the proposed method is superior to prior art approaches because it prevents a mathematical result that returns a negative value for the concentration of one or more of the elements.

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., le-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 a particularly preferred 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 (xo) for a proportion matrix (x) are calculated using the direct solver (equation 2). In the particularly preferred 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 (x ₀ id).

At step 206, the cost function (F) (equation 1) is evaluated with the initial values (x ₀ id) 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 (x _new ) is then calculated by subtracting a gradient factor from the current solution for the proportion matrix (x _old ). The gradient factor is defined as the product of the gradient (dF/dx) and a step size factor (a). The step size factor (a) is preferably small so that the incremental change between iterations is well controlled. In a particularly preferred embodiment, the step size factor (a) is set at 0.01. The value of the step size factor (a) 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 (x _new ) includes a negative entry. If an entry (i) within (x _new ) 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., (x _new (i) = Xoid (i))- At step 216, (x ₀ w) is then set for the subsequent iteration as equal to the determined value of (x _ne w). If, on the other hand, no entries (i) within (x _new ) are negative, the method moves directly to step 216 without the intervening step 214 and (x ₀ id) is updated for the subsequent iteration as the value of (x _new ). 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 a particularly preferred embodiment, the decision step 218 queries whether the cost function returned a result less than lxlO ^"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 a particularly preferred 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 FIG. 6. For a measured spectra 300, the conventional solution 302 yielded concentrations of 5.9501% potassium, -2.4817 ppm uranium and 2.3967 ppm thorium, respectively. The negative proportion of uranium erroneously skewed and exaggerated the presence of thorium. In contrast, the curve 304 and solution generated by the iterative method 200 more accurately reflects concentrations of 5.4274% potassium, 0.0001 ppm uranium and 0.3523 ppm thorium. This illustrates the benefits realized through the use of the iterative method 200 with non-negative constraints of the preferred embodiments.

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.