KR101292659B1 - Method for interpreting well log - Google Patents

Abstract

PURPOSE: A geophysical well-logging data analysis method reduces the time and costs required to analyze a borehole's well-logging data by analyzing the data without setting an antecedent assumption. CONSTITUTION: Raw geophysical logging data is obtained from a geophysical logging borehole (S110). As linear trends of the geophysical logging data are removed, the data becomes normalized. An empirical mode separation technique is applied to the data in order to divide the data into mode functions. C (t) is obtained by summing up values of all the mode functions and subtracting shortest wave value from the sum (S120). The laplacian of the C(t) is obtained as an attribute value LC(t) (S140). Strata boundaries are analyzed with the minimum value or the maximum value of the LC(t) or by using zero-crossings (S150). [Reference numerals] (AA) Start; (BB) End; (S100) Physical logging data; (S120) C(t) = The sum of a mode function - c1; (S130) LC(t) : The laplacian value value of C(t); (S140) Display the physical logging data and attribute value (LC(t)); (S150) Identify a sheet boundary

Description

Interpretation of physical logging data {METHOD FOR INTERPRETING WELL LOG}

The present invention relates to a method for analyzing physical logging data, and more particularly, to a method for simply and reliably analyzing physical logging data collected from boreholes.

Various geophysical exploration techniques are known for understanding the physical characteristics of the crust itself or for obtaining underground buried underground resources.

Among them, geophysical exploration through boreholes is known to be the most precise as well as universal.

At this time, in the case of the ground borehole, the drilling of several hundred meters to several thousand meters underground, and in the case of the seabed borehole after passing through several thousand meters of the seabed, and then again drilling the seabed rock layer.

Given these perforation depths (also called depths) of the boreholes and the data obtained through the massive amount of core samples obtained at the same time, it is necessary to quickly interpret the geological data collected from the boreholes to obtain desired geological information. have.

The desired stratum information may include changes in physical properties of the rocks or strata.

The geological data is also referred to as physical logging data or, more simply, logging data.

Well logging is the process of recording cracks, changes in sediment or stratification around the borehole from various data obtained during drilling of the borehole.

Data obtained by tabulating data acquired by a logging operation according to depth is called a well log.

Geo-physical well logging during logging is usually recorded by inserting a logging device (sonde, also known as sonde) into the borehole, then obtaining a response to the artificially generated physics and comparing them at each depth. Doing.

In this case, it is preferable not only to record the response result for the physical phenomenon generated artificially, but also to record the natural physical phenomenon unique to each layer (rock boundary) in the borehole at the same time.

Here, the physical phenomena include temperature, density, sound waves, electrical conductivity, specific resistance according to electrical properties, porosity of core sample, neutron, natural radiation decay data, and the like.

On the other hand, as methods for analyzing and analyzing physical log data as described above, methods used up to now can be broadly classified as follows.

(1) A method that has been traditionally used to interpret physical log data without applying any signal processing techniques.

(2) A method for detecting various periodicity that is difficult to detect when only physical log data is visually analyzed by performing Fourier or wavelet transformation on the physical log data.

(3) Noting that the obtained physical logging data exhibits fractalness, the physical logging data is not used as it is, but the holder exponent profile (Hurst exponent) or the same amount of the physical logging data ( Holder exponent profile analysis to identify heterogeneity of rock and strata.

Here, the method of (1) takes a lot of time to interpret for reasons such as the effect of noise included in the physical logging data and the subtle change of the signal itself apart from the noise. The subjectivity of the interpreter can have a lot of influence.

Therefore, in order to reduce the time required for such analysis / analysis and to control the ambiguity of analysis / analysis, various transformation techniques (for example, the Fourier transform or wavelet transform described above), in particular, physical logging data, Various techniques based on fractality have been proposed.

Among them, Fourier transform and wavelet transform are effective in obtaining the ground information when it is clear that the physical log data has stationarity and linearity, but since the actual physical log data is very complex, It is rare to meet the conditions and there are limitations in application.

Therefore, in view of the fact that the physical logging data may have non-stationary and non-linear characteristics as described above, in order to more effectively grasp the geological information, the prior art in the present invention. In Non-Patent Document 2 (Gaci et al., 2010) cited in the literature, sonic logging data are analyzed from physical logging data using a local regularity analysis technique or a local Hurst exponent analysis technique. Interpreted

However, in the case of the above-described local normal analysis technique, it has been performed under the assumption that the physical log data to be analyzed / analyzed will have fractal or multi-fractal properties.

In addition, although many physical logging data are assumed to satisfy the fractal or multiple fractal conditions, such estimates have limitations in that not all physical logging data are uniformly fractal.

Currently, the analysis and modeling of physical logging data is very widely used to obtain geologic parameters related to petrological features. Because they reflect the effects of motion, they are usually non-stationary and nonlinear.

Therefore, it was necessary to analyze a time series of data representing abnormal and nonlinear characteristics by a novel method to derive a model that better reflects geological reality.

Regarding a preferred embodiment of the present invention, the following non-patent literature may be cited as the prior art literature.

(Non-Patent Document 1) Local regularity analysis of strata heterogeneities from sonic logs, S. Gaci, N. Zaourar, M. Hamoudi, and M. Holschneider, Nonlin. Processes Geophys., 17, 455-466, 2010

quez, Geof

sica Internacional 49 (2), 55-67 (2010)(Non-Patent Document 2) Muti-scale analysis of well-logging data in petrophysical and stratigraphic correlation, E. Coconi-Morales, G. Ronquillo-Jarillo, JO Campos-Enr

quez, Geof

sica Internacional 49 (2), 55-67 ( 2010 )

(Non-Patent Document 3) Empirical mode decomposition as a filter bank, Patrick Flandrin, IEEE SIGNAL PROCESSING LETTERS, VOL 11., NO. 2, February 2004

Accordingly, the present invention can be interpreted by reflecting the unique characteristics of each strata based on only the physical logging data, excluding a priori assumption that was required in the conventional physical logging data analysis / analysis technique. It is an object to provide a method for interpreting data.

That is, according to the present invention, as the depth increases in the earth's crust, the above-described changes in physical properties such as porosity, density, and the presence of cracks and changes in sedimentation are sudden changes in measured values such as maximum, minimum, or inflection points in physical logging data. Although these features can be assumed to occur in the vicinity of where they occur, these features are often difficult to recognize without ambiguity when only the raw data, ie physical logging data, are examined, so that they can be signal-processed so that they can be interpreted better. It aims to provide a way to.

The problems to be solved by the present invention are not limited to the above-mentioned problem (s), and another problem (s) not mentioned can be understood by those skilled in the art from the following description.

In order to solve the above problems, according to a preferred embodiment of the present invention, (A) obtaining the first physical logging data from the borehole physics; (B) removing and normalizing the linear trend of the physical logging data; (C) applying an empirical mode decomposition technique that decomposes the normalized physical logging data into several mode functions; (D) a first filtering step, comprising: summing the plurality of mode functions, and then subtracting the shortest wavelength component to obtain C (t); And (E) a second filtering step, wherein the obtained Laplacian of C (t) is obtained as an attribute value LC (t).

Here, the physical logging data may include one or more of temperature, density, sound waves, electrical conductivity, specific resistance according to electrical properties, porosity of the core sample, neutron, and natural radiation decay data.

Further, in the step (D), C (t) can be obtained by subtracting the shortest wavelength component and the next shortest wavelength component from the sum of the several mode functions.

In addition, the step (E), preferably further comprises the step of analyzing the stratum boundary through the maximum value, the minimum value, and zero crossing of the LC (t).

Meanwhile, the present invention may also be provided in the form of a computer-readable recording medium having stored thereon a program capable of performing a method for interpreting physical logging data according to a preferred embodiment.

The details of other embodiments are included in the detailed description and the accompanying drawings.

Advantages and / or features of the present invention and methods for achieving them will become apparent with reference to the embodiments described below in detail in conjunction with the accompanying drawings.

However, the present invention is not limited to the embodiments disclosed below, but will be implemented in various forms, and only the present embodiments are intended to complete the disclosure of the present invention, and the general knowledge in the art to which the present invention pertains. It is provided to fully convey the scope of the invention to those skilled in the art, and the invention is defined only by the scope of the claims.

Throughout the specification, the same reference numerals refer to the same components, it should be understood that the size, position, coupling relationship, etc. of each component constituting the invention may be exaggerated for clarity of the specification.

According to the analysis method of the physical logging data of the present invention, unlike the conventional method, since physical logging data can be directly analyzed without prior assumptions, the time and cost required for the analysis / interpretation of the physical logging data are remarkably improved.

1 is a schematic flowchart illustrating a method of analyzing physical logging data according to a preferred embodiment of the present invention.
2 is a schematic cross-sectional view of a central US Appalachian basin from which physical logging data is obtained, which is used in a method for analyzing physical logging data according to a preferred embodiment of the present invention.
3 is a graph showing some mode functions to which EMD is applied in a method for analyzing physical logging data according to a preferred embodiment of the present invention.
Fig. 4 is a diagram showing LC (t) and rock boundary analyzed according to a method for analyzing physical logging data according to a preferred embodiment of the present invention.

Hereinafter, embodiments of the present invention will be described in detail with reference to the accompanying drawings.

First, in the present invention, among the data recording the various physical phenomena described above, the most suitable physical logging data is natural radiation decay data, in particular, natural gamma ray decay data, It turns out that it used as preferable physical logging data of this invention.

In the present invention, natural gamma-ray collapse data is used as physical log data because a certain amount of natural gamma rays are emitted from a rock or a specific stratum, and when data obtained by measuring the natural gamma ray by time series are used, And that it would be very useful in understanding the way in which the sedimentation environment of the stratum varies with the depth.

At this time, the rock or the strata which can obtain the natural gamma ray decay data in abundance may include a sail in the sedimentary rocks.

1 is a schematic flowchart illustrating a method of analyzing physical logging data according to a preferred embodiment of the present invention.

Referring to FIG. 1, according to a preferred embodiment of the present invention, a method of analyzing physical logging data includes: preparing a physical logging data step S100, an EMD processing step S110, a C (t) function obtaining step S120, LC (t) function acquisition step S130, and physical logging data and attribute value LC (t) display step S140.

In Fig. 1, the rock boundary identification step (S150) is a method of analyzing physical logging data according to a preferred embodiment of the present invention from the graph obtained in the physical logging data and attribute value (LC (t)) display step (S140). It is a step of obtaining a result by comparing the validity with the original data, that is, the physical logging data.

Physical Logging Data Preparation

In the physical logging data preparation step (S100), physical logging data acquired in the US Central Appalachian Basin and published by USGS in the United States was prepared.

Prior to the two filtering steps described below, linear trends were removed from the physical logging data, and also normalized to show the values with varying values in different ranges so that the physical logging data had values of 0 and 1 o'clock. It was.

After that, values formed through the processes to be described later from all data are also normalized to be conveniently compared in the same range (value between 0 and 1).

Fig. 2 schematically shows a cross section of a place where the physical log data is acquired. The contents thereof will be described later.

EMD processing steps

The EMD processing step (S110) is a step of performing EMD processing on the physical logging data.

The EMD stands for Empirical Mode Decomposition and is preferably understood as an empirical mode decomposition technique.

In the present invention, by combining the filters of the filter bank formed by the EMD, the newly filtered data is generated from the existing physical logging data and used as a method of dividing the stratum boundary (or rock boundary).

Here, it should be noted that the EMD technique is a method proposed by Huang of NASA.

The EMD technique is basically a method of decomposing a given data into several intrinsic mode functions (IMFs) (hereinafter referred to as "mode functions" in the present invention).

A method of decomposing a given data x (t) into a sum of a plurality of mode functions using EMD is as follows.

(1) identify all extremes of x (t) (including extrema, both maximum and minimum)

(2) Envelope formed by the maximum value and the envelope formed by the minimum value are obtained by using cubic spline interpolation method.

(3) The local average value m1 is obtained from the two envelopes.

(4) Define h1 by x (t) - m1 = h1.

(5) This process is repeated, Huang named this process sifting, and when the standard deviation of h formed in each iteration step is about 0.2 ~ 0.3, the process is stopped and finally The value formed by is assigned to the first mode function c1.

(6) x (t) - c1 = r1. Repeat the above steps (1) to (5) for r1 to obtain c2, c3, etc. continuously.

(7) The process of finding the mode functions stops when c becomes a monotonic function, that is, a function that always increases or decreases without a maximum / minimum value in a given interval.

C formed up to this step form a set of mode functions, and x (t) = Σc (i) + r, that is, the original data can be expressed as the sum of all the mode functions and the remaining r.

The advantage of EMD is that it does not decompose the data into previously defined functions irrespective of the given data, such as Fourier analysis or wavelet method, but rather decomposes and analyzes the data into several mode functions formed entirely from the data.

Thus, EMD can be referred to as an adaptive data processing method.

In an ideal situation, the individual mode functions represent values of components of the various periods or wavelengths inherent in the data, and in the atmospheric sciences and related fields, this is what happens in different periods of day, month, month, year, etc. We have used the mode functions obtained with EMD to interpret them.

However, there are always many unwanted noises in the data, and mode mixing, which has long wavelength components in the shortest wavelength mode function and short wavelength components in the long wavelength mode function, has been a problem. Methods have also been developed.

In the present invention, it is possible to sufficiently suppress the noise that was present in the original data by applying a method of excluding specific wavelength component (s) from a combination of individual mode functions only without geologically analyzing individual mode functions, There was no need for a separate process to solve the mode mixing problem.

In the present invention, the physical logging data x (t) can be decomposed into several mode functions by the EMD technique as described above, each mode function being the last longest wavelength component cn from the shortest wavelength component c1 of the original data (where , n is a positive integer greater than 1).

The inventors of the present invention have found that most natural gamma ray decay data can be decomposed into about 12 mode functions as a result of careful examination.

The mode function may be expressed as a set when decomposed, and thus, the set of mode functions according to the preferred embodiment of the present invention may be represented as A = {c1, c2, ..., c12}.

As mentioned in Non-Patent Document 3, mode functions that are members of A can be regarded as filters of various wavelength components, and the present invention has suggested an optimal method for filtering physical logging data using them appropriately.

In other words, the inventors of the present invention found that the effect of long-wave noise is reduced in the process of removing the linear trend from the physical logging data through various analysis processes. The disturbing component is identified as the shortest wavelength component, and as will be described later, the shortest wavelength component is removed to form a filtered signal C (t), and the Laplacian of this value is again taken to show a sharper change pattern than the original data. Attribute value LC (t) was derived.

Acquiring C (t) Function

Acquiring the C (t) function (S120), first, decomposes the physical logging data x (t) into a mode function of c1 to c12 by the EMD technique, and then c1 + c2, ..., + c12, That is, a step of obtaining the sum of each mode function, and obtaining a value obtained by subtracting the shortest wavelength component c1, or the shortest wavelength component and the next shortest wavelength component c1 + c2 from the sum.

As a result of performing numerical experiments several times, a better result was obtained when excluding only c1 from the sum of the mode functions, and C (t) mentioned in the present invention excludes only the shortest wavelength component c1.

However, it should be noted that even if we subtract both c1 and c2 from the sum of the mode functions, a relatively good result is obtained, we can subtract both c1 and c2 from the sum of the mode functions.

Thus, C (t) is C (t) = Σc (i)-c1 = x (t)-c1-r, i.e., the value of subtracting c1 from the sum of all the mode functions or the original data x (t) In (because x (t) = Σc (i) + r), it can be defined as the value obtained by subtracting c1 and the remaining r after the EMD process.

C (t) is more advantageous for applying mathematical operators by removing the shortest wavelength mode function c1 to form a smooth curve that changes much smoother than the original data, and also removes the linear trend and Since there is no r remaining after the EMD, it is an optimal filtering result.

LC (t) function acquisition step

The LC (t) function obtaining step S130 is a step of applying a second finite difference differential operator or Laplacian to the C (t) function obtained from the C (t) function obtaining step S120.

The resulting new value may be defined as LC (t), which may be referred to as an attribute value in accordance with a preferred embodiment of the present invention.

Here, the Laplacian is a second order finite difference differential operator, and when the Laplacian is used, data corresponding to a component having a long wavelength is erased and not emphasized, while data corresponding to a component having a shorter wavelength is emphasized more clearly. Effects can be obtained, and because of such characteristics, it is utilized in the field of image processing.

Display of physical logging data and attribute values (LC (t))

In the physical logging data and attribute value LC (t) display step (S140), as shown in FIG. 4, the depth (unit: ft) is displayed on the horizontal axis and normalized on the vertical axis. It is a step to display.

Rock boundary identification stage

The dark boundary boundary identification step S150 is based on the new physical property data, not based on the initial physical logging data, from the graph displayed in the physical logging data and the attribute value LC (t) display step S140. It can be seen that the strata boundaries can be more easily identified based on the defined LC (t).

Specifically, as shown in FIG. 4, the top of the formation named Sunbury shale and Berea sand, or rock boundary, that is, the stratum boundary is the local maximum / minimum value of LC (t). It can be seen that this is at the depth of appearance.

In Bedford shale, it can be seen that the stratum boundaries can be identified by the depth at which zero crossings occur.

Next, FIG. 2 is a schematic cross-sectional view of the US Central Appalachian Basin which obtained physical logging data used in the method for analyzing physical logging data according to a preferred embodiment of the present invention.

More specifically, the physiologic logging data were obtained from # 1 Windbigler borehole, Ohio, USA, at a depth of 4900 feet, where ft is in metric terms, 30.48 The data of natural gamma ray collapse up to cm) are measured in time series.

As described above, the physical logging data used in the preferred embodiment of the present invention is natural gamma ray decay data, and the natural gamma ray decay data is indicated by a curve shown next to a depth line shown vertically in FIG. 2. .

Note that FIG. 2 does not show the entire cross section of the # 1 Windbigler borehole, but only a portion of the upper layer up to about 650 feet in depth.

3 is a graph showing some mode functions to which EMD is applied in a method for analyzing physical logging data according to a preferred embodiment of the present invention.

In FIG. 3, the physical strata data is displayed in red at the bottom, and the mode function of the physical strata data is displayed in blue.

Fig. 4 is a diagram showing LC (t) and rock boundary analyzed according to a method for analyzing physical logging data according to a preferred embodiment of the present invention.

It can be seen from FIG. 4 that the strata boundaries of Sunbury shale, Berea sand, and Bedford shale are clearly marked.

Specifically, each strata boundary of the Sunbury shale layer (shown vertically) at depth 402 feet, the Berea sand layer (shown vertically) at 405 feet, and the Bedford shale layer (shown vertically) at 424 feet. They each exactly match the maximum, minimum, and zero crossings of LC (t) at the same depth (horizontal axis).

At this time, the maximum value, the minimum value, and the zero crossing are named as strata boundary indicators, respectively.

Next, Table 1 shows the depth of 4900 feet of the # 1 Windbigler borehole shown in FIG. 2, the depth at which the top of a total of 35 strata existing within the depth of 4900 feet (indicated by rock boundary strata), and the attributes described above. In the present invention, the strata boundary indicators appearing on LC (t) are tabulated.

In Table 1, the attribute value LC (t) indicating each strata boundary is a local maximum value (indicated by 'max') / minimum value (indicated by 'min') and zero crossing ('zero'). Are indicated respectively.

Also, in Table 1, the error (in feet) is zero if the depth at which the strata boundary indicators (maximum / minimum value, zero crossing) of the attribute value LC (t) and the depth of the actual strata boundaries coincide. Can be displayed.

The actual data were acquired at 0.5 foot intervals and can be defined as 'distance' to indicate errors.

For example, if the stratum boundary indicator of LC (t) appears 0.5 ft and 1 ft above the depth of the strata boundary above (denoted as' + '), +1 distance or +2 distance, or downward ('- If you appear below 0.5 feet and 1 foot, it can be displayed as -1 distance or -2 distance.

Here, 1) even if the actual strata boundary appears to be close to the maximum value (max) / minimum value (min), it was determined that the distance exceeding 2 distances (1 foot) does not coincide with the strata boundary. ) There may be cases where no distinctive features appear in LC (t) at the stratum boundary, and in both cases it is determined that the stratum boundary cannot be determined from the calculated properties and is marked as ND ("not defined"). It was.

In Table 1, it should be noted that the distance is relative distance.

From Table 1, the indicators of the strata boundary are expressed as 'max', 'min' or 'zero' by the attribute value LC (t), and these values and the error values (+1, +0.5, 0, -0.5, and -1.0), out of a total of 35 rock boundary strata, 33 strata boundaries were identified except for the two strata indicated by ND.

This level of success can mean approximately 94% accuracy.

Finally, the present invention may be implemented in the form of a computer-readable recording medium having stored thereon a program capable of performing a method for interpreting physical logging data according to a preferred embodiment of the present invention.

While the present invention has been described in connection with what is presently considered to be practical exemplary embodiments, it is to be understood that the invention is not limited to the disclosed embodiments.

Therefore, the scope of the present invention should not be limited to the described embodiments, but should be determined by the scope of the appended claims and equivalents thereof.

While the present invention has been particularly shown and described with reference to exemplary embodiments thereof, it is to be understood that the invention is not limited to the disclosed exemplary embodiments, but, on the contrary, Modification is possible.

Accordingly, the spirit of the present invention should be understood only in accordance with the following claims, and all equivalents or equivalent variations thereof are included in the scope of the present invention.

S100: physical logging data preparation step
S110: EMD processing step
S120: C (t) function acquisition step
S130: LC (t) function acquisition step
S140: Physical Log Data and Attribute Value (LC (t)) Displaying Steps
S150: rock boundary identification step

Claims (5)

(A) obtaining initial physical logging data from a physical exploration borehole;
(B) removing and normalizing the linear trend of the physical logging data;
(C) applying an empirical mode decomposition technique that decomposes the normalized physical logging data into several mode functions;
(D) a first filtering step, comprising: summing the plurality of mode functions, and then subtracting the shortest wavelength component to obtain C (t); And
(E) a second filtering step, comprising: obtaining the obtained Laplacian of C (t) as an attribute value LC (t);
The step (E) further comprises the step of analyzing the stratum boundaries through the maximum, minimum, or zero crossing of the LC (t),
How to interpret physical logging data.
The method of claim 1,
The physical logging data may include one or more of temperature, density, sound waves, electrical conductivity, specific resistance according to electrical properties, porosity of the core sample, neutrons, and natural radiation decay data.
How to interpret physical logging data.
The method of claim 1,
In the step (D)
Wherein C (t) is obtained by subtracting the shortest wavelength component and the next shortest wavelength component from the sum of the several mode functions,
How to interpret physical logging data.
delete A computer-readable recording medium having stored thereon a program capable of performing the method for analyzing physical logging data according to any one of claims 1 to 3.

Images