METHOD AND SYSTEM FOR RETRIEVING SEISMIC DATA FROM A SEISMIC SECTION IN BITMAP FORMAT

Abstract

A computer implemented method for retrieving seismic data from a seismic section provided in a bitmap format. The method includes reconstructing a two dimensional matrix of seismic interpolated data in which a value at a given pixel in the matrix is proportional to local density of wiggles in the seismic section with an added value of a previously calculated pixel. The method may be implemented in either of computer hardware configured to perform said method and computer software embodied in a non-transitory, tangible, computer-readable storage medium. Also disclosed are corresponding computer program product and data processing system.

Description

BACKGROUND

The present invention relates to seismic data. More specifically the present invention relates to retrieving seismic data from a seismic section provided in the form of a scanned seismic profile in a binary image format (bitmap).

Seismic data has been collected for several decades and a vast part of the older seismic data is only available in paper format of seismic sections. According to conservative estimates there are several millions of paper seismic sections owned by oil companies and similar institutions, the value of which is huge and sometimes immeasurable, especially when surveys are no longer possible in certain areas.

With the introduction of modern seismic workstations and user friendly interpretation software, most of these paper seismic sections have become useless, unless they can be converted into digital format (e.g., into standard SEG-Y format). In current conversion technology, the conversion process is slow, due to large amount of manual intervention, and involves relatively high costs and long throughput times. The slow process and the high costs discourage many oil companies from converting the legacy paper seismic sections into SEG-Y. The conversion is currently done on a project by project basis, driven by opportunity or to fulfill legal obligations.

A seismic section typically includes seismic traces, which represent the responses of the elastic wavefield to the density of earth layers and acoustic velocity across these layers as acoustic energy travels from a source through these layers to one or more receivers. A seismic section paper is typically a large sheet of about 2-3 meters by 1 meter. It is composed of two main parts: metadata and traces. The metadata is a graphical and a textual data describing mainly acquisition parameters. The traces area, which is the largest part, is composed of a large series of a traces plotted one next to the other. A trace is a series of amplitudes at time steps (from top to bottom), around an imaginary base line. The traces appear as vertical wiggling lines (referred to as “wiggles”). Seismic sections sometimes include horizontal time lines (or depth lines if conversion from time to depth had been performed). In some seismic traces the positive or negative sides of the wiggles are blackened—areas to the right (for positive) or to the left (for negative) of the imaginary base line about which the line wiggles confined between the line and that base line appear as black blobs.

A snippet of a larger paper scanned image of a seismic section is shown in FIG. 1. Each vertical line represents a vertical function (and is called a wiggle). It is evident that digitizing this graphical representation is not an easy task. It is evident that the vertical lines in some cases overlap and it is hard to distinguish and determine the correct line. Furthermore, typically the right hand side of a wiggle has amplitude overlaps, which are characterized by black areas (blobs), that represent the area confined between the imaginary zero axis (also called “the base line”) and the function line when the value of the line is above zero.

This makes it extremely difficult to obtain information on the positive side (or negative, depending on the “blobbed” side) of the line functions, and it is, therefore, practically impossible to follow a function line that crosses a blob.

Another problem is image quality and in particular discontinued lines (mainly associated with scan quality).

SUMMARY

According to embodiments of the present invention there is provided a computer implemented method for retrieving seismic data from a seismic section provided in a bitmap format. The method includes reconstructing a two dimensional matrix of seismic interpolated data in which a value at a given pixel in the matrix is proportional to local density of wiggles in the seismic section with an added value of a previously calculated pixel.

Furthermore, in accordance with embodiments of the present invention, there is provided a computer program product stored on a non-transitory tangible computer readable storage medium for retrieving seismic data from a seismic section provided in a bitmap format. The computer program includes code for reconstructing a two dimensional matrix of seismic interpolated data in which a value at a given pixel in the matrix is proportional to local density of wiggles in the seismic section with an added value of a previously calculated pixel.

Furthermore, in accordance with embodiments of the present invention, there is provided a data processing system for retrieving seismic data from a seismic section provided in a bitmap format. The system includes a processor; a computer usable medium connected to processor. The computer usable medium contains a set of instructions for retrieving seismic data from a seismic section provided in a bitmap format wherein the processor is designed to carry out a set of instructions to reconstruct a two dimensional matrix of seismic interpolated data in which a value at a given pixel in the matrix is proportional to local density of wiggles in the seismic section with an added value of a previously calculated pixel.

BRIEF DESCRIPTION OF THE SEVERAL VIEWS THE DRAWINGS

The subject matter regarded as the invention is particularly pointed out and distinctly claimed in the concluding portion of the specification. The invention, however, both as to organization and method of operation, together with objects, features, and advantages thereof, may best be understood by reference to the following detailed description when read with the accompanying drawings in which:

FIG. 1 is a snippet of a larger paper scanned image of a typical seismic section saved on paper or provided in the form of an image.

FIG. 2 shows isolated blobs of the seismic section of FIG. 1, after the seismic trace was subjected to morphological processing in which wiggles and time lines were removed, according to embodiments of the present invention.

FIG. 3 shows the seismic section of FIG. 1, after the time lines had been removed, according to embodiments of the present invention.

FIG. 4 is a plot of the positive values of a trace of the seismic section of FIG. 1, according to embodiments of the present invention.

FIG. 5 is a plot of the negative values of the trace of the seismic section whose corresponding positive values are plotted in FIG. 4, according to embodiments of the present invention.

FIG. 6 is a plot of the combined negative and positive values of the trace (of FIG. 4 and FIG. 5), according to embodiments of the present invention.

FIG. 7 is a plot of the trace plotted in FIG. 6 after performing smoothing, according to embodiments of the present invention.

FIG. 8 illustrates a grey-level image representing a two-dimensional matrix of recovered and interpolated seismic data corresponding to the seismic section of FIG. 1, according to embodiments of the present invention.

FIG. 9 shows a straight line drawn across a grey level image derived by performing a method for recovering seismic data from a seismic section according to embodiments of the present invention.

FIG. 10 is a plot of the seismic values present along the line shown in FIG. 9.

FIG. 11 shows a recovered seismic section in which only wiggles appear, which was extracted from and superimposed on a corresponding original seismic section (shown in grey for enhanced viewing), according to embodiments of the present invention.

FIG. 12 is a portion of a seismic section, with visible base lines (typically not present in such plots) indicating considered parameters in a mathematical approach according to embodiments of the present invention.

FIG. 13 is a flow chart of a method for retrieving seismic data from a seismic section, according to embodiments of the present invention.

FIG. 14 illustrates a computerized system for retrieving seismic data from a seismic section according to embodiments of the present invention.

DETAILED DESCRIPTION

As will be appreciated by one skilled in the art, aspects of the present invention may be embodied as a system, method or computer program product. Accordingly, aspects of the present invention may take the form of an entirely hardware embodiment, an entirely software embodiment (including firmware, resident software, micro-code, etc.) or an embodiment combining software and hardware aspects that may all generally be referred to herein as a “circuit,” “module” or “system.” Furthermore, aspects of the present invention may take the form of a computer program product embodied in one or more computer readable medium(s) having computer readable program code embodied thereon.

Any combination of one or more computer readable medium(s) may be utilized. The computer readable medium may be a computer readable signal medium or a computer readable storage medium. A computer readable storage medium may be, for example, but not limited to, an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any suitable combination of the foregoing. More specific examples (a non-exhaustive list) of the computer readable storage medium would include the following: an electrical connection having one or more wires, a portable computer diskette, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or Flash memory), an optical fiber, a portable compact disc read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the foregoing. In the context of this document, a computer readable storage medium may be any non-transitory, tangible medium that can contain, or store a program for use by or in connection with an instruction execution system, apparatus, or device.

A computer readable signal medium may include a propagated data signal with computer readable program code embodied therein, for example, in baseband or as part of a carrier wave. Such a propagated signal may take any of a variety of forms, including, but not limited to, electro-magnetic, optical, or any suitable combination thereof. A computer readable signal medium may be any computer readable medium that is not a computer readable storage medium and that can communicate, propagate, or transport a program for use by or in connection with an instruction execution system, apparatus, or device.

Program code embodied on a computer readable medium may be transmitted using any appropriate medium, including but not limited to wireless, wireline, optical fiber cable, RF, etc., or any suitable combination of the foregoing.

Computer program code for carrying out operations for aspects of the present invention may be written in any combination of one or more programming languages, including an object oriented programming language such as Java, Smalltalk, C++ or the like and conventional procedural programming languages, such as the “C” programming language or similar programming languages. The program code may execute entirely on the user's computer, partly on the user's computer, as a stand-alone software package, partly on the user's computer and partly on a remote computer or entirely on the remote computer or server. In the latter scenario, the remote computer may be connected to the user's computer through any type of network, including a local area network (LAN) or a wide area network (WAN), or the connection may be made to an external computer (for example, through the Internet using an Internet Service Provider).

Aspects of the present invention are described below with reference to flowchart illustrations and/or block diagrams of methods, apparatus (systems) and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and/or block diagrams, and combinations of blocks in the flowchart illustrations and/or block diagrams, can be implemented by computer program instructions. These computer program instructions may be provided to a processor of a general purpose computer, special purpose computer, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions/acts specified in the flowchart and/or block diagram block or blocks.

These computer program instructions may also be stored in a computer readable medium that can direct a computer, other programmable data processing apparatus, or other devices to function in a particular manner, such that the instructions stored in the computer readable medium produce an article of manufacture including instructions which implement the function/act specified in the flowchart and/or block diagram block or blocks.

The computer program instructions may also be loaded onto a computer, other programmable data processing apparatus, or other devices to cause a series of operational steps to be performed on the computer, other programmable apparatus or other devices to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide processes for implementing the functions/acts specified in the flowchart and/or block diagram block or blocks.

Flowchart/s and block diagram/s in the figures illustrate the architecture, functionality, and operation of possible implementations of systems, methods and computer program products according to various embodiments of the present invention. In this regard, each block in the flowchart or block diagrams may represent a module, segment, or portion of code, which comprises one or more executable instructions for implementing the specified logical function(s). It should also be noted that, in some alternative implementations, the functions noted in the block may occur out of the order noted in the figures. For example, two blocks shown in succession may, in fact, be executed substantially concurrently, or the blocks may sometimes be executed in the reverse order, depending upon the functionality involved. It will also be noted that each block of the block diagrams and/or flowchart illustration, and combinations of blocks in the block diagrams and/or flowchart illustration, can be implemented by special purpose hardware-based systems that perform the specified functions or acts, or combinations of special purpose hardware and computer instructions.

According to embodiments of the present invention it is proposed to employ image processing techniques, in which a bitmap image of the paper or image seismic data is converted into a grey-level image. This approach does not involve pattern recognition or any decision making, and does not necessitate locating base lines.

The description that follows below refers to performing reconstruction of seismic data according embodiments of the present invention on a seismic section which has wiggles showing both positive and negative values, but without blackened areas (blobs). Treatment of a blobbed seismic section is described further below.

Seismic sections are plot or display of seismic data along a line. The vertical direction usually refers to the arrival time (of echo) but sometimes relates to depth. Thus, two dimensional functions of x, the distance from a certain starting line S and of t, the time are considered, as depicted in FIG. 12 Seismic sections are also usually visually represented as wiggles, a bunch of one dimensional time functions along the vertical time dimension. Each wiggle swings about a typically invisible base-line. Base-lines, which have been made visible in FIG. 12, but which are usually invisible, represent zero amplitude of the wiggles which are actually functions of time t along a certain time-line x. Thus, wiggles are a conventional representation of the two dimensional function of y=f(x,t), where y is the distance of the curve from its base-line x.

A mathematical approach is hereinafter presented, according to embodiments of the present invention to reconstruct a two dimensional matrix of seismic data from the seismic section in bitmap format (22 and 24 in FIG. 13).

A two dimensional function z(x,t) defines the distance from a curve to the starting line S (see FIG. 12). Thus, equation (1) relates between those functions.

z=x+y  (1)

It is convenient to consider y(x,t) and z(x,t) as continuous functions of x and t, rather than as defined only for discrete values of x.

From (1) it follows that dy=dz−dx, where dx may be expressed in terms of dz.

According to embodiments of the present invention the seismic traces of the obtained seismic section (10 in FIG. 13) are approximated by a grey level image. A transformation may be used which translates the geometrical features of the traces into radiometric counterparts. For example, when considering densely clustered traces, bare eyes would not separate one trace from another; in this case one is expected to see a grey levels image instead. Those grey levels are roughly proportional to the distance between adjacent traces.

According to embodiments of the present invention a complementary grey scale parameter {tilde over (g)} is defined as {tilde over (g)}=1−g, so that black corresponds to {tilde over (g)}=1, and white corresponds to {tilde over (g)}=0, where g is the standard grey scale (where black is assigned 0 and white is assigned 1). It can easily be verified that grey level scale is inversely proportional to the distance between adjacent traces.

To describe seismic traces mathematically one may use a comb function. First, a regular comb is defined that corresponds to the case where the signal is mute for every time t. Hence, y≡0, and the distance between adjacent traces is always equal to Δx.

$\begin{array}{cc}\mathrm{Comb}\left(x,\Delta x\right)\equiv \sum _{k=-\infty }^{\infty }A·\delta \left(x-k\Delta x\right)& \left(2\right)\end{array}$

It is convenient to add to the definition of delta function an additional equality:

$\begin{array}{cc}{\int }_{-b}^a\delta \left(x\right)x=1,{\int }_0^a\delta \left(x\right)x=0.5,{\int }_{-b}^0\delta \left(x\right)x=0.5& \left(3\right)\end{array}$

For any a, b>0

Using this convention and applying a smoothing operator:

$\begin{array}{cc}{S}_d\left(\psi \left(x\right)\right)\equiv \frac{1}{2d}{\int }_{x-d}^{x+d}\psi \left(x^{\prime }\right)x^{\prime }& \left(4\right)\end{array}$

to comb

$\begin{array}{cc}\mathrm{Comb}\left(x,\Delta x\right)\equiv \sum _{k=-\infty }^{\infty }A·\delta \left(x-k\Delta x\right)& \left(2\right)\end{array}$

produces for 2d=n·Δx, where n>0 is an integer, a uniform grey level {tilde over (g)}₀:

$S_d\left(\mathrm{Comb}\left(x,\Delta x\right)\right)\equiv \frac{A·n}{n·\Delta x}=\frac{A}{\Delta x}={\tilde{g}}_0,$

where A≡{tilde over (g)}₀Δx, and since the smoothing operation aperture, 2d, spans always over n delta functions.

Next another train of delta functions is defined, this time not equidistant, to describe the actual trace picture, where the distance between traces may vary (along each line where t=constant). This means y should not be equal to zero (y≠0), the integration variable is z, as in Error! Reference source not found., z(x)=x+y(x).

$\begin{array}{cc}\mathrm{Train}\left(z\right)\equiv \sum _{k=-\infty }^{\infty }A·\delta \left(z-z_k\right)& \left(5\right)\end{array}$

where z_k=z(k·Δx)

Applying the smoothing operator

$\begin{array}{cc}{S}_d\left(\psi \left(x\right)\right)\equiv \frac{1}{2d}{\int }_{x-d}^{x+d}\psi \left(x^{\prime }\right)x^{\prime }& \left(4\right)\end{array}$

to Train(z) at z, where 2d=n·Δz,

$\begin{array}{cc}\Delta z=z\left(x+\Delta x\right)-z\left(x\right)\begin{array}{c}\tilde{g}\left(z\right)=S_d\left(\mathrm{Train}\left(z\right)\right)=\\ =\frac{1}{2d}{\int }_{z-d}^{z+d}\mathrm{Train}\left(z^{\prime }\right)z^{\prime }=\\ =\frac{{\tilde{g}}_0\Delta x}{n·\Delta z}·m\\ =\frac{{\tilde{g}}_0\Delta x}{\Delta z}·\frac{m}{n}\end{array}& \left(6\right)\end{array}$

Note that Δz and hence the filter aperture 2d varies, changing for each z. This time one cannot expect that these are exactly n delta functions within the aperture, 2d, of the smoothing operator, since the distances between successive delta functions are not equal. Therefore, within the smoothing operator aperture there can be m≠n delta functions. Although in average it can be expected that m and n are close to each other, one cannot reduce m and n as a common factor in

$\begin{array}{cc}\begin{array}{c}\tilde{g}\left(z\right)=S_d\left(\mathrm{Train}\left(z\right)\right)=\\ =\frac{1}{2d}{\int }_{z-d}^{z+d}\mathrm{Train}\left(z^{\prime }\right)z^{\prime }=\\ =\frac{{\tilde{g}}_0\Delta x}{n·\Delta z}·m\\ =\frac{{\tilde{g}}_0\Delta x}{\Delta z}·\frac{m}{n}.\end{array}& \left(6\right)\end{array}$

Since Δx and Δz in

$\begin{array}{cc}\begin{array}{c}\tilde{g}\left(z\right)=S_d\left(\mathrm{Train}\left(z\right)\right)=\\ =\frac{1}{2d}{\int }_{z-d}^{z+d}\mathrm{Train}\left(z^{\prime }\right)z^{\prime }=\\ =\frac{{\tilde{g}}_0\Delta x}{n·\Delta z}·m\\ =\frac{{\tilde{g}}_0\Delta x}{\Delta z}·\frac{m}{n}\end{array}& \left(6\right)\end{array}$

suggest a derivative it is convenient to choose n as small as possible. Setting n=1 and 2d=Δz provides that there is always one delta function, that is m=1. Therefore, the final result is

$\begin{array}{cc}\tilde{g}\left(z\right)=\frac{{\tilde{g}}_0·\Delta x}{\Delta z}.& \left(7\right)\end{array}$

Thus for infinitely small Δx, Δz

$\begin{array}{cc}\mathrm{dx}=\frac{\tilde{g}\left(z\right)}{{\tilde{g}}_0}\mathrm{dz}& \left(8\right)\end{array}$

From equation Error! Reference source not found. and the last result it follows that

$\begin{array}{cc}\mathrm{dy}=\mathrm{dz}-\mathrm{dx}=\frac{{\tilde{g}}_0-\tilde{g}\left(z\right)}{{\tilde{g}}_0}\mathrm{dz}=\frac{\Delta \tilde{g}}{{\tilde{g}}_0}\mathrm{dz}& \left(9\right)\end{array}$

where Δ{tilde over (g)}={tilde over (g)}₀−{tilde over (g)}

Thus

$\begin{array}{cc}\mathrm{dy}=\frac{\Delta \tilde{g}}{{\tilde{g}}_0}\mathrm{dz}& \left(10\right)\end{array}$

Without derogating generality, the proportionality factor in

$\begin{array}{cc}\mathrm{dy}=\frac{\Delta \tilde{g}}{{\tilde{g}}_0}\mathrm{dz}& \left(10\right)\end{array}$

may be set so that the derivative of y is equal to Δ{tilde over (g)} as in

dy=Δ{tilde over (g)}·dz  (11)

dy=Δ{tilde over (g)}·dz  (11)

Note, in

dy=Δ{tilde over (g)}·dz  (11)

Δ{tilde over (g)}={tilde over (g)}₀−{tilde over (g)}(z) is a function of z, rather than of x. Therefore, a numeric integration according to

dy=Δ{tilde over (g)}·dz  (11)

may be computed by a formula:

y(z_n)=Δ{tilde over (g)}·Δz+y(z_n−1)  (12)

where z_n, z_n−1 are arbitrary geometrical points on an horizontal axis, that are close one to another, and Δĝ(z) is computed between these points.

Moreover, a value y(z_n) as computed by

y(z_n)=Δ{tilde over (g)}·Δz+y(z_n−1)  (12)

refers to a signal not at point z_n, but at a point, x=z_n−y(z_n). For simplicity it may be assumed, that the signal does not change much with the point of observation x, since geological structures are also supposed to change slowly in the horizontal direction, and it may be also assumed that y does not exceed a couple of Δx. Therefore, change between x and z introduces errors of second order on dx and one may use Error! Not a valid bookmark self-reference. instead of

y(z_n)=Δ{tilde over (g)}·Δz+y(z_n−1)  (12)

y(x_n)=Δ{tilde over (g)}·Δx+y(x_n−1)  (13)

where x_n is the location of the n-th base-line and Δ{tilde over (g)}=Δ{tilde over (g)}−(x) is the deviation of the image grey level from {tilde over (g)}₀ at a point x between x_n−1 and x_n.

There is a problem with numeric integration in

y(z_n)=Δ{tilde over (g)}·Δz+y(z_n−1)  (12),

since it represents an IIR filter with a pole at |z|=1, where z is a complex variable of the z transform. Therefore, one way to cure the problem is to introduce an α<1 that multiplies the old value y(x_n−1):

y(x_n)=Δ{tilde over (g)}·Δx+α·y(x_n−1)  (14)

Otherwise the integration will amplify low frequency noises.

However, this coefficient, α, should be very close to 1, so that previous values of y at x_n−1 will be (almost) summed up. If α is much smaller than 1, the IIR filter of

y(x_n)=Δ{tilde over (g)}·Δx+α·y(x_n−1)  (14)

quickly “forgets” previous values y(x_n), so that the integration will be inaccurate.

There are other, more efficient ways to cure the problem of low frequency noise.

First, some averaging along vertical axis, t, may help to decrease the noise amplitude, at least to leave only the noise components with low frequencies both in x and in t direction. Second, a trace signal used to be a band-limited, such that it does not contain too low and too high frequencies (in t direction). Therefore, using an appropriate band-pass filter eliminates low frequencies noise and allow α=1.

A simple way to compute the grey level image corresponding to the binary wiggles image would involve using a uniform weight function

$\begin{array}{cc}\tilde{g}\left(z,t\right)=\frac{\Delta x}{2d}\sum _{i=\frac{\left(z-d\right)}{\Delta x}}^{i=\frac{\left(z+d\right)}{\Delta x}}{\mathrm{bw}}_i& \left(15\right)\end{array}$

where bw is a black (1) or a white (0) pixel, and 2d is the averaging support.

Convolution

$\begin{array}{cc}\tilde{g}\left(z,t\right)=\frac{\Delta x}{2d}\sum _{i=\frac{\left(z-d\right)}{\Delta x}}^{i=\frac{\left(z+d\right)}{\Delta x}}{\mathrm{bw}}_i& \left(15\right)\end{array}$

should preferably not be with a uniform kernel and may even be two dimensional. The kernel support should preferably be carefully chosen. In particular, while using convolution

$\begin{array}{cc}\tilde{g}\left(z,t\right)=\frac{\Delta x}{2d}\sum _{i=\frac{\left(z-d\right)}{\Delta x}}^{i=\frac{\left(z+d\right)}{\Delta x}}{\mathrm{bw}}_i& \left(15\right)\end{array}$

it may be desired that

$\frac{2d}{\Delta x}$

be an integer number, otherwise the grey level image may contain high frequency noise even in the case of y≡0, where wiggles are exactly at a distance Δx from each other. Thus, uniform kernel in

$\begin{array}{cc}\tilde{g}\left(z,t\right)=\frac{\Delta x}{2d}\sum _{i=\frac{\left(z-d\right)}{\Delta x}}^{i=\frac{\left(z+d\right)}{\Delta x}}{\mathrm{bw}}_i& \left(15\right)\end{array}$

may be disadvantageous, since in general it introduces high frequency noise for each choice of the kernel support. More advantageous would be a convolution kernel that has larger weights at its center and whose weights gradually decrease to zero at the support ends. Since convolution with the uniform kernel may be easily implemented, it is used as a standard convolution. Applying a convolution

$\begin{array}{cc}\tilde{g}\left(z,t\right)=\frac{\Delta x}{2d}\sum _{i=\frac{\left(z-d\right)}{\Delta x}}^{i=\frac{\left(z+d\right)}{\Delta x}}{\mathrm{bw}}_i& \left(15\right)\end{array}$

with 2d=Δx twice in series one gets a triangular kernel with width 2Δx. The next step is to revise an approximate equation

where z_n, z_n−1 are arbitrary geometrical points on an horizontal axis, that are close one to another, and Δ{tilde over (g)}(z) is computed between these points.

Moreover, a value y(z_n) as computed by

y(z_n)=Δ{tilde over (g)}·Δz+y(z_n−1)  (12)

refers to a signal not at point z_n, but at a point, x=z_n−y(z_n). For simplicity it may be assumed, that the signal does not change much with the point of observation x, since geological structures are also supposed to change slowly in the horizontal direction, and it may be also assumed that y does not exceed a couple of Δx. Therefore, change between x and z introduces errors of second order on dx and one may use Error! Not a valid bookmark self-reference. instead of

y(z_n)=Δ{tilde over (g)}·Δz+y(z_n−1)  (12)

where Δ{tilde over (g)} is computed in some intermediate point between x_n and x_n−1; much better approach is to replace equation

where z_n z_n−1 are arbitrary geometrical points on an horizontal axis, that are close one to another, and Δ{tilde over (g)}(z) is computed between these points.

Moreover, a value y(z_n) as computed by

y(z_n)=Δ{tilde over (g)}·Δz+y(z_n−1)  (12)

refers to a signal not at point z_n, but at a point, x=z_n−y(z_n). For simplicity it may be assumed, that the signal does not change much with the point of observation x, since geological structures are also supposed to change slowly in the horizontal direction, and it may be also assumed that y does not exceed a couple of Δx. Therefore, change between x and z introduces errors of second order on dx and one may use Error! Not a valid bookmark self-reference. instead of

y(z_n)=Δ{tilde over (g)}·Δz+y(z_n−1)  (12)

by the exact equation:

y(x_n)=y(x_n−1)+Av(x_n−1,x_n,Δ{tilde over (g)})·Δx  (16)

where Av is an averaging functional:

$\begin{array}{cc}\mathrm{Av}\left(x_{n-1},x_n,u\left(x\right)\right)=\frac{1}{\Delta x}{\int }_{x_{n-1}}^{x_n}u\left(x^{\prime }\right)x^{\prime }& \left(17\right)\end{array}$

Note that the value that the averaging functional

$\begin{array}{cc}\mathrm{Av}\left(x_{n-1},x_n,u\left(x\right)\right)=\frac{1}{\Delta x}{\int }_{x_{n-1}}^{x_n}u\left(x^{\prime }\right)x^{\prime }& \left(17\right)\end{array}$

produces coincides with the value of yet another convolution with the uniform kernel at the point

$\frac{x_{n-1}+x_n}{2}.$

Therefore, if the grey level image is obtained with a double convolution (i.e. triangular kernel), then in order to perform integration

y(x_n)=y(x_n−1)+Av(x_n−1,x_n,Δ{tilde over (g)})·Δx  (16)

one has to convolve the grey level once more with uniform kernel and pick a value of convolution at a point

$\frac{x_{n-1}+x_n}{2}.$

In other words, one can compute the grey level image with 3 sequential convolutions with uniform kernel. The resulting kernel looks like an inverted parabola at its peak and as a straight parabola at its tails.

Note, that a computation according to equation

y(x_n)=y(x_n−1)+Av(x_n−1,x_n,Δ{tilde over (g)})·Δx  (16)

is still approximation to the real signal, since a smoothing of the binary trace image is not uniquely defined operation, as it depends on the form of the smoothing kernel and the size of its support.

Thus, using the mathematical approach described hereinabove one may obtain a two-dimensional matrix of retrieved and interpolated seismic data from a seismic section provided in bitmap form (e.g. scanned from paper). This matrix is in fact a grey-scale image of the retrieved seismic data.

The mathematical approach described above assumed that the seismic section only included wiggles.

However typically seismic sections may also include time lines and blobs. Therefore, in order to retrieve seismic data from such seismic sections the time lines and blobs should be addressed too.

According to embodiments of the present invention it is therefore suggested to disregard the time lines. This may be achieved, for example, by removing the time lines (14, 12 in FIG. 13) from the seismic section for example, by employing morphological processing which remove the time lines form the scanned seismic section, but may also be achieved in alternative ways. FIG. 3 shows the seismic section of FIG. 1, after the time lines had been removed, according to embodiments of the present invention.

When considering a seismic section with blobbed parts (16 in FIG. 13)—typically the positive sides of the wiggles, but may also be the negative sides of the wiggles—one may treat the blobbed sides separately from the non-blobbed sides of the wiggles.

According to embodiments of the present invention, it is suggested to regard the blobs separately by isolating the blobs from the seismic section. For example, this can be done in morphological processing.

FIG. 2 shows isolated blobs of the seismic section of FIG. 1, after the seismic section had been subjected to morphological processing in which wiggles and time lines were removed, according to embodiments of the present invention.

After isolating the blobs pixel values of the entire image of isolated blobs are saved in a two-dimensional matrix where a first constant value (e.g. 0.5) is assigned to pixels located in blobs (18 in FIG. 13).

FIG. 4 is a plot of the positive values of a trace (vertical line) of the seismic section of FIG. 1, according to embodiments of the present invention, retrieved from an image with isolated blobs. It should be mentioned that as scans of seismic data may include artifacts attributed to noise, dirt or other accidental factors that may affect the scan, artifact removal techniques may be used to enhance the final results. For example, one may employ morphological processing such as “salt and pepper” removal (black artifact pixels in white areas and white artifact pixels in blob zones). Other artifact removal techniques may be employed too.

FIG. 5 is a plot of the negative values of the trace of the seismic section whose corresponding positive values are plotted in FIG. 4, according to embodiments of the present invention. These values represent the values derived in the mathematical approach described hereinabove, according to embodiments of the present invention (20 in FIG. 13.

FIG. 6 is a plot of the combined negative and positive values of the trace (of FIG. 4 and FIG. 5, and see 22 in FIG. 13), according to embodiments of the present invention.

FIG. 7 is a plot of the trace plotted in FIG. 6 after performing smoothing, according to embodiments of the present invention. Smoothing (26 in FIG. 13) may be useful in order to limit the bandwidth of the retrieved signal, as, in fact, the positive sides of the trace the blobbed areas were assigned the value of 1, whereas in the negative sides actual values were retrieved. Smoothing is performed across the entire two-dimensional matrix.

A band pass filter may also be performed on the data (30 in FIG. 13), in a one-dimensional manner (along the time coordinate axis).

FIG. 8 illustrates a grey-level image representing a two-dimensional matrix of retrieved and interpolated seismic data corresponding to the seismic section of FIG. 1, according to embodiments of the present invention. In compliance with the sampling theorem (also known as the Nyquist-Shannon theorem), as the retrieved signal is band-limited, it is possible to reconstruct the entire field of the seismic section and obtain seismic values at any given pixel in that field, and thus all seismic values of the two-dimensional matrix forming the grey-level image shown may be interpolated.

FIG. 9 shows a straight line drawn across a grey level image derived by performing a method for retrieving seismic data from a seismic section according to embodiments of the present invention.

Once the entire seismic data field has been retrieved it is possible to retrieve seismic profiles at any given direction in the two-dimensional field.

FIG. 10 is a plot of the seismic values present along the line shown in FIG. 9. This is the seismic profile along that line.

FIG. 11 shows an extracted seismic section in which only wiggles appear, which was extracted (32 in FIG. 13) from and superimposed on a corresponding original seismic section (shown in grey for enhanced viewing), according to embodiments of the present invention. Since the original traces in the original seismic section represent signals received by physical receivers whose locations are known it may desired, when reconstructing the seismic section, to present traces which correspond to the actual traces found in the original scanned seismic section. However, any other way of trace presentation may be derived form the two-dimensional matrix of reconstructed seismic data.

The fact that the two-dimensional matrix of reconstructed seismic data covers the entire field makes it possible to retrieve traces in any desired manner.

For example, it may happen that the scanned seismic section was not properly aligned and is therefore tilted. If this is determined the recovered seismic section may present traces along tilted lines in the recovered field.

The scanned image of the seismic section may be distorted, for example, due to non-parallel feeding of the paper on which the seismic section is provided into the scanner. In such a case and upon determining the distortion, the recovered seismic section may include traces which represent profiles along distorted lines whose distortion corresponds to the determined distortion of the scan.

FIG. 14 illustrates a computerized system for recovering a seismic section according to embodiments of the present invention. The system, such as, for example, a computer device, may include an input/output unit 60, which may include, for example, a keyboard, a pointing device, a touch screen, a printer and a monitor. The system also includes processing unit 62, which may include a single processor or a plurality of processors, storage medium 64, which may include, for example, a hard-disk, flash memory, floppy disk, disk-on-key, on which a computer executable program may be stored, which includes instructions to perform a method according to embodiments of the present invention. A communication unit 66 may be provided for communicating with another system across a communication line or a network over communication wire or wires or wirelessly.

Claims

1. A computer implemented method for retrieving seismic data from a seismic section provided in a bitmap format, the method comprising reconstructing a two dimensional matrix of seismic interpolated data in which a value at a given pixel in the matrix is proportional to local density of wiggles in the seismic section with an added value of a previously calculated pixel.

2. A computer implemented method as claimed in claim 1, comprising removing time lines from the binary image.

3. A computer implemented method as claimed in claim 1, comprising assigning a constant value to pixels which form a blub in the bitmap.

4. A computer implemented method as claimed in claim 3, comprising combining the pixels assigned with the constant values with reconstructed two dimensional matrix.

5. A computer method as claimed in claim 1, comprising smoothing the seismic data of the two dimensional matrix.

6. A computer implemented method as claimed in claim 1, comprising extracting a seismic section that includes traces.

7. A computer method as claimed in claim 1, comprising producing a grey-level image from the two-dimensional matrix.

8. A computer program product stored on a non-transitory tangible computer readable storage medium for retrieving seismic data from a seismic section provided in a bitmap format, the computer program including code for reconstructing a two dimensional matrix of seismic interpolated data in which a value at a given pixel in the matrix is proportional to local density of wiggles in the seismic section with an added value of a previously calculated pixel.

9. A computer program product as claimed in claim 8, wherein the computer program includes code for removing time lines from the binary image.

10. A computer program product as claimed in claim 8, wherein the computer program includes code for assigning a constant value to pixels which form a blub in the bitmap.

11. A computer program product as claimed in claim 10, wherein the computer program includes code for combining the pixels assigned with the constant values with reconstructed two dimensional matrix.

12. A computer program product as claimed in claim 8, wherein the computer program includes code for smoothing the seismic data of the two dimensional matrix.

14. A computer program product as claimed in claim 8, wherein the computer program includes code for extracting a seismic section that includes traces.

15. A computer program product as claimed in claim 8, wherein the computer program includes code for producing a grey-level image from the two-dimensional matrix.

16. A data processing system comprising: a processor; a computer usable medium connected to processor, wherein the computer usable medium contains a set of instructions for retrieving seismic data from a seismic section provided in a bitmap format wherein the processor is designed to carry out a set of instructions to reconstruct a two dimensional matrix of seismic interpolated data in which a value at a given pixel in the matrix is proportional to local density of wiggles in the seismic section with an added value of a previously calculated pixel.

17. A data processing system as claimed in claim 16, wherein the set of the instructions include instructions to remove time lines from the binary image.

18. A data processing system as claimed in claim 16, wherein the set of the instructions include instructions to assign a constant value to pixels which form a blub in the bitmap.

19. A data processing system as claimed in claim 18, wherein the set of the instructions include instructions to combine the pixels assigned with the constant values with reconstructed two dimensional matrix.

20. A data processing system as claimed in claim 16, wherein the set of the instructions include instructions to smooth the seismic data of the two dimensional matrix.

21. A data processing system as claimed in claim 16, wherein the set of the instructions include instructions to extract a seismic section that includes traces.

22. A data processing system as claimed in claim 16, wherein the set of the instructions include instructions to produce a grey-level image from the two-dimensional matrix.