Seismic data processing

Abstract

A method of horizon interpretation comprises: (i) obtaining seismic data for a subsurface region; and (ii) analysing the seismic data to interpret one or more horizons within that data. Step (ii) comp

Description

The present invention relates to the field of seismic data processing. In particular, it relates to a method of processing seismic data to identify or interpret horizons within or from the data. It is known to collect seismic data for a subsurface region and then process or analyse that seismic data in order to identify horizons or surfaces in the subsurface region. In this context, a horizon is, or represents, a boundary or surface between two layers or regions of different materials (e.g. different types of rock, water, gas, hydrocarbon, etc.) in a subsurface region. By locating and identifying such horizons in seismic data, a three-dimensional (3D) model of the subsurface region can be created, using the identified horizons. Existing horizon interpretation is typically performed by a computer executing a "point-to-point" exercise, followed by interpolation and smoothing of the result. This is generally a non-reversible process that gives mediocre results in terms of the accuracy of the identified horizons. One such example of this is disclosed in US 5056066 A. US 5056066 A discloses a method for tracking seismic events (horizons). The method tracks such events in a two-dimensional slice of a 3D seismic data volume. The method designates a starting data point on the seismic event and tracks the event through the data grid by sequentially establishing tiles of data, defined by data points of the grid, about the starting point. Each data point is then tested to see if it meets an acceptance criterion for the seismic event. If it does, such data points of the tile are stored as identifying the seismic event. Next, each of the data points of a previously accepted tile are used as a starting data point about which a new tile is defined. The process is repeated until no more tiles are available for testing. The accepted data is then displayed so as to distinguish the accepted data from other data of the grid thereby identifying the seismic event (horizon). However, there is a desire for a more accurate method of horizon interpretation. According to a first aspect, there is provided a method of horizon interpretation, the method comprising: (i) obtaining seismic data for a subsurface region; and (ii) analysing the seismic data to interpret one or more horizons within that data. Step (ii) comprises modelling the horizon as a surface with at least two different kinds of forces acting on the surface, and determining the horizon for which: the net forces acting on the surface are zero, and / or an energy function representing the horizon is minimised. Obtaining seismic data for a subsurface region may comprise obtaining such data from one or more memories or data storage media. Alternatively or additionally, it may comprise performing a seismic survey to obtain such data. The seismic survey could be performed in any known way. Optionally, after the seismic data has been obtained, it may be (pre)processed, e.g. to remove noise, before step (ii) is performed. Once a horizon has been determined, preferably data representing that horizon is stored in a memory, e.g. in the form of an array. The determined horizon may be displayed graphically, e.g. on a screen, e.g. in two or three dimensions. The surface preferably comprises an array of (data) points, with each point preferably corresponding to a location in the seismic data. As such, a horizontal spacing between the points may be constant. The at least two different kinds of forces acting on the surface may act in opposing directions. The at least two different kinds of forces acting on the surface may comprise a first kind of force and a second kind of force, the first kind of force preferably being an internal force to the surface and / or the second kind of force preferably being an external force to the surface. The first kind of force may represent a force that acts to flatten the surface (e.g. reduce any curvature in the surface) and / or the second kind of force may represent a force that acts to pull the surface towards maxima in the seismic data. Preferably, the at least two different kinds of forces acting on the surface act only in a vertical direction. The surface may be modelled as a blade spring. As such, the first kind of force preferably representing a force that acts to flatten the surface could be in the form of a force from a blade spring. Preferably, at least one (and preferably both) of the at least two different kinds offerees acting on the surface are modelled as linear spring forces, e.g. with spring constants. The method may further comprise tuning a magnitude of at least one of the at least two forces acting on the surface. Such tuning may be performed for the surface as a whole and / or separately / differently for different locations within the surface. For example, in locations where data quality is low, a magnitude of at least one of the at least two forces acting on the surface (e.g. a force acting to pull the surface towards a maxima in the seismic data) may be tuned to be lower. This is because in regions where data quality is low, the data may be less accurate and so “pulling” the surface towards a maxima in the seismic data may be less desirable. Preferably, the method further comprises determining whether the inclusion of any cuts or gaps in the horizon (e.g. between points in the surface) would result in a reduction in an energy function representing the horizon and, if so, which cuts or gaps in the horizon would minimise the energy function. This could involve iteratively testing pairs of data points to determine if and where any cuts would minimise the energy function. Such cuts or gaps in the horizon could indicate the location of faults in the subsurface region. As such, the locations of such cuts or gaps in horizon are preferably stored in a memory, e.g. as representing cuts or gaps in the horizon. Alternatively, such cuts or gaps could be used to identify a boundary or edge of a horizon. The energy function may be determined based on energy in the springs in the model of the surface. Preferably, the energy function includes an energy value representing an energy cost (e.g. a “breakup energy”) when a cut is made in the horizon. As such, identifying an excessive number of (false) faults in the region may be avoided, as the energy cost associated with this would prevent this from minimising the energy function. The method preferably comprises tuning a magnitude of the energy value representing an energy cost when a cut is made in the horizon, e.g. such that a compromise is met between identifying true faults and misidentifying false faults. Step (ii) may comprise modelling the horizon as surface with at least three different kinds of forces acting on the surface, wherein at least one kind of force acting on the surface represents a force pulling the surface towards a known object in the subsurface, such as a well marker. The method may further comprise tuning the magnitude of that force, for example as a function of horizontal position. For example, a force constant such as a “spring constant” representing the strength of the force pulling the surface towards a known object may be adjusted as a function of the horizontal distance to the known object. For example, the force constant could be proportional to the inverse of the horizontal distance to the known object. According to a further aspect, there is provided a method of modelling a subsurface region, the method comprising performing the method described herein (with any of its optional or preferred features) to interpret at least one horizon from seismic data, and using the interpreted horizon as at least one input for modelling the subsurface region. According to a further aspect, there is provided a method of determining a trajectory or location of a well, the method comprising performing the method described herein (with any of its optional or preferred features) to interpret at least one horizon from seismic data and / or model a subsurface region, and using the interpreted horizon and / or modelled subsurface region to determine the determine the trajectory or location of the well. According to a further aspect, there is provided a method of drilling a well, the method comprising performing the method described herein (with any of its optional or preferred features) to determine the trajectory or location of the well and then drilling the well with that trajectory or location. According to a further aspect, there is provided a computer program product comprising computer readable instructions that, when executed by one or more computers, cause the one of more computers to perform any of the methods described herein (with any of its optional or preferred features). According to a further aspect, there is provided a system comprising one or more computers and one or more memories, the system being configured to perform the method described herein (with any of its optional or preferred features). Thus, a method may be performed that can provide improved and more accurate interpretation of horizons (and faults) within seismic data. More accurate knowledge of the location of such horizons can have various applications and benefits. For example, horizons identified from seismic data may be used: - to build models or 3D representations of reservoirs; in the decision-making process for the drilling of a well, e.g. for deciding, where, how far, and / or in which direction(s) / with which trajectory to drill a well; - to locate a hydrocarbon reservoir more accurately, e.g. in a location where there is already a known / existing oil field; - to locate a fault more accurately, e.g. such that better, more informed decisions may be made. Preferred embodiments of the invention will now be described by way of example only and with reference to the accompanying figures, in which: Fig. 1 is a chart showing an exemplary slice of seismic data and an enlarged portion thereof; Fig. 2 is a chart showing an exemplary vertical section of seismic data.; Fig. 3 is a chart showing an exemplary vertical section of synthetic seismic data with noise with an initial interpretation of a horizon; Fig. 4 is a chart showing the exemplary vertical section of synthetic seismic data with noise of Fig. 3 with an improved interpretation of the horizon; Fig. 5 shows a model of a horizon surface with “curvature forces” acting on each point; Fig. 6 shows a model of a horizon surface with “seismic forces” acting on each point; Fig. 7 shows a model of a horizon surface with a force pulling it to the level of a well marker; and Fig. 8 shows a model of a horizon surface from seismic data in three dimensions. A method of interpreting horizons in seismic data is provided. This method (or at least step (ii) thereof may be provided in the form of computer executable instructions (a computer program(s)) stored on a software product(s), and may be executed by a computer(s). The method involves the following steps (although further steps may also be performed as described herein): (i) obtaining seismic data; (ii) analysing the seismic data to interpret one or more horizons in the seismic data; and (iii) using the one or more horizons as interpreted in step (ii) in a model and / or as the basis for a further decision. At step (i), seismic data is obtained for a subsurface region. The subsurface region is one about which it is desired to have information (or better / more accurate information) about one or more horizons or surfaces that may be present in the subsurface region. The step of obtaining seismic data may involve simply obtaining or retrieving the seismic data for one or more data storage means or memories. Altematively or additionally, the step of obtaining seismic data could involve performing one or more seismic surveys to detect, record and store the seismic data in a memory. Such seismic surveys could be performed in any known way for performing such a survey. The seismic data may then be processed (or pre-processed) according to known techniques, e.g. for noise removal. The seismic data is then ready for further analysis and, specifically, the method of horizon interpretation described herein. At step (ii), analysing the seismic data is analysed to interpret (determine / identify / locate) one or more horizons in the seismic data. In this method, horizon interpretation may be performed as a kind of (or in a similar / equivalent way to an) “energy minimisation” problem. In other words, the horizon surface is modelled as a physical surface that is being constrained in its movement by internal forces holding the surface together (referred to as “curvature forces” herein), and at the same time being pulled on by external forces (referred to as “seismic forces” herein), making it possible to tune the degree of smoothness. The horizon can be identified or interpreted by finding the surface that has zero net force acting on it, such that it is in equilibrium, and / or that minimises its (potential) energy, the two being essentially equivalent. In this method, the horizon is modelled as a stiff mesh with elastic properties, similar to a steel mesh used to reinforce concrete. Such a mesh is elastic so it will resist bending. However, if the forces are strong enough, it will bend or even break, depending on the situation. As such, this mesh model can be used to model smooth or flat horizons, e.g. by removing noise from the seismic data, and also breaks or discontinuities (e.g. faults in a subsurface) in a horizon. Such a model is physical, not just mathematical, which means that terms from physics, such as forces, energy etc, may be used to describe it The model (or its parameters) may then be tuned to suit particular needs of horizon interpretation. In this model, the horizon is considered to take the form of a physical object or surface that consists of a set of points forming that surface. Each point is a pixel in a map view. The lateral distance between points is given by the distance between seismic survey inline (IL) datapoints and crossline (IL) datapoints. Thus, for a typical seismic data set, the points are separated by 12.5 m, although this may vary. This is illustrated in Fig. 1, in which the chart on the left shows a horizontal slice of seismic data and the chart on the right shows an enlarged portion of that data. In the enlarged portion on the right, a single (data) point is indicated as a square pixel. A typical seismic survey has a few thousand IL and XL data point coordinates, which means that a horizon in an extensional tectonic setting (no overlap in time) would typically contain somewhere in the range of a million to a billion data points. Fig. 2 shows an exemplary vertical section of seismic data. The (red, white and blue) shading indicates the seismic data. A horizon 1 is marked over the seismic data with a solid black line (for ease of illustrating the principles, this chart is in two dimensions so the horizon is a line not a surface). The horizon 1 contains around 100 points of seismic data. The line 1 has been drawn from point to point so that it appears as a continuous line. The goal is for the horizon to represent the underlying reflection that the seismic is imaging, not to fit perfectly to every small wiggle in the seismic data, since there is noise in the seismic data. As such, it is desirable to smooth out the horizon to minimise the impact of noise. However, it is important for this smoothing not to be too extreme and remove true deviations, or even discontinuities / faults, in the data, whose existence and location it is important to determine. Such smoothing can be achieved by assigning each (possible / candidate) point of the horizon an amount of “energy” (i.e. a quantity which can be considered as if it were their energy in this model). This energy is based on the position of the horizon point relative to the seismic data (specifically relative to the position of a maximum in the seismic data, that could indicate the location of a reflector surface), relative to the other (neighbouring) points in the horizon, and, in some cases, relative to other objects in the subsurface as well, such as wells. The points can then be allowed to change their vertical position iteratively so that a final vertical position for each point can be found, which is believed to correspond most closely to the true location of the horizon. The final vertical position of each point is the position that gives the lowest “energy” to the system (i.e. the grid or mesh of points, taking into account any “forces” acting on it) as a whole. There are various ways in which the energy functions modelling the system can be chosen and tuned, which may depend on the desired behaviour of the system. In one example, energy terms from physics are used, although this is not a requirement. In the model, the horizon points can only be influenced by the energy functions that are used to describe the system. Thus, if it is desired to smooth the interpreted horizon (krieging), this is done with the energy functions. This principle is illustrated in Figs. 3 and 4. Fig. 3 shows an example vertical section of synthetic seismic data with noise. An initial interpretation of a horizon 2 is indicated by a series of adjacent straight lines, connecting the data points that have been interpreted as forming the horizon 2. However, a small fault 3 is clearly visible in this section but it has not been accounted for in the horizon 2, which just continues smoothly over it. The way in which the initial interpretation of the horizon is performed is not critical to the rest of the method. For example, it may simply require that it is possible to measure or detect the location of a peak in the seismic data. As the method is able to perturb the surface by adding points, it is not strictly necessary to have any initial interpretation. All that is required is a single point on a peak. However, it can be beneficial to start from an initial interpretation using a method such as described in US 5056066 A or from a patch (e.g. as described in GB 2583910B) as this can help to avoid the need for performing many time-consuming perturbations. On the other hand, it may also be that it is not desired to let the method add or remove points. In such cases, the method may involve just adjusting the time / depth of existing points and a prior interpretation is not needed. In any case, it is conceptually easier to explain and understand the rest of the method when starting from an existing (initial) horizon interpretation so for clarity of description and ease of understanding, an initial interpretation is indicated here. The objective of the horizon interpretation method is to find the true or most accurate representation of the horizon or surface, as indicated by the horizon 4 in Fig. 4. This horizon 4 is smooth because noise in the seismic data has been accounted for but it also has a cut (discontinuity) at the fault 3. As such, it constitutes an accurate representation of the true surface. This objective may be achieved with (only) two parameters: one that governs how smooth the output horizon is; and one that governs how large faults should be in order to cut (i.e. form a break in) the horizon. Prior knowledge of faults (e.g. whether any are present and where) is not required. The method may detect whether there are any such faults in the horizon and their location(s), as described below. However, in some cases, if there is some a priori knowledge about this, it can be used to improve the method. This is also described below. Finding, with an automated method, a horizon that both is smooth where it should be and includes any true faults is not possible with known horizon interpretation techniques. In such known techniques, the horizon is first snapped to the peak in the seismic data. Then, to account for any noise in the data, a smoothing process is applied which typically does not take the seismic data itself into account when it is performed. This smoothing process does not have any knowledge about where to cut the horizon (e.g. the presence or location of any fault), so in order to achieve a smooth output (without noise), it will simply smooth the horizon across any faults that may be present, leading to poor results. Alternatively, the faults may be pre-identified, with cuts formed in the horizon at the pre-identified points. However, there is no known method that itself can perform both smoothing and the identification and locating of any faults. The present method, which can both smooth and identify and locate any faults in a horizon, achieves this with a model involving two “forces”: a “seismic force” that pulls each horizon point towards the position of the maximum peak in the seismic data (for a given pair or horizontal coordinates); and - a “curvature force” that pulls each horizon point towards a horizon surface with zero curvature (i.e. to lie on a straight line, or flat surface, with all the other points forming the horizon). In the model, the system is modelled as a physical surface. The only way the surface can be influenced or affected is by giving it energy. The system is of course not really a physical surface, rather it is a set of locations in the data, which together can be taken to form a surface. The “energy functions” used in the model may be selected in various ways. Suitable types of energy functions are preferably functions with few parameters, functions that have physical analogues, and / or functions that are easy to understand in terms of how they should be adjusted (tuned) to get the desired results. Fig. 5 shows a model of a horizon surface 5 in two dimensions. In this model, each circle or point 6 represents a single data point of the horizon surface 5. The distance x is the horizontal trace separation (typically 12.5 m). This distance x is fixed. There are only vertical forces, as indicated by the single-ended arrows in Fig. 5, acting on the data points 6. These (internal) forces are pulling the points 6 towards a flat, zero-curvature surface 7. With no external forces acting on the horizon 5, the horizon 5 wants to minimise its energy and be a zero-curvature surface 7. To describe this in a little more detail, the horizon surface 5 connecting the points 6 can be considered as a blade spring (e.g. like a plastic ruler that can be bent elastically). Each point 6 can be considered to have force acting on it (indicated by the three single-headed arrows in Fig. 5) by the “blade spring” horizon surface 5, pulling it towards the zero-curvature surface 7. Considering this as a static scenario, the sum of the forces acting on the points 6 has to be zero. This means that the “blade spring” pushes the two outer points 6 up, each with half the force that is pushing the central point 6 down. Such forces may be referred to as “curvature forces”. The fundamental force set up by curvature (the curvature force Fc) for a point i is given by: Fc. = -k.Ax I I where Ax is the deviation from the zero-curvature situation and k is the spring constant. Thus: Fc. = -k. (x -1 / 3 (x. + x. + x. J) or i i ' i ' i-1 i 1+1" Fc.= 1 / 3 k. (X.-2X. + X. J Other exponents of x would be possible, but a force that is linear with displacement is a classical spring and gives a linear system, and is thus a good place to start. The sum of forces acting on a system is zero for a stationary system: a spring that is pulling a point down, is at the same time pulling its neighbours up. Each neighbour thus experiences an opposite force with half the magnitude. Thus, the total force experienced by a point or element i is: Fc. = 1 / 3 k.(x„ -2x. + x. J - 1 / 6 k. <(x. „ -2x. + x) - 1 / 6 k. <(x. -2x. + x ) I I ' 1-1 I 1+1' 1-1 ' 1-2 1-1 f 1+1 ' I 1+1 t+Z (Eq. 1) It is possible to work with different values of k. (spring constants, equivalent to a “penalty” for a point being away from a straight line) for each point. However, if we consider only one value for all points, i.e. k. = k (same for all), a simpler scheme is provided. In this blade spring model, the spring (point) needs to be connected to a neighbour on each side to result in a force. Thus, k. and kn are both zero. The full schema becomes (using a factor of 6 to eliminate fractions): 6 Fc1 = - k (x1 -2x2 + xj 6Fc2 = 2k (x1-2x2 + x^ - k (x2-2x3+x^ = k (2x1 -5x2+4x3~xJ 6 Fc. = 2 k (x^ -2x. + - k (x. 2 -2x.^ + x} - k (x. -2xM + x^ = k (- x+4x- 6x + 4x.. -x.J, i >2 with similar adjustments for elements n and n-1 as for 1 and 2. There are no forces if there are only one or two elements as they naturally can only form a point (one element) or a straight line (two elements). It can be beneficial to keep the formula in the form of Eq. 1 above for calculation purposes. Now let us consider the modelling of the so-called “seismic forces”. Fig. 6 illustrates a model of the “seismic forces”, illustrated by a series of coil springs 13, pulling the horizon 11 towards the seismic data points, whose locations are indicated by the thick black line 12. For simplicity and linearity, this is modelled as a series of coil springs 13 although other models of these forces are also possible. In this model, the seismic data (indicated by line 12) is “pulling” each point 14 of the horizon 11 towards the position at which it would end up in a “snap-to-peak” exercise, i.e. if the horizon 11 were to correspond exactly to the peaks in the seismic data. The strength of the pull (spring constant) on each point 14 can be set to vary. In an example, it may be set to be proportional to the amplitude of the peak for each horizontal location. The “seismic force” acting on the horizon 11 thus becomes: Fs. = ks. * (x. - xs.) where xs is the location of the peak in the seismic data and ks.is the amplitude of the seismic data (or 1 if the amplitude of the seismic data is taken into account). As only the equilibrium position is being modelled, the sum of the “curvature” and “seismic” forces is zero everywhere. Thus, the total system is fully described by a set of n equations with n unknowns. Fs + Fc. = 0 I I As the ks. terms are obtained from the seismic data, there is only one tuning parameter in this set-up, namely the curvature force constant k. k can be selected to fit a desired compromise between fitting to the seismic data (snap to fit) and obtaining a low-curvature solution (reducing the impact of noise). If k has a value of zero, this results in a classical “snap-to-peak” horizon output. On the other hand, if k is infinite (which is easiest to set up by setting all values of ks.to zero, which is effectively the same thing, as it is the ratio k / ks that matters) this results in the horizon being a straight line. The preferred, selected, or optimum, value of k may be a matter of taste and data. Generally, it is desired for the surface to appear connected (except in any areas with breaks or discontinuities). This means that the ratio ks / k is generally (i.e. in areas without breaks or discontinuities) less than 1 as it is more important or desirable for a point to have the same (or close to the same) time / depth as that of its neighbour(s) than to have a time / depth equal to that of the seismic data. On the other hand, ks may vary laterally. In one example, the equations set out above are solved iteratively (e.g. using a computer program) by calculating forces and moving the points in each timestep until convergence is reached. In each timestep, the elements are moved a small distance proportional (e.g. with a proportionality constant equal to 1) to the sum of the forces acting on them. The distance should be small so as to avoid getting a numerically unstable method. Since it is only the ratio k / kthat really matters, we are free to choose the proportionality constant to be 1 and to just tune k to not get numerical instability. We then get a starting value for the set of kss by finding the medial seismic amplitude and setting the initial proportionality constant for the kss to be such that the median ksdivided by k is 1 / 1000. This gives a solution that is a quite tight fit to the seismic. The kss are then reduced by reducing the proportionality constant to make the solution smoother. For simplicity and ease of understanding, the above explanation was in the context of two dimensions. However, it is straightforward to extend the same principles to a system in three dimensions. An example of a horizon surface 8 made up of a number of points 6 and connections 9 therebetween, identified in three dimensions, is illustrated in Fig. 8. As well as fitting a horizon to data whilst smoothing it to eliminate the effect of noise (by using the “seismic” and “curvature” forces, as described above), it is also important to determine whether there are any breaks or faults in the horizon and their location. A fault would be represented in the model by a gap 10 or broken connection between two points 6 in the horizon. It is easy to define broken connections by treating these the same way as edges, i.e. setting k and k.+1 to 0 to break the bond between them. Then, no forces are transmitted between the points i and i+1, and the vertical distance between them can be anything. For a given set of points, and information on which connections are “active” (i.e. not broken), it is possible to solve for the optimum (lowest energy) horizon surface. Breaking a bond is an operation similar to adding or removing a point. The information about which bonds are broken may be provided as external output. This information may be fixed, e.g. so that the method cannot connect or disconnect. Alternatively, as described herein in relation to energy perturbation, the method itself may perturb the system to find the best solution, in which case a “breakup energy” is assigned to each connection. Opting for a linear force in x, (F = kx) the system is linear and can be solved by matrix inversion. For other choices, the system can be solved numerically. The result is a surface with low curvature, that adapts to the seismic data and is able to handle sharp breaks such as faults. As such, it provides a very good representation of real-world layer boundaries, and should be well-suited for horizons that are used as input to e.g. geomodels. By adjusting the parameters, such as k and / or ks., in the model, the fit to seismic data can be tuned on a location-by-location basis. For example, close to faults, the seismic reflector position data can often be less reliable due to the imaging conditions in such locations. As such, it may not be desirable for the horizon interpretation to follow the fault surface is that is imaged there, which can extend steeply upwards or downwards, in some cases. If it is actually desired to have a (more) planar horizon surface in regions adjacent to a fault (assuming their location is known), the seismic “spring constant” may be scaled down by a factor that could depend, for example, on the distance to the fault. A fault would correspond to a broken bond. In one example, the effective ks could be set to 0 at all points within a particular radius from a broken bond. For points within such a radius, this would mean that the seismic data for such points is ignored when determining the horizon. In another example, ks could be multiplied by a factor that is less than 1 and depends on the distance to the broken bond (the factor being lower for points closer to the broken bond, and increasing as points are further from the broken bond). In other cases, if it is desired to force the horizon surface to pass through a specific location, e.g. a well marker, an additional spring force may be introduced, acting on the horizon, that pulls the horizon surface to that level. Such a force is illustrated in Fig. 7 by the springs 24 pulling the upper four points 23 towards the well level 20. The line 20 indicates the level of a well. The horizontal location of the well is indicated by line 21. Positions without adjustment to the well level 20 are indicated by the lower four circles 22. Positions with adjustment to the well level 20 are indicated by the upper four circles 23. At the horizontal location of the well 21, the spring is not visible as it is a spring of infinite strength pulling the point right up to the level 20 of the well. In this example, a spring constant for such a force, i.e. pulling the points towards the well level 20, could be selected that varies proportionally to 1 / d, where d is the horizontal distance to the point from the well location 21, i.e. indicated by the intersection of lines 20 and 21. Alternatively, the same effect could be achieved by scaling down each of the values of k. Note that all of these operations are fully reversible. As such, if seismic data input is changed, for example, a new solution of the horizon can be found instantly, based on the same (tuned) parameters, as there is no additional workflow involved. As the above models contain only springs, the only energy in the system is the potential energy stored in the springs. It is thus equivalent to consider the situation from an energy minimization perspective. 2 The potential energy stored in a spring is k Ax The derivative of this with respect to Ax is k Ax, which is the force expression. As such, finding the minimum energy (where the derivative is zero) is the same as finding the horizon with zero net force. At equilibrium, the sum of the forces is zero and the energy is at its minimum; these are fully equivalent ways of looking at the system. Put another way, a method could seek to find a solution where the sum of the forces is zero or where the energy is at its minimum as these are in effect the same thing, i.e. the same solution (horizon) will be found. To find the state of minimum energy of the system, the method can involve determining all possible states or horizons, with all possible cuts or faults, calculating the energy of each of these states, and selecting the state with minimum energy as the best representation of the horizon. Taking this “minimum energy” approach has an advantage: it provides the opportunity to introduce new forms of “energy” that might be difficult to express as a force. An example of this is “breakup energy”. More specifically, with this approach the degrees of freedom may be increased to let the system “decide by itself’ whether the points 6 should be connected by blade springs (i.e. to form a continuous surface 5) or not. With no “breakup energy” in the model, all points would choose not to be connected to anything, as they could then snap directly to the seismic data (to the xs. position) and reach zero potential energy. As this would not be helpful (there would be no smoothing to remove noise), to prevent this from happening, an “energy cost” for breaking up may be introduced. Similar to excitation energy for electrons, breaking up can be chosen to be reversible. If a pair of adjacent points is broken, the system gets more energy but it can release that energy again by reconnecting. By adjusting how much energy is deemed needed to break up a pair of points, the model can be tuned to create small faults everywhere when there is a small amount of strain in the curvature force, to connect across everything, or to reach an intermediate situation with a realistic determination of the location of any faults and smooth surfaces elsewhere. By adjusting the breakup energy, the situation can be adjusted continuously between the two extremes. The chosen value of the breakup energy may be dependent on the data and the user. First, the breakup energy could be set as infinite to have everything connected. The breakup energy could then be reduced until a certain fraction (e.g. 1%) of the bonds are broken. An assessment could then be made about whether there are too few or too many faults and the breakup energy could then be adjusted (increased or decreased) accordingly. In some cases, there may be some a priori information about which points should or should not be connected. For example, if data quality is varying across an area, the breakup energy could be chosen to vary accordingly. More specifically, in any regions where data quality is low, this could be an indication that there are faults in such regions and so the breakup energy could be set as lower for those regions, to permit the creation of a fault in the horizon more easily in such locations. Data quality can be assessed of measure in various ways. For example, data quality can be assessed using an amplitude continuity measurement, doing a cross-correlation between the trace (1D timeseries in time / depth direction) between neighbours, or by more sophisticated methods (e.g. using e.g. seismic tiles). The principles described above could also be used to achieve the same goal as a horizon auto-tracker. An auto-tracker is a method where one or a few points are initially selected, and then a computer picks more points following some rules. Such methods can typically use cross-correlation to find the time location of a point located next to an existing point and use the correlation coefficient between the existing point and the potential new point as a criterion for accepting or not accepting the point to the horizon. However, the methods described herein can be used to perform such tracking, resulting in surfaces that are instantly smooth. This is done by perturbing an existing horizon by adding a new point to it, checking if this results in lower energy and accepting it if the answer is yes. To include the correlation between two traces in the measurement, which is sensible, the breakup energy can be modified by including the cross-correlation coefficient in the term. If the value of the coefficient is close to 1, a higher breakup-energy is obtained and thus the system really wants to include the new point. A threshold can be applied so that the breakup energy “score” becomes zero for a correlation less than a certain amount, e.g. 0.7, resulting in connections never happening. It may be natural to perform the tracking by following the steepest descent down the energy landscape, which would be to include first the points that show strong correlation and that lie on a straight line out from the existing points. If there are areas of poor data quality, so that some points in between show poor correlation and have low amplitude, the system can be perturbed more by more points simultaneously. As it can be favourable to have many connections (fewer faults), it may be better for the system to include a low-amplitude trace in the middle rather than for it to create a cut or gap there. As the “energy” of the system can be determined according to any criteria / function, it may also be decided to e.g. correlate traces to their next-but-one neighbours (x against x x ) and release more energy from the system with a high correlation than a low. This could be way to help include a few scattered bad data points in a horizon. Note that the method described here is a way to find the objectively best horizon interpretation, given the objective function. The results, unlike those for most auto-trackers, will be the same no matter from which point it starts. Thus, the term “tracking” is slightly misleading for this method. There is no “adjust your parameters and your seeds and re-track and then smooth” as is commonly done with known auto-trackers. Rather, the approach is to just change parameters and see what the result is now, and repeat. In some cases, it may be advantageous to release energy if neighbouring grid points have a high correlation also if they do not connect. This would keep this independent of the breakup energy, but possibly only for one such neighbouring grid point per point, as this could prevent the surface from “wandering off” without connecting. A mechanism or rule (such as a geometric criterion related to only having one subtle fault in one location) could be employed to avoid spurious results if tracking is allowed above small faults and not just around them. Once a horizon has been interpreted or defined, according to the methods described herein, data representative of that horizon is stored in a memory. In some cases, the horizons identified by the methods described herein are stored in the form of arrays. This can be advantageous as it allows calculations involving the horizons to be performed quickly, and the horizons can be easily displayed as a map. However, the methods described here are more general and can also be used to describe multi-z surfaces such as thrusts and reverse faults. In such cases, it may be advantageous to represent and / or store the horizons as / with vaex dataframes, for example. A multi-z surface is a surface that has two or more z-values for the same x,y coordinate pair. This situation happens in compressional settings, when a hanging-wall block is pushed over a footwall block. In the methods described herein, the effect of the “curvature force” is to pull each point towards the mean position in a local area. The “basic” expression for the curvature force on the central point is: 1 / 3 Fc .(x. -2x. + x. J i' i-i i i+r If the movement of this force is simulated in a timestep, and a timestep is chosen that is 1 / Fc, the position at the next timestep is then x’. = 1 / 3 Fc. (x.^ -2x. + xM). This is exactly the effect of a running smoother. The system can thus be visualized as an iterative approach of doing a running-sum smoother to smooth the horizon, then “unsmoothing” it a bit by moving each point towards the location of the seismic maximum amplitude, then smoothing again, and so on. This will converge to a solution that is a trade-off of smoothness and “snap to seismic”, and, provided that the edges are handled properly in the running-sum operation, achieve exactly the same results. The tuning parameter in this formulation is how big “a bit” is in the statement above. Following the methods described herein, it is possible to “find” a smooth horizon in noisy data and at the same time cut it at the faults without a priori knowledge of where the faults are. Since the system decides itself whether to add or remove points to / from a horizon, this is effectively an auto-tracker and not just a smoother of existing horizons. The method also provides a way to identify subtle faults in the data. Using 3D data, as would usually be the case, gives fundamentally better results than 2D since the horizons are smooth(er) in 3D whereas the noise is not. It is thus possible to work with higher noise levels in 3D than in 2D. With the methods described herein the smoothness of the horizon can be controlled in a reversible way. The methods can also allow modelling of other subsurface features, such as wells, as well as adding or dropping points to make a horizon bigger or smaller. The methods described herein may provide improved and more accurate interpretation, determination or identification of horizons from seismic data. More accurate knowledge of the location of such horizons can have various applications and benefits. Step (iii) of the method involves using the one or more horizons as interpreted in step (ii) in a model and / or as the basis for a further decision For example, horizons identified from seismic data may be used: - to build models or (3D) representations of reservoirs; in the decision-making process for the drilling of a well, e.g. for deciding, where, how far, and / or in which direction(s) / with which trajectory to drill a well; - to locate a hydrocarbon reservoir more accurately, e.g. in a location 5 where there is already a known / existing oil field; to locate a fault more accurately, e.g. such that better, more informed decisions may be made. Accordingly, the method may then involve drilling a well according to a decision that has been made, e.g. along a particular trajectory, or to a particular 10 point.

Claims

1. A method of horizon interpretation, the method comprising:(i) obtaining seismic data for a subsurface region; and(ii) analysing the seismic data to interpret one or more horizons within that data;wherein step (ii) comprises modelling the horizon as a surface with at least two different kinds of forces acting on the surface, and determining the horizon for which:the net forces acting on the surface are zero, and / or an energy function representing the horizon is minimised.

2. A method as claimed in claim 1, wherein the surface comprises an array of points, each point corresponding to a location in the seismic data.

3. A method as claimed in claim 1 or 2, wherein the at least two different kinds of forces acting on the surface comprise a first kind of force and a second kind of force, the first kind of force being an internal force to the surface and the second kind of force being an external force to the surface, wherein the first kind of force is preferably representing a force that acts to flatten the surface and / or the second kind of force is preferably representing a force that acts to pull the surface towards maxima in the seismic data.

4. A method as claimed in any preceding claim, wherein the surface is modelled as a blade spring.

5. A method as claimed in any preceding claim, wherein at least one of the at least two different kinds of forces acting on the surface are modelled as linear spring forces.

6. A method as claimed in any preceding claim, further comprising tuning a magnitude of at least one of the at least two forces acting on the surface.

7. A method as claimed in any preceding claim, the method further comprising determining whether the inclusion of any cuts or gaps in the horizon result ina reduction in an energy function representing the horizon and, if so, which cuts of gaps in the horizon minimise the energy function.

8. A method as claimed in claim 7, wherein the energy function includes an energy representing an energy cost when a cut is made in the horizon.

9. A method as claimed in claim 8, further comprising tuning a magnitude of the energy representing an energy cost when a cut is made in the horizon.

10. A method as claimed in any preceding claim, wherein step (ii) comprises modelling the horizon as surface with at least three different kinds of forces acting on the surface, wherein at least one kind of force acting on the surface represents a force pulling the surface towards a known object in the subsurface.

11. A method of modelling a subsurface region, the method comprising performing the method of any preceding claim to interpret at least one horizon from seismic data, and using the interpreted horizon as at least one input for modelling the subsurface region.

12. A method of determining a trajectory or location of a well, the method comprising performing the method of any preceding claim to interpret at least one horizon from seismic data and / or model a subsurface region, and using the interpreted horizon and / or modelled subsurface region to determine the determine the trajectory or location of the well.

13. A method of drilling a well, the method comprising performing the method of claim 12 to determine the trajectory or location of the well and then drilling the well with that trajectory or location.

14. A computer program product comprising computer readable instructions that, when executed by one or more computers, cause the one of more computers to perform the method of any of claims 1 to 12.

15. A system comprising one or more computers and one or more memories, the system being configured to perform the method of any of claims 1 to 12.

Images