Geophysical Data Processing Systems

FIELD OF THE INVENTION

This invention relates to methods, apparatus, and computer program code for processing geophysical data, more particularly gravity gradiometer data from a potential field survey to provide a representation of the underlying geology of the surveyed region. Embodiments of the techniques we describe are particularly useful for processing data from airborne surveys, in particular gravity field surveys.

BACKGROUND TO THE INVENTION

A potential field survey is performed by measuring potential field data which, for a gravity survey, may comprise one or more of gravimeter data (measuring gravity field) or gravity gradiometer data (measuring gravity field gradient), vector magnetometer data, true magnetic gradiometer data, and other types of data well-known to those skilled in the art. A common aim of a geophysical potential field survey is to search for signatures which potentially indicate valuable mineral deposits.

The goal of most potential field surveys is to deduce the sub-surface geology. It is well known to the geophysics industry that the high spatial frequency potential field signals are dominated by near surface anomalies. Further, in most situations, the terrain and the near surface geology also give rise to substantial low spatial frequency signals. In general, the low frequency signals from the near surface anomalies are extremely hard to distinguish from the low frequency signals coming from deeper geology. Accordingly, it is possible to postulate a geology having only a near surface mass distribution which matches the measured signal implying that there are no anomalies left in the measured signal coming from deep geologies. Fortunately most such mass distributions can be ruled out using a range of tools including geological knowledge,

known typical density distributions, prior information, similar known structures elsewhere in the world etc.

There remains a need, however, for improved techniques for processing geophysical data from such surveys in order to identify the underlying geology.

SUMMARY OF THE INVENTION

According to a first aspect of the invention there is therefore provided a method of processing geophysical data including at least measured potential field data from a potential field survey of a surveyed region of the earth to provide a representation of the geology of said surveyed region, the method comprising: generating a first model of said surveyed region by fitting data predicted by said first model to said measured data for a specified frequency range; predicting full range potential field data for all measured frequencies using said generated first model; comparing said full range predicted data to said measured potential field data to provide full range residual data representing a difference between the full range predicted data and the full range measured data, and interpreting said full range residual data to provide a representation of said geology of said surveyed region.

hi general, the goal when processing geophysical data is to try and find structure below near surface geology. Accordingly, an Earth Model at least down to and including the depth of interest needs to be generated. The signals from shallow structures are relatively larger than those from the deeper sources, thus it is particularly important to be able to construct a good Shallow Earth Model (SEM). A SEM usually consists of either or both of two elements, the topography (which can be assumed to be known to a high accuracy) and the depth to bedrock (which can only be inferred). The depth of an SEM depends on the geology of the surveyed region and may be a few hundred metres to several kilometres. The SEM maybe in the range of 100m to lkm.

The high spatial frequencies of potential field data are dominated by the near surface sources, whereas the low frequency signals are generated by both the shallow and deep sources. The low frequency signals are thus not unique to the near surface features. Accordingly, the method preferably comprises generating said first model for high frequencies wherein said first model may be a shallow earth model. The specified frequency range depends on the nature of the surveyed region and may be determined by the parameters of the first model. For example, if the first model is a SEM and the SEM has a depth of approximately 1 km, the specified high frequency range represents frequencies having a wavelength of less than 1 km.

Said first model may be generated by comparing the data predicted by said first model and the measured data to provide initial residual data, filtering said initial residual data by frequency and minimising said filtered initial residual data. Fitting parameters (e.g. density) of the first model to the high frequency data will provide a more accurate shallow earth model (SEM) which is not contaminated by any low frequency signals. The idea of fitting parameters to a pre-defined bandwidth of a given data set is known as a focused inversion. If the method attempted to fit the SEM to the full bandwidth of the signal, the SEM parameters would be of inferior quality due to their attempt to represent the low frequency content in the signal that originates from deeper sources.

According to another aspect of the invention there is provided a method of processing geophysical data including at least measured potential field data from a potential field survey of a surveyed region of the earth to provide a representation of the geology of said surveyed region, the method comprising: generating a first model of said geology using a focussed inversion; predicting full range potential field data using said generated first model; comparing said full range predicted data to said measured potential field data to provide full range residual data representing a difference between the predicted data and the measured data, and interpreting said full range residual data to provide a representation of said geology of said surveyed region.

For both aspects, said full range residual data may be interpreted using standard techniques. For example, the interpreting step may comprise generating a second model of said geology by fitting data predicted by said second model to said full range residual data, and generating a representation of said geology of said surveyed region using said second model. Alternatively, a full second model may not be constructed and said full range residual data may be interpreted to identify key geological features.

Said second model may be calculated by minimising the differences between the full range residual data and data predicted by the model given a set of parameters, e.g. by using well known fit measures such as least squares. Where the first model is a shallow earth model, the second model may be a model of the deeper or underlying geology of the surveyed region.

For both the first and second models, if the geology being modelled consists of regions of unknown but constant densities, the model parameters are these densities. The prediction function pred(x) may be stated in matrix form as

pred(x) = Ax

where A is a matrix of constant numerical values and x is the set of densities.

hi both aspects, the potential field data preferably comprises gravity gradiometer data, particularly for deriving the SEM. In its simplest form a gradiometer is two accelerometers separated by some distance. There will be a very small change in gravity between the two sensors and it is these differences which are measured by all gravity gradiometers. However, the initial accelerations that contaminate conventional gravity are common to both sensors and are therefore not seen for a properly balanced gravity gradiometer. This fundamental characteristic means that gravity gradiometry is the ideal choice for high dynamic environments and allows for a much higher bandwidth measurement to be recorded.

Gradiometry has superior performance over almost all frequencies usually encountered in a geophysical survey and most particularly at higher frequencies (e.g. for

wavelengths less that approximately 20000 metres). Accordingly, gradiometry is particularly useful for deriving a SEM.

With current technology, the crossover at which other measurement systems perform as well as, or better than, gradiometry appears to be at 23km wavelength. Beyond 34km wavelength a gravity gradiometer still recovers long wavelength signals, albeit at a reduced accuracy compared with conventional gravity. At these longer wavelengths it is difficult to quantify the smallest signal present in a gradiometer survey. Nevertheless, it is possible to derive a model for the sub-surface geology using gradiometry, although other potential field data may be used to derive this model.

According to another aspect of the invention, there is provided a method of processing geophysical data including at least measured gravity gradiometer data from a potential field survey of a surveyed region of the earth to provide a representation of the underlying geology of said surveyed region, the method comprising: processing said gravity gradiometer data to generate a shallow earth model, removing data predicted by said shallow earth model from said measured gravity gradiometer data to leave residual data; processing said residual data to generate a model of sub-surface geology, and generating a representation of said underlying geology of said surveyed region using at least one of said models.

The method may further comprise generating the shallow earth model by fitting data predicted by said shallow earth model to said measured data for high frequencies. The data predicted by said shallow earth model may comprise predicted data across the whole bandwidth of the measured data. As explained above, this ensures that the SEM is an improvement over known SEMs.

hi each aspect, the method may comprise performing terrain correction by subtracting- off particularly the higher frequencies (which dominate the power spectrum). A signal with wavelength λ falls off with height z as exp(-kz) where k = 27rλ(from which it can be seen that longer wavelengths are less attenuated) and the wavelength scale corresponds to a signature expected given a target's size and depth. This technique is

described in more detail in the applicant's pending UK patent application GB0715031.1 which is incorporated by reference in its entirety.

In embodiments, in particular where the field comprises a gravity field, the processing to generate a model comprises identifying spatial features by identifying one or more of maxima, minima and points/lines of inflection in the filtered or residual potential field data. Thus, for example, with the gravity gradient tensor the on-diagonal components G _xx , G _yy , and G _zz , (which are differential signals) are processed to determine inflection points or changes in slope since these generally correspond to geologically significant features of the surveyed region such as an interface between two different types/densities of rock. For off-diagonal components, in particular G _2x and G _zy (which emphasise symmetries in the x- and y- direction respectively, maxima and/or minima are preferably identified; for G _xy points are preferably identified by locating pairs of dipoles since these tend to identify corners of a subterranean body. Off-diagonal elements G _Z1 (where i is x or y) tend to emphasise symmetries in the i-direction. On- diagonal components G _x* and G _yy are always zero along respective axis x and y of the source and since the choice of axis is often arbitrary, optionally the co-ordinates system may be rotated about one or more axis to potentially identify further geologically useful information. In embodiments the co-ordinate axis may be selected to maximise the apparent useful geological information. Similarly magnetic data may be processed to identify maximum/minimum inflection points/lines.

Referring again to the above-mentioned maxima, minima and lines of inflection, multiple sets of spatial features maybe identified for the filtered or residual potential field data, for example from different vector or tensor components of a surveyed gravity field, magnetic field or gravity gradient field and/or from survey data other than from a potential field survey. Other survey data which may be combined with the data obtained from a potential field survey includes (but is not limited to): topographic information, for example determined by lidar, spectral or more preferably hyperspectral imagery, gas saturation data, chemical analysis data (from soil sampling), and other soil survey data.

In preferred embodiments the potential field survey is conducted from a moving platform such as an aircraft.

The method may further comprise generating fault polygon data from the 3D representation of the underlying geology of the surveyed region. As the skilled person will understand, fault polygon data comprises data representing a location of one or more geological faults on a surface or horizon of the surveyed region. Such a fault polygon may be defined by a plurality of corner points and/or edges, in particular defining a loop, for example on the earth's surface. Such a fault polygon may thus define horizons of stratigraphic layers. Fault polygons are useful because it is generally undesirable to drill at the location of a fault because typically the geology changes to either side of a fault and thus by drilling through the fault one may miss the desired oil or mineral.

The invention also provides a method of prospecting for oil or a mineral from the earth, the method including conducting a potential field survey according to an aspect or embodiment of the invention as described above to generate a representation of the underlying geology of the surveyed region, and then using this representation to locate the desired oil or mineral. A further aspect of the invention also provides oil or mineral extracted using this technique.

The invention further provides processor control code to implement the above-described methods, in particular on a data carrier such as a disk, CD- or DVD-ROM, programmed memory such as read-only memory (Firmware), or on a data carrier such as an optical or electrical signal carrier. Code (and/or data) to implement embodiments of the invention may comprise source, object or executable code in a conventional programming language (interpreted or compiled) such as C, or assembly code, code for setting up or controlling an ASIC (Application Specific Integrated Circuit) or FPGA (Field Programmable Gate Array), or code for a hardware description language such as Verilog (Trade Mark) or VHDL (Very high speed integrated circuit Hardware Description Language). As the skilled person will appreciate such code and/or data may be distributed between a plurality of coupled components in communication with one another, for example distributed across a network.

The invention further provides a data processing system configured to implement embodiments of the above-described methods, to determine one or more parameters relating to physical properties of the Earth's interior from processed geophysical data. Such a data processing system may comprise: data memory for storing measured potential field data and plot data for representing the underlying geology of the surveyed region, program memory storing processor control code as described above; and a processor coupled to said data memory and to said program memory to load and implement said control code.

BRIEF DESCRIPTION OF THE DRAWINGS

These and other aspects of the invention will now be further described, by way of example only, with reference to the accompanying figures in which:

Figure 1 shows an aircraft with flight survey data, and an example of a data processing system configured to implement an embodiment of a method according to the invention;

Figure 2 shows a flow diagram of a procedure for processing measured potential field data to implement an embodiment of a method according to the invention;

Figure 3 is a map of the simplified geology of an area being surveyed;

Figure 4 is a photograph of an example of terrain encountered when conducting a survey;

Figure 5 shows the densities and stratigraphy for the area of Figure 3;

Figure 6 is a map of the surface geology of the area of Figure 3, the map created according to the invention;

Figure 7 is an example of an interface depth constraint added to the model used to create the map shown in Figure 6;

Figure 8 is the resultant shallow earth model (SEM) of the area of Figure 3 together with a balanced cross-section and section through the (SEM);

Figure 9 shows the surface location and geology for a well, and

Figure 10 shows the vector G _zz highlighting the structural variations at the Debolt level for the sub-area shown on Figure 9.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

Potential field surveys

When we refer to a field, in particular a gravity field, this is not limited to a vector field but includes scalar and tensor fields, a potential field and any derivatives deriving from the potential field. Potential field data includes, but is not limited to, gravimeter data, gravity gradiometer data, vector magnetometer data and true magnetic gradiometer data. Elements and representations of a potential field maybe derived from a scalar quantity.

Consider an airborne potential field survey such as a gravity survey, flown on a grid pattern defined by orthogonal sets of parallel lines (flight paths) on a two-dimensional surface which is draped over the underlying terrain. When looking for underlying anomalies the nearby mass has a dominating effect and to provide an accurate representation of deep features a good representation of surface features is desirable so as to be able to perform terrain correction by subtracting-off particularly the higher frequencies (which dominate the power spectrum). A signal with wavelength λ falls off with height z as eχp(-kz) where k = 2ττ/λ (from which it can be seen that longer wavelengths are less attenuated) and the wavelength scale corresponds to a signature expected given a target's size and depth.

For gravity, the relevant potential is the gravity scalar potential, φ(r), defined as

where r, p{r'), G are respectively, the position of measurement of the gravity field, the mass density at location r' , and the gravitational constant. The gravitational acceleration, which is how a gravitational field is experienced, is the spatial derivative of the scalar potential. Gravity is a vector in that it has directionality. It is represented by three components with respect to any chosen Cartesian coordinate system as:

Each of these three components varies in each of the three directions and the nine quantities so generated form the Gravity gradient tensor:

G =

Although there appear to be nine components of the gravity gradient tensor, there are only five independent components. Firstly, the tensor is symmetric as the order of differentiation of a scalar quantity does not matter (implying that G _xy = G _yx ). Secondly, the sum of the diagonal terms equals zero (Laplace's equation). The ability to measure five spatially independent gravity components has obvious advantages over conventional gravity measurements, which only recovers the vertical component (G _z )

There is a relationship between the depth (and shape) of a buried object and the wavelength (and amplitude) of the detected signal. In general, a measured quantity - say a component of the gravity vector or of the gravity gradient tensor will be a

summation of the form shown below. Here we use gg as notation for the measured quantity, for example G _zz .

SS calculated \ measure) 2-i mass-element \ measure ' mass-element ) all-masses

hi the above equation F is called a Greens function (see for example, RJ. Blakely, "Potential Theory in Gravity and Magnetic Applications", Cambridge University Press, 1995, at page 185, incorporated by reference) and r _mass . _e ι _ement defines the location of the mass element (for example the centre of gravity or some other defined point).

The functions F are standard functions, essentially, the influence a source (mass element) of unity mass or density and defined shape would have at the relevant (measurement) point. The source may be a point source, sphere or ellipsoid but, in practice is more often a prism, which may be irregular. For example, if the presence of a particular geological layer or, say, geological anomaly, e.g. a kimberlite pipe, is suspected a shape can be defined to take account of this. A number of textbooks list Greens functions for simple shapes; functions for more complex source geometries can be found in the literature. Also the source influence superposes so that if a complex shape can be discretised into a plurality of simpler shapes then the Greens functions for the discrete shapes can be added together. This in principle allows numerical values for the Greens function of any arbitrary shape to be determined, although in practice relatively simple shapes are generally preferable. By way of example, the Green's function F for a rectangular prism (Blakely, ibid, at page 187), has 8 terms each of which corresponds to a vertex of the prism.

Referring now to Figure 1, this shows an example of an aircraft 10 for conducting a potential field survey to obtain data for processing in accordance with a method as described above. The aircraft 10 comprises an inertial platform 12 on which is mounted a gravity gradiometer 14 which provides potential field survey data to a data collection system 16. The inertial platform 14 is fitted with an inertial measurement unit (IMU) 18 which also provides data to data collection system 16 typically comprising attitude data (for example, pitch, roll and yaw data), angular rate and angular acceleration data,

and aircraft acceleration data. The aircraft is also equipped with a differential GPS system 20 and a LIDAR system 22 or similar to provide data on the height of the aircraft above the underlying terrain. The aircraft 10 may also be equipped with other instrumentation 24 such as a magnetometer, TDEM system and/or hyperspectral imaging system, again feeding into the data collection system. The data collection system 16 also has an input from general aircraft instrumentation 26 which may comprise, for example, an altimeter, air and/or ground speed data and the like. The data collection system 16 may provide some initial data pre-processing, for example to correct the LIDAR data for aircraft motion and/or to combine data from the EVIU 18 and DGPS 20. The data collection system 16 may be provided with a communications link 16a and/or non-volatile storage 16b to enable the collected potential field and position data to be stored for later processing. A network interface (not shown) may also be provided.

Data processing to generate map data for the potential field survey is generally (but not necessarily) carried out offline, sometimes in a different country to that where the survey data was collected. As illustrated a data processing system 50 comprises a processor 52 coupled to code and data memory 54, an input/output system 56 (for example comprising interfaces for a network and/or storage media and/or other communications), and to a user interface 58 for example comprising a keyboard and/or mouse. The code and/or data stored in memory 54 may be provided on a removable storage medium 60. In operation the data includes data collected from the potential field survey and the code comprises code to process this data to generate map data, for example in accordance with the procedure shown in Figure 2, described below.

Referring now to Figure 2 this shows an example of a procedure for implementing on a data processor which may, in embodiments, comprise a general purpose computer system, for processing data from a flight survey in accordance with the previously described techniques. Thus, at step S200 the procedure inputs the measured potential field data, for example gravity gradiometer, and associated 3D position data. Optionally at step S200a, some pre-processing may be applied, for example to remove anomalies and/or to decrease (or increase) or select the data to be processed.

At step S202 the procedure fits a model to the high frequency content. Because the high spatial frequencies of the signal are dominated by the near surface sources, fitting parameters (e.g. densities) of a model to these frequencies will provide for a shallow earth model (SEM). The idea of fitting parameters to a pre-defined bandwidth of a given data set is known as a focused inversion and is achieved by applying:

m.m\L[Filter(pred(x) — m)] }

where Filter is a spatial filter tailored to remove frequencies outside the range of interest, i.e. low frequencies in this example, L is a fit measure (e.g. standard least squares Z ₂ norm), pred(x) are the predicted signals from the model given the set of parameters x, and m represents the set of measurements.

Standard algorithmic methods can be applied to adjust the model parameters x in such a way to achieve the minimum. The solution x will be that which fits the high frequencies in the measurements m and is thus a more accurate solution for the near-surface geology.

Optionally terrain correction data may be included in the SEM as indicated at stage (S204).

After solving for the unknown model parameters in a focused inversion, at the next stage (S206), the entire bandwidth of the signal due to the model is calculated and then removed from the measurements at step S208 by applying,

residual = m - pred(x) .

The residual signal represents only the signals from the deeper sources and becomes the preferred data set for interpreting the deeper unknown geology. Thus at step 210, a model for the sub-surface geology may be generated, e.g. by applying the standard method of minimising the residual of the model predictions to the measurements, for example by solving

mm{L[pred(y) - residual] } ,

where L is a fit measure (e.g. standard least squares L ₂ norm),pred(y) are the predicted signals from the model given the set of parameters y, and residual represents the residual set of measurements after the predicted signal from the SEM is subtracted.

A forward calculation can be performed using these parameters to generate a desired map of the potential field for the sub-surface geology.

In both prediction steps 202 and 210, it may be possible to simplify the calculation of the model. For example, if the near surface area to be modelled consists of regions of unknown but constant densities, the model parameters x are these densities. In this case, the inversion simplifies to a linear case where the prediction function pred(x) can be stated in matrix form as

pred(x) = Ax

where A is a matrix of constant numerical values.

Such a formulation would be useful in determining the unknown densities in a terrain model since the geometry of the model is often accurately known from lidar surveys.

For sub-surface geological models, the geometry can be inferred from other sources; well data, seismic surveys, etc. Accordingly, the model parameters y maybe unknown but constant densities as for the near surface model.

Figure 3 shows a geological map of an area being surveyed. The thickest and most stratigraphically complete part of the carbonate/elastics sequence occurs in the west of the region. Figure 3 shows well data indicating oil and gas. High exploration within the eastern margin is evidenced by the density of well data. There is a sharp western boundary to the loci of well data plotted with only a few wells drilled west of this

margin. This locality corresponds to the eastern limit of cordilleran deformation and an increase in topographic relief.

Figure 4 is a photograph of the terrain encountered in this region which shows the outcropping geology. Crt is the Cretaceous elastics, Utrs the Upper Triassic Carbonates and LTrs the Lower Triassic Carbonates. Clearly the terrain is mountainous and thus represents a difficult challenge for many traditional processing methods. With various lithologies outcropping at the surface (including high density carbonates), a simple tri- model density distribution would not correct across all the area sufficiently and a detailed shallow earth model is required to take into account all the variations.

A gradiometer measures the rate of change of density within the Earth's subsurface and the key density interfaces giving rise to the gradient signal in the area are the Cretaceous, Upper Triassic and Lower Triassic units. Although there may be subtle variations within these lithologies, the system effectively measures the aggregate density of each layer. Figure 5 shows the strata and associated densities for the area. The main exploration target for the area is the Debolt level Mississipian carbonates at a depth of 2km subsea plus 1.5km of terrain. In order to image the Carboniferous Debolt level structures, the SEM would have to accurately model all of the upper strata down to a detachment surface within the PermoTriassic elastics. Thus, for this example, an SEM is 3.5km thick whereas in other geologies, an SEM may only be a few hundred or even a few tens of metres thick.

In order to build a high resolution SEM, an accurate representation of the terrain is necessary. A detailed digital terrain model (DTM) is first generated using LiDAR data. The next stage is to obtain the surface geology map to determine the outcropping units and map them onto the DTM. In certain areas, the surface geology maps were found to be inaccurate and magnetic gradiometry data is used to remap the carbonate contacts since slightly ferrous units above and below the carbonates produce a weak magnetic signal at surface outcrop.

The resultant surface does not take into account the depth extent of the various layers and the map is enhanced by building thirty geological cross-sections for each key

interface spanning across the entire survey area, e.g. as shown in Figure 7. The resultant cross-sections provide a series of depth constraints for each interface which are integrated into the newly mapped DTM as shown in Figure 6.

With all the information gathered, the next stage of building the SEM is to derive the resultant surfaces prior to removal from the full earth model. A simple gridding operation applied to the surfaces would not correctly honour the topographic variations and, by using GIS (geographic information systems) analysis methods, the LiDAR data is used to generate dip and plunge information. The resultant surfaces are then assigned their respective density variations and the gradient response of each key interface forward calculated and removed form the total Earth model. The remaining density distribution following subtraction would now represent the density variations arising from the Debolt level. Figure 8 shows the resultant SEM, with each layer clearly displayed.

Figures 9 and 10 show how the new methodology described above could add additional understanding to previously explored targets. As shown in Figure 9, a well was drilled based on the surface anticline and 2D seismic data that appeared to image the main structure. Despite chasing borehole dip data to sidetrack towards the anticline axial plane, the well was too far off structure and therefore abandoned. Figure 10 shows the sub-surface map generated according to the invention, i.e. with the SEM data removed. Figure 10 shows that the actual structural target is positioned east of the surface anticline. This is because the Debolt anticline sits beneath a detachment surface located within younger Triassic strata. Consequently, the seismic drilling line and placement for the well only clipped the northern limit of a wrench- fault bounded anticline. Applying the methodology described above, the well location would probably have been changed.

No doubt many other effective alternatives will occur to the skilled person. It will be understood that the invention is not limited to the described embodiments and encompasses modifications apparent to those skilled in the art lying within the spirit and scope of the claims appended hereto.