GB2428827A - Gravity survey data processing - Google Patents

Abstract

A method of processing potential field measurement data, such as from an airborne gravity survey, in which a plurality of relationships between the measurement and the positions, in three dimensions, associated with each measurement are determined, each relationship relates a measurement to a function of the position, preferably to estimate the coefficients of a harmonic expansion of the function.

Description

I

Gravity Survey Data Processing This invention relates to improved techniques for processing potential field measurement data from airborne surveys such as gravity surveys, and to improved techniques for data acquisition which are enabled by the improved data processing techniques.

Conventionally airborne potential field surveys such as gravity surveys are flown on a grid pattern. The grid is defined by orthogonal sets of parallel lines (flight paths) on a two-dimensional surface which is draped over the underlying terrain. The draped surface satisfies a minimum height constraint (defined by the closest the aircraft is permitted to fly to the ground) and by a constraint on. the maximum rate of climb/descent of the aircraft, typically around three percent. This approach suffices for flat telTarn but for hilly or mountainous terrain the surface on which the aircraft flies can vary by as much as two or three kilometres from, say, the bottom of an underlying valley to the top of the mountains/survey area so another approach is needed.

ft is useful to be able to collect potential field data, in particular gravity data, from close to the ground ie. a low height. In a gravity survey nearby mass provides both high and low (spatial) frequency data, whereas the influence of a deeper mass is seen primarily at lower frequencies only. When looking for underlying anomalies the intervening 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 subtract-off particularly the higher frequencies (which dominate the power spectrum). For example, a signal with wavelength ? falls off with height z as exp('-kz.) where k = 2irJ from which it can be estimated that a signal component of wavelength 200 metres from a mass at a depth of 100 metre has fallen to approximately 1/20 of its initial value at the earth's surface (and is progressively further attenuated with increasing height), whereas it can be seen that longer wavelengths are much less attenuated. Generally the size and position of a survey is broadly chosen according to a wavelength scale corresponding to a signature expected given the target's size and depth.

It will be appreciated from the foregoing discussion that it is generally desirable to be able to perform a flight survey at a low height, but in practice the aircraft's limitations and gridded flight plan can impose significant restrictions. The gridded flight plan is necessary because conventional techniques for processing gravity survey data rely on a constant height assumption. Broadly, the assumption is that for a deep source the aircraft's height is approximately constant, the error from this assumption supposedly being only a small correction. Furthermore conventional gravity survey data processing techniques rely on regularly spaced data points, normally a power of two, so that a (fast) Fourier transform can be applied, this defining the requirement for orthogonal sets of parallel flight paths. The flight paths are required to be in a common surface since existing techniques assume, for example, that where two paths cross they cross at the same height. A further problem with existing flight paths arises when the surveyed area is not precisely rectangular, for example because of the local terrain. In this case in order to be able to apply the conventional techniques the data points are "padded", for example by interpolation or extrapolation to generate a regular set of data points over a rectangular area. The wavelengths (or more particularly, wavenumbers) used for the Fourier analysis are then determined by the maximum x andy (length and breadth) dimensions of the now padded rectangular area.

In view of the above drawbacks of conventional techniques, improved potential field survey data processing techniques and survey flight patterns are desirable According to a first aspect of the present invention there is therefore provided a method of processing measured potential field data from a potential field survey of the earth to determine map data for mapping said field, the method comprising: inputting said measured potential field data, said measured potential field data comprising data defining a plurality of potential field measurements and associated positions, each said position defining a position of a said potential field measurement in three dimensions; determining a plurality relationships between said potential field measurements and said positions, each said relationship relating a said potential field measurement to a function of a said associated position in three dimensions multiplied by a field mapping parameter; and determining a substantially self-consistent set of said field mapping parameters for said plurality of relationships to thereby determine said map data.

The potential field data may comprise gravimeter data (measuring gravity field), gravity gradiometer data (measuring magnetic field gradient), vector magnetometer data, true magnetic gradiometer data, or other similar data. The field mapping parameters (coefficients or parameters which allow the characterisation of the potential field) may be used to generate similar potential field data and in embodiments, broadly speaking, the parameters are determined by selecting parameter values which define potential field data such as a gravity field or gravity gradient field which most closely approximates the measured potential field data. This may comprise, for example, rninimising a mean squared error or another similar optimisation.

In embodiments the relationships comprise a set of simultaneous equations which may be represented, for example, by a matrix, and the determining comprises solving these equations. Preferably the equations are overconstrained and the solving comprises solving by a noise-reducing procedure such as a least squares optimisation. Preferably the relationships are such that they have the form of a real (gravity) potential field so that noise, which will not generally have this form, is in effect discarded by the optimisation procedure. Thus, preferably, the (three-dimensional) position function satisfies Laplace's Equation. In some representative examples of the method there may be of order 10K parameters and TOOK positions.

The map data may simply comprise the set of field mapping parameters (the potential field coefficients or parameters) since this data may be used to generate a range of potential field data (for example gravity field or gravity gradient field data) depending upon the desired height, region of coverage and/or surface to be mapped. However the method may further comprise determining the relevant field, typically on a surface, from the set of field mapping parameters to generate a map, that is performing a forward calculation of potential field data from the set of field map parameters. As will be mentioned again later, embodiments of the method generate the same type of potential field data as that measured, at least implicitly in determining the field mapping parameters through a procedure which aims to minimize or otherwise optimise differences between calculated and measured values. That this is valid whether potential field data comprises, say, gravimeter data or gravity gradiometer data, can be seen by application of one of Green's Theorems. Furthermore knowledge of one type of potential field data such as a gravity field or gravity gradient field, allows the other to be determined by integration/differentation, as well as the scalar potential.

In one embodiment of the method an equivalent source representation of the potential field data is employed. In this embodiment a surface density or (mass) values for a plurality of surface mass elements are found which together produce a gravity (or other) field which best approximates the measured gravity (or other) field. The surface is generally, but need not be planar and, for example, may be chosen to approximate the surface of the earth in the surveyed region. Once values for these mass elements has been found a simple forward calculation can be employed to predict the scalar potential or gravity field or gravity gradient field, for example to generate any desired map of the surveyed region. Typically a forward calculation is made to determine one or more

components of the field over a flat mapping plane.

In outline, in embodiments of this method the square of the difference between the measured field and that determined from the equivalent source masses (whether gravity or gravity gradient is measured/calculated) is a function of the mass errors. Thus by partial differentation of this function with respect to, separately, the mass of each surface element a set of simultaneous equations is formed in which the only variables are the masses of each surface mass element (providing the measurement position with respect to the surface mass element, r- ,, is known). Conveniently this set of simultaneous equations can be represented by a matrix, and solved for values of the surface mass elements. Since the number of measurements is generally much greater than the number of simultaneous equations (i.e., the number of mass elements), for example at least five or ten times greater, the set of equations is overconstrained. This is useful because if this were not the case the values of the surface mass elements would ft to the noise; instead an optimisation procedure, for example a least squares fit, is used to determine the values of a surface mass elements, thus reducing the noise. Or moreover because a physical model gravitational field is used the equations obey Laplace's Equation so that it is the gravitational field rather than the noise which tends to be fitted.

In a second embodiment of the method a modified Fourier basis functions technique is employed. In this method a two-dimensional Fourier expansion is used to represent the measured data but with a height-dependent coefficient (e') multiplying each two- dimensional component (eH) this representation is still completely general and obeys Laplace's Equation (when differentiated the multipliers of the height-dependent coefficient and of the 2D components cancel). The expansion can be truncated at a point dependent upon the desired accuracy of the representation, for example determined by the number of dependents made. Other expansions can be employed, for example spherical harmonics, which are useftil for larger areas, in particular where the curvature of the earth becomes significant.

Broadly speaking coefficients for the expansion are determined in the following way.

At position r the measured potential field (gravity or gravity gradient) is known and, for a particular pair of values k,,, and k, the value of e'e''e' is also known, Thus in the Fourier expansion (for example, refer to Equation (3) later) every data point can, in effect, be treated separately and, as with the first embodiment described above, a set of simultaneous equations is generated which can be solved for the Fourier expansion coefficients. These coefficients can then be employed to determine (by forward calculation) the gravity or gravity gradient on any two-dimensional surface, that is to provide a map. It will be recognised that with this (and the preceding) embodiment the data can be randomly sampled since the data points are effectively disconnected - all that is needed is a set of measurements and associated x, y, z,) positions in order to be able to determine the gravity or gravity gradient field on any other set of points, in particular on a surface or plane, that is a map. Furthermore it can be seen from this that, by conventional methods, there is no need for padding of measured data to build up a rectangle. Conceptually the data points can be treated as an arbitrarily shaped portion of a two-dimensional region of maximum dimensions determined by maximum wavelengths for the Fourier expansion in the x and y directions (although, as can be seen, there is no need for the measured data to lie on a two-dimensional surface). It can be seen, therefore, that values for k,, and k can be chosen independent of the survey dimensions (here k being wavenumber, that is 2ir/? whereas with a conventional survey the maximum wavelength must correspond to either the length or breadth of the rectangular surveyed area).

In another aspect the invention provides a data processing system for processing measured potential field data from a potential field survey of the earth to determine map data for mapping said field, the system comprising: data memory for storing said measured potential field data, said measured potential field data comprising data defining a plurality of potential field measurements and associated positions, each said position defining a position of a said potential field measurement in three dimensions; program memory storing processor control code; and a processor coupled to said data memory and to said program memory, to load and implement said control code, said code comprising code for controlling the processor to: input said measured potential field data; determine a plurality relationships between said potential field measurements and said positions, each said relationship relating a said potential field measurement to a function of a said associated position in three dimensions multiplied by a field mapping parameter; and determine a substantially self-consistent set of said field mapping parameters for said plurality of relationships to thereby determine said map data.

In a further related aspect the invention provides a method of processing measurement data from an airborne gravity survey to provide data for a gravity map, the measurement data comprising a plurality of gravitational potential field measurements each with an associated measurement position in three dimensions, the method comprising using said measurement data to estimate coefficients multiplying a function of said three- dimensional position in a harmonic expansion to thereby determine a representation of a gravity field in an area of said survey for said gravity map.

The invention further provides a method of processing measurement data from an airborne gravity survey to provide data for a gravity map, the measurement data comprising a plurality of gravitational potential field measurements each with an associated measurement position in three dimensions, the method comprising using said measurement data to estimate a plurality of mass elements each multiplying a function of said threedimensional position in an equivalent source representation of said gravity field to thereby provide said data for said gravity map.

The invention further provides a method of mapping a gravity field, the method comprising: inputting measurement data from an airborne gravity survey, the measurement data comprising a plurality of gravitational potential field measurements each with an associated measurement position in three dimensions; determining a plurality relationships between said potential field measurements and said positions, each said relationship relating a said potential field measurement to a function of a said associated position in three dimensions multiplied by a field mapping parameter; and mapping said gravity field by determining a substantially self-consistent set of said field mapping parameters for said plurality of relationships.

Embodiments of the above described methods further comprise surveying the field to obtain the measured potential field data by flying a set of paths which are not constrained to be parallel or define a rectangular grid pattern.

As previously mentioned, embodiments of the above method may be employed, if desired, with effectively randomly sampled data. Whilst random sampling may not be practically convenient nevertheless it will be appreciated that because the methods we describe are generally unconstrained the flight survey paths may be chosen so as to obtain better initial data, in particular by flying as low as possible constrained only by general safety considerations and the aircraft's capabilities. Thus, for example, a broadly conventional set of paths may be supplemented by additionally flying along the lines of valleys and/or modified by, say, curving one or more flight paths to fit around a mountain.

Although there is no need for paths to cross at substantially the same height (a "crossover"), nevertheless crossovers are useful for reduction of noise and low- frequency drift because they provide two closely correlated measures of similar regions.

However, unlike conventional surveys, there is no need for the paths to cross at the same height. Furthermore a variety of crossing patterns is also available, for example in a valley two piecewise linear or serpentine paths could zig-zag in opposite directions to provide a succession of crossovers along the paths.

According a further aspect of the invention there is therefore provided a method of flying an airborne potential field survey, the method comprising flying an aircraft along a set of paths and measuring potential field data at points on said paths, wherein said set of paths has one or more of the following features: two paths cross at heights differing by at least 50 metres; in a region of said survey paths in the same general direction are non-parallel by greater than 5 degrees; said paths include curved paths; said paths of said set of paths, taken together, do not substantially lie in a surface; said paths of said set of paths, taken together, define a surface wherein at least one of said paths defines one of two orthogonal directions in said surface such that said surface has a rate of change of height with distance in the other orthogonal direction of greater than 5 percent.

In embodiments of the above described methods the two paths cross at heights differing by more than 50, 100, 150 or 200 metres. The paths may be linear, piecewise linear or curved and, in general, adjacent flight paths may be non-parallel by more than two degrees, three degrees, five degrees, ten degrees or more. As previously mentioned the paths, taken together, need not substantially lie in a surface, in particular because there is no "substantially constant height" constraint. However one convenient method of constructing a set of paths for such an airborne survey is to build a conventional survey and then modify this so as to more closely approach the ground in particular regions, such as valleys where better coverage is desirable.

A typical method of constructing a conventional survey is to perform a two-dimensional drape of a surface on which the survey is to be performed over the underlying terrain.

At high points of the underlying terrain the surface height is typically determined by a minimum permitted flying height for the aircraft (for safety), and then the surface height is reduced subject to aircraftdetermined constraints on rate of climb/descent, typically of order two percent to three percent. Because the paths are constrained by a rectangular grid these constraints apply in two orthogonal directions. By contrast with embodiments of the above described method such a constraint need only be applied in one direction (one dimension rather than two). In effect a one-dimensional drape may be employed to determine the flight paths. Thus in embodiments of the method we describe the set of paths, taken together, defines a surface in which lies a said path defining one or two orthogonal directions, in the other direction the surface being permitted to have a rate of change of height with distance of greater than a permitted climb/descent rate of the aircraft, for example greater than three percent, five percent, ten percent, 20 percent or more.

The invention also provides a data carrier carrying aircraft navigational data for such a set of airborne potential field flight survey paths.

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 andlor data may be distributed between a plurality of coupled components in communication with one another.

Further aspects of the invention are as follows: A method for conducting potential field surveys from the air whereby the aircraft is equipped with a range of geophysical measurement equipment including one or more potential field measurement instruments, for example vector gravinleter, gravity gradiometer, magnetometer, magnetic gradiorneter or other and in which the aircraft flies an irregular set of disconnected, non level, non straight flight lines in a manner which covers the survey area, preferably with the lowest possible surface to plane separation for each line flown and with approximately uniform coverage over the survey area (such lines being consistent with safe operation of the aircraft). Preferably in this method the earth potential field quantity of interest is determined using algorithms which work with the measured potential field measurement instrument data, collected along the survey lines as it is without the need to either level the data, to bring data to a common horizontal plane or to grid the data, such data including the measurement of the potential field quantity itself and the (instantaneous) position and optionally attitude

of the potential field measuring instrument(s).

A method of conducting potential fields surveys from the air as described above but with a flight pattern altered slightly with the goal of achieving an acceptable number of "crossovers" per unit area, here a "crossover" being a point where lines flown at different headings of the aircraft nominally cross at the same point in space.

A method for conducting potential field surveys from the air as described above where some lines from a set of essentially parallel but not necessarily straight lines, are deliberately made to meander so as to force a larger number of crossovers with other lines in that set which are flown parallel.

A method for conducting potential field surveys from the air as described above where all lines are flown with the goal of forcing a maximum (practicable) number of crossovers over the total survey area (given a constraint, say, on total path length/flying time). None of the survey lines will necessarily then be parallel to any other survey line, none of the lines will necessarily be straight and in general many of the lines will not traverse the whole survey from one side to the other.

A method for conducting potential field surveys from the air as described above but in which the aircraft flies a more conventional survey pattern, such as 2-D drape or constant height survey pattern or some other pattern which involves lines and tie lines which are flown with the intention that they intersect at every crossover. (For example the navigational data may define an exact intersection although in practice the pilot will generally only achieve this within a tolerance such as I Om or 20m).

A method for conducting potential field surveys from the air as described above in which the data is "levelled" before processing. Levelling here is a generic term which covers techniques, which include one or more of the following: noise reduction, removal of low frequency drift, matching low frequency content of neighbouring lines, referencing data to a fixed height plane and the like. The data may also be gridded before processing.

A method for conducting potential field surveys from the air as described above where the first stage of the analysis comprises the removal or adjustment of data which cannot realistically be caused by either geological structures under investigation or the geornorphology (the terrain) of the survey area. The effect of the terrain may also be used to correct the data at any stage of the processing before deduction of the gravity potential of the earth.

A method for conducting potential field surveys from the air as described above where an accurate DEM (digital elevation model) is produced using a combination of LIDAR (Laser radar) and an IMU (Inertial Measurement Unit) in conjunction with DGPS (Differential Global Positioning System) to correct the LIDAR data for the plane motion. The DEM and DGPS data may also be used to correct the measured potential field data for the terrain. Likewise aircraft acceleration, attitude, angular rate and angular acceleration data may also be used to correct the output data of the potential field instrumentation. Any onboard or remote sensor can be used to provide the position and motion information for the aircraft and/or the potential field instrumentation.

Preferably therefore the plane is fitted with any of a range of additional standard airborne geophysical survey instrumentation such as instrumentation for: GPS, DGPS, altimeter, altitude measurement, pressure measurement, hyperspectral scanner, an electromagnetic measurement (EM), a Time Domain Electromagnetic system (TDEM), a vector magnetometer, accelerometer, gravimeter, and other devices including other

potential field measurement devices.

A method of conducting potential field surveys from the air as described above in which the outputs from instrumentation on the survey plane are corrected using instrumentation in a fixed or movable base station, for example according to best practice at the time. Such equipment may include OPS and magnetic instrumentation and high quality land gravimeters.

A method for conducting potential field surveys from the air as described above in which data collected according to any of the above methods is combined with any ground based or satellite based survey data to help improve the analysis, such data including terrain, spectral, magnetic or other data.

These and other aspects of the invention will now be further described, by way of

example only.

Background

We first describe some background, helpful for understanding the invention.

Potential field data includes, but is not limited to, gravirneter data, gravity gradiometer data, vector magnetometer data and true magnetic gradiorneter data. Such data is characterised mathematically by a series of relationships which govern how the quantities vary as a function of space and how different types of measurement are related. Elements and representations of a potential field may be derived from a scalar quantity.

For gravity, the relevant potential is the gravity scalar potential, (r), defined as cI)(r)= flJ 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 force, which is how the gravitational field is experienced, is the spatial derivative of the scalar potential. Gravity is a vector in that it has directionality as is well known - gravity acts downwards. It is represented by three components with respect to any chosen Cartesian coordinate system as: ( \ (&1)(r) DcJ)(r) a(r) ax ay az Each of these three components varies in each of the three directions and the nine quantities so generated form the Gravity gradient tensor: ô OcD(r) 8 ô(r) 0 O'1(r) G G G OxOx aay OXOZ G - G G G 8 5((r) 8 acP(r) a 0(r) - Oy Ox Oy Oy Oy Ozx G G, G 8 OcI(r) 8 OcIr) a aci(r) Oz ax 5z Oy Oz Oz The mathematical theory of potential fields is well established - the fundamental equations and relationships follow from analysis of the properties of the scalar potential function, its derivatives, its Fourier transforms and other mathematical quantities. From one of Green's theorems, in the event that any of the spatial

derivatives of the scalar potential (including the scalar potential itself) are known over a closed surface, then the value of that spatial derivative is known at all points within the volume enclosed by that surface. A corollary to this is that once this quantity is known at all points, by differentiation and integration, all other derivatives of the scalar potential can be derived including the scalar potential itself. Thus the scalar potential and all of its derivatives are effectively known at all points within a volume when just one of its derivatives is known over a surface closing that volume. This indicates that a full measurement of any component of any of the derivatives of the scalar potential allows the calculation of any other component of any derivative of the scalar potential.

Following from it also does not matter, at least in theory, which quantity is measured, the choice of instrumentation comes down simply to which instrument measures the desired quantity with the largest signal to noise.

Differentiation of the gravity scalar potential above (after some issues associated with where r -+ 0) ultimately yields: V (r) 2 a2 a2 = -4,rGp(r) ax ay- which, in regions where there is no matter reduces to Laplace's equation, an important fundamental relationship in gravity: ô2 (r) ô2(r) 52r)0 ox2 Harmonic functions satisfy Laplace's equation and they have many properties which may be utilised in the analysis of data collected from potential field surveys.

Conventionally, data is analysed and processed using a range of techniques all of which work with the data collected from the survey as a starting point but which thereafter alter both the data and/or its format so the values associated with the measured quantities all appear on a regular 2-D grid which is on a horizontal, fixed altitude analysis plane, the processes used being known as "levelling" and "gridding".

Gridding is a standard geophysics numerical processing technique but one which is handled slightly differently by everyone who uses it. The basis idea is * to break the survey area up into rectangular cells whose sides are usually aligned to the principle directions flown for the survey, * to replace the actual measurement data with data which is entirely "equivalent" to the measured data but which is now assigned values at points in the middle of each cell.

Such data is called "gridded" data. There is much debate and there are many ways to ensure that the data "invented" as being representative of the measured data is indeed "equivalent" to the measured data. There are also many ways to choose the dimension of each cell but they all relate to the average separation of lines flown in the two orthogonal directions. Where survey lines are not flown in a regular manner, as proposed in this patent, the concept of gridding, although not necessary, could still be used but there is much more scope for choice on the cell size and orientation.

Once the data is in this format, it is much more tractable mathematically and this is what is used for all cuiTent potential field surveys. There is no requirement on any of the data involved in these processes to satisfy Laplace's equation. The data is treated as set of numbers and statistical and other methods are applied to the numbers to give the best estimate of the potential field on the horizontal analysis plane.

Data is generally reduced to be a 2-D Fourier series so, in addition to the requirement that the data is at fixed altitude, there is the requirement that each line of data has 2" data points in order to make the Fast Fourier transform method work, and that the data is collected in orthogonal directions.

The general form of the 2-D Fourier series representation of any quantity is well known and for gravity gradient can be represented as a sum of 2-D spatial sine waves in the form g (x, y) = g_, (k,,,, k,, )euIdhhi(ent)? in j where the wavenumbers k,,,, k,, are related to the size of the survey, L,, L, in the x, y directions respectively by 2r, 2ir r r This expansion is only valid at constant height and we show explicitly below the wavenumber and height dependence of the coefficients g(k,,,,k,,), which collectively form the 2-D Fourier transform of the gravity gradient.

Detailed description of preferred embodiments

Data does not have to be handled as above if one uses known properties of the gravity or magnetic potential fields. The discussion below will be based upon gravity but can

easily be extended to magnetic fields.

First we will state some observations relating to the gravity field: 1. The gravitational field outside a body can be modelled in all respects as if it came from matter situated entirely in a vanishingly thin layer at the surface of the body and which accurately follows the surface of the body. Such a layer defines an equivalent source - i.e. a source of gravity which produces substantially (theoretically exactly) the same gravity signature as does the body itself. There are many ways in which equivalent sources can be defined, they can be as above or they can be strictly horizontal, they can be all or partly above or below the surface of the earth or other but they all have the same property; namely, they generate the same gravitational field as is generated by the earth.

For further information reference may be made to R.J. Blakely, "Potential Theory in Gravity and Magnetic Applications", Cambridge University Press, 1995.

2. The gravity field outside a body can also be uniquely written as a modified 2-D Fourier basis function series which is valid at all positions outside the body.

This function is explicitly harmonic and is reproduced below. Its form is completely general and can represent any possible mass distribution within the body.

In preferred embodiments we measure G as a function of position, r,,,eas,, rc, using a gravity gradiometer and work with this without needing to generate the other elements of the gravity gradient tensor. This can be used to generate a representation of the underlying mass distribution. Here we describe two techniques which may be employed, an equivalent source technique and a modified Fourier basis functions technique. The above mentioned Green's theorem shows that the underlying mass distribution can be derived, in principle from G, although in the techniques we describe Green's theorem is not used explicitly.

Equivalent Source Method In this embodiment the surface of the survey area is broken up into small pieces, typically of order 50m on a side - these are called platelets (later, mass elements). It is easy to forward calculate the gravity from each platelet, Blakely (ibid), the mass of which is adjusted until the best overall fit to the measured data is obtained. This mass detemiination is done via a standard least squares fitting procedure. The fit is obtained by matching the data at the true measurement position to the gravity field generated by the proposed equivalent source at the identical true measurement positions. This process is mathematically rigorous and does not introduce any artificial adjustments to the data in order that it conforms to a horizontal rectangular survey.

Once this fit is obtained, the fit is deemed to be the primary data set. All subsequent analysis to determine geological structure compares and minimises differences between the gravitational field that any given geological structure would generate with that generated by the equivalent source. One real advantage of the technique is that the best fit conies from a mass distribution, albeit a synthetic one, and therefore the best fit solution will automatically satisfy Laplace's equation. This is unlike the conventional method as described above, which produces a numerical best fit but which does not impose the added restriction that the data has to satisfy Laplace's equation, i.e. that it could come from a real mass distribution.

The equivalent source method does not have to use a surface conformal to the topography, it can use sources which cover any surface which can be at constant altitude, above or below the earth's true surface, can cut through the earth's true surface and so on. The choice of a surface following the topography is likely to produce less variation in mass of the individual platelets but the overall result is not, in principle, dramatically affected by any reasonable choice of surface.

Likewise, the platelets can be of any size or shape and they do not even need to be equal in size - indeed, it helps the efficiency of the analysis if their size and geometry is allowed to vary depending upon how fast the topography is varying in each region.

The mathematical complexity of the process is determined purely by the number of platelets used and how many of these platelets are used in the analysis at each position on the survey. One advantage of this tecimique is that for some components of gravity or gravity gradient, it is possible to use just platelets in the region of the data point and this dramatically reduces the complexity of the analysis. It is clear that this analysis is one in the space domain.

Once the equivalent source is generated, it is then possible to predict any derivative of the gravity scalar potential on any surface by direct forward calculation. This process is helpful both from an analysis and from a visualisation perspective.

In more detail, given the masses of each surface element, a straightforward calculation is used to predict what value would be obtained for the measured quantity at every measurement point, be that quantity a component of the gravity vector or of the gravity gradient tensor. In general, this will be a summation of the form shown below, Here we use gg as notation for the measured quantity which, as noted above, is G in some preferred embodiments.

ggcakulate, (mecsure) = ass-elemelu F(r,,,eas,gre - ass_eje,,() Equ (I) all-masses In the above equation F is called a Greens function (Blakely, ibid, at page 185, incorporated by reference) and r,,,aeIe,,,e,,t defines the location of the mass element (for example the centre of gravity or some other defined point). The function F is standard and known to those skilled in the art (and niay be deduced for may different mass element geometries). For example the Green's function F for a rectangular prism is defined in Blakely, ibid, at page 187, and is incorporated by reference; this has 8 terms each of which corresponds to a vertex of the prism.

The important thing about this relationship is that function "F" is known and it is the masses of each element which are not. The calculated value, in general, will not be the same as the measured value and one constructs the following sum S = [(caicuiaiei (r,,,eas,irc) - gg,,:asnrd (mcosure))2] Equ (2) all -nieasure,neisis The quantity S in effect defines the total square error S (noise and modelling errors), and by minimising this a best fit ggca1cIlatc( can be found. This can be done by standard procedures, for example using Matlab (RTM), By extension to the above argument, for the quantity S the only unknowns are the actual masses of each mass element, and of course the measurement noise associated with the measurement at each point. The best fit, that is the best estimate of the values of these mass elements, is obtained by requiring S to be a minimum with respect to variations in the masses of each element. When the correct coefficients are obtained, the quantity S is merely the sum of the squared noise terms associated with each measurement and this is the minimum value which can be obtained by S. By differentiating S with respect to the mass of each mass element, that is, the first derivative of S with respect to,,,ase/epi,c,,1 must be zero for every mass element. This ultimately reduces to the generation of a set of simultaneous equations in which the only variables are the masses of each surface mass element and where the number of equations equals the number of mass elements involved in the fitting procedure.

Conveniently this set of simultaneous equations can be represented by a matrix, and solved for values of the surface mass elements. Since the number of measurements is generally much greater than the number of simultaneous equations (i.e., the number of mass elements), for example at least five or ten times greater, the set of equations is overconstrained. This is useful because if this were not the case the values of the surface mass elements would fit to the noise; instead an optimisation procedure, for example a least squares fit, is used to determine the values of a surface mass elements, thus reducing the noise. Moreover because a physical model gravitational field (real mass elements) is used the equations obey Laplace's Equation so that it is the gravitational field rather than the noise which tends to be fitted..

Once the mass elements have been found a forward calculation (ie summing the effects of the mass elements) allows other scalar potential components to be derived and, by differentiation other components of G. The deduced values of G can be compared with a geological model (referred to as an "interpretation") to determine the underlying geological structure.

Modified Fourier Basis Functions Method In a second embodiment of the method a modified Fourier basis functions technique is employed. In this method a two-dimensional Fourier expansion is used to represent the measured data but each two-dimensional coefficient A,,,,, of each component (A,,,,,e1k1eikh1)) of the 2-D Fourier transform of the gravity field is rewritten as a product of a position independent coefficient, (C,,,,,) multiplying a height-dependent coefficient (e'm), That is, A,,,,, = Explicitly, one can write the gravity

field of any mass distribution in the form

gg(x,y,z) = E C,,,,,e' e'e" Ik11I,J,,2 + k) in ii In the language of Fourier transforms, A,,,,: is the 2-D Fourier transform of the gravity field. This representation is completely general and each component of the expansion representing the gravity field explicitly obeys Laplace's Equation (when differentiated the multipliers of the height-dependent coefficient and of the 2D components cancel).

What is now made explicit is the height dependence of each Fourier coefficient and it is also clear that the representation of the gravity field afforded by knowledge of the C,,,,, coefficients allows the complete determination of the gravity field at any point in space where there is no mass. The expansion implied by this modified Fourier series can be truncated at a point dependent upon the desired accuracy of the representation.

At position r the measured potential field (gravity or gravity gradient) is known and, for a particular pair of values k,,, and k,,, known as wavenumbers, the value of each of the modified Fourier series basis fijnctions, e1Te1e, is also known. Each wavenumber describes the harmonic spatial variation of the Fourier basis function, for example, every coefficient with a given wavenumber k,11 = 2r/2,, will be periodic in space with a period 2,,,. A further property of the Fourier method is that, for each of the independent directions, every wavenumber will be a integer multiple of a minimum wavenumber which corresponds to a wavelength broadly equivalent to the size of the survey Other expansions can be employed, for example spherical harmonics, which are useful for larger areas, in particular where the curvature of the earth becomes significant.

Broadly speaking, coefficients for the expansion are determined in a manner virtually identical to that described above for the equivalent source method. Once one has chosen the maximum wavelength of interest for each direction, and of the orientation and origin of the coordinate system for measurement, the value of the basis fUnctions can be determined at every point in space. Once this is done, the same difference squared function, S, is calculated for the measured quantity and the best values for the coefficients C,,,,, are determined by nhinimising such function with respect to variations in C,,,,,. Thus, as with the first embodiment described above, a set of simultaneous equations is generated which can be solved for the Fourier expansion coefficients, These coefficients can then be employed to determine (by forward calculation) the gravity or gravity gradient on any two- dimensional surface, that is to provide a map. it will be recognised that with this (and the preceding) embodiment the data can be randomly sampled since the data points are effectively disconnected - all that is needed is a set of measurements and associated cx, y, z,) positions in order to be able to determine the gravity or gravity gradient field on any other set of points or indeed on a plane, that is a map Furthermore it can be seen from this that there is no need for padding of measured data to build up a rectangle. Conceptually the data points can be treated as an arbitrarily shaped portion of a two-dimensional region of maximum dimensions determined by maximum wavelengths for the Fourier expansion in the x and y directions (although, as can be seen, there is no need for the measured data to lie on a twodimensional surface). It can be seen, therefore, that values for k,, and k can be chosen independently of the survey dimensions (here k being wavenumber, that is 2it/? whereas with a conventional survey the maximum wavelength must correspond to either the length or breadth of the rectangular surveyed area).

In more detail, the concept behind the Fourier series route is to solve the problem of Equation (1) in the Fourier domain rather than in the space domain In what follows, we will presume that g is the measured quantity (also referred to above as Ga).

The gravity field from any body or mass distribution can be written in the completely general form: g (x,y, )= g'(/c, , k,, )mxe1uYe_Ikmnk Equ (3) Hi H which is always valid and where = .Ik,,,2 +k112 In Equation (3) the coefficients g, (k, ,k,1) are formally independent of height and wavenumber, they depend only upon the geology under investigation. Note the formal similarity of Equation (3) to the general form of the 2-D Fourier series representation given above. The expressions are identical, which they have to be, only if k k i- k k -jk,,Iz 5:z'. in' nJ6zz ". us' ;jP What we see is the strong dependence upon height and of wavenuniber of the standard Fourier transforms, that is, of the coefficients g,, k,,, , k,,), This is what makes it impossible to do a standard Fourier transform of data collected on a variable height survey. It also illustrates why it is impossible to "upward" continue data on a point by point basis; upward continuation only works if the precise frequency content of the quantity is known. It can also be seen that each component of the expansion of Equation (3) has a different height dependence.

However, the modified Fourier basis function method we describe does not attempt to perform a Fourier transformation per se - it is, as described above, not a mathematically sound approach.

The process here is that every data point is treated on a similar basis and that the quantity e1kmeikt 1Hkn,uk can be assigned a value which is accurate once the measurement position is known. This is because we chose the k-values as described later, and r(x,)',z,.) is known. Essentially maximum and minirntim values for k and k are selected, for example by notionafly putting a rectangle around (or slightly inside) the survey area - that is which approximately represents the survey area - and the dimensions of the rectangle define the maximum wavelength (k-values). The minimum values may be chosen, for example, dependent upon the approximate average height of the survey, for example a multiple of the average survey height between 0.5 and 3. This is because the length scale of the variations depends on the survey height and, broadly speaking a minimum value is selected to represent variations which are expected or seen at the average survey height at greater than a threshold significance, taking into account computing power. To illustrate the concept, the exp(- kz) term above at z = survey height has attenuated by approximately e6 for wavelengths equal to the survey height.

We then have a set of equations based on Equation (3) equal to the number of data points used whereby the measured values of the gravity gradient, including noise, and the quantity e''e'e" are known at every data measurement point. From that point, the only unknowns are the coefficients g' (k, , k,,) and the fitting process clearly then becomes only one of obtaining the best estimates of the coefficients.

Referring back to Equation (3), this gives the calculated gzz, and by analogy with Equation (1), a version of Equation (2) can be derived. This provides a set of equations equal to the number of unknowns, which in this case are values of the coefficients g'.

This is done in the same way as for the equivalent source method, i.e. the coefficients are determined by a least squares fitting procedure. The important issue is that, in this form, the coefficients g,, (k1 k) explicitly have no height or wavenumber dependence.

This technique has a similar advantage to the equivalent source distribution in that the solution is harmonic because each term in the expansion of Equation (3) is explicitly harmonic, Instrument noise will not be harmonic so, this technique, even in isolation, has potential to remove much of the noise in the process.

This technique, like the equivalent source method, does not need regular gridded data, One additional surprising advantage of the technique is that it is not, like the Fourier series method itself, tied in the wavelengths which can be used for the analysis. This is a subtle result depending upon the validity of the modified Fourier Series method in completely general terms over whatever area might be of interest. The process involves choosing the maximum wavelengths of interest in the two orthogonal directions of the survey plane, which leads to a maximum wavelength of interest in each direction, 2xniax' j'max. This wavelength manifests itself as a minimum wavenumber in each direction given by - 2n' 2,r "xmin - kyin = 2 xmax ymQX All wavenumbers in the expansion of Equation (3) are integer multiples of these. That is: k,,, = in Ic mm, k. ky mm It is further of advantage that this method does not have to have a power of 2 as the number of data points to be used in the analysis of each line, as is needed for the fast Fourier transform method. This avoids the need for data padding to achieve such a situation.

Unlike the equivalent source method, where the analysis is done in the space domain, Fourier series analysis is done in the frequency domain. By frequency is meant spatial frequency which is a measure of how many periods a given component has per unit measure of distance and this is represented by the wavenumbers introduced above - the larger the wavenumber, the larger the number of periods per unit of distance and the shorter the repeat distance. The proposed technique does not explicitly work in either domain, it merely generates a set of equations from which the best estimates of the unknowns can be calculated.

Fourier methods, in practical tenns, usually require measurement over a rectangular area and the minimum wavenumbers are the longest wavelengths which can be used to represent the data - the longest wavelengths are explicitly the size of the survey in each direction. Where a survey does not conform to this criterion, data is usually created so that the revised data set does meet it in a process called padding.

The modified Fourier basis function method has many advantages compared to conventional Fourier analysis.

1. It is not necessary to perform any Fourier transformations at all - we know what the standard format is going to be and the coefficients in the expansion of Equation (3) are to be determined via a best fit process rather than calculating the Fourier transforms explicitly.

2. The wavenumbers to be used do not have to be integer related to the size or shape of the survey. This is because we know the form of the Fourier transformation and we are determining the coefficients via a best fit process as just described. 3.

Perhaps the most striking result is that the survey data can be collected at a non constant height. 2-D Fourier transforms as used conventionally, must be performed at a constant height and the Fourier transforms introduced in Equation (3) are strong functions of height. This functionality is explicitly expressed in Equation (3)and the coefficients g.., (k,,? ,k,1) are no longer height dependent.

Summary

In all, embodiments of the methods we describe allow the full use of the actual measured data and both are mathematically rigorous; they do not corrupt the data by making it fit to a regular constant height grid, the solutions generated are harmonic and the data does not have to have a power of two as the number of data points in two orthogonal directions.

The process is to use all data from every data point and try to generate an equivalent surface mass distribution or a set of modified Fourier series coefficients which best fit the data over the area of the survey.

This is fundamentally different to the normal processing method which assumes the data is from a fixed height, even if it is not, and generates new data from the measured data for the sole purpose of generating data on a regular 2-D grid. This "gridded" data set assigns values to the measured data which are not measured and are assigned to positions where the aeroplane did not collect data. Although this is clearly incorrect, it is what is done because of the advantages afforded by the various processing algorithms used to analyse the data.

The errors generated by such a process are likely to be more serious if the fixed height approximation is poor in a given situation. They are also serious if one is trying to fly as close to the ground as possible to maximise the signal from small features relatively close to the surface of the earth - these features giving rise to short wavelength components in the modified Fourier series expansion. Finally, the errors are less serious when the instrument noise is large as the errors involved in the extrapolations are then smaller than those associated with the instrument.

For conventional airborne gravity it was not necessary to fly very close to the earth because the instruments were relatively insensitive and additionally were set up to be sensitive only to long wavelengths. Surveys could be then be flown at more or less constant altitude or at a height where the data could be reasonably well assumed to be corning from a constant altitude.

However, for any survey using the new high sensitivity gravity gradiometer systems which are also capable of resolving much shorter wavelengths, it is desirable to fly as close to the ground as possible in order to maximise the signal, particularly that contained in the short wavelength components. Such surveys are not likely to adequately satisfy the constant altitude criterion but all current analysis methods proceed as if they do.

It is, however, possible to utilise the equivalent source or modified Fourier transformation method to process the data and combine this with a flying method which

is completely novel for potential field surveys.

In its most general form, the method is to fly a completely disconnected set of survey lines, collect the data from the lines and, using both the data and its absolute position of acquisition, obtain best fit parameters to either the equivalent source platelet masses or the components of the modified Fourier transform method.

The survey lines do not need to be straight, they do not need to be horizontal, they do not need to be flown in orthogonal directions, lines in approximately orthogonal directions do not need to intersect. Ideally the flight lines will follow the terrain as closely as possible without any other constraints. Of course, any other flying methodology which sitsbetween this survey flying method and those involving constant altitude orthogonal intersecting flight lines is possible and any survey could involve elements of either flying method. It is the power and mathematical correctness of the analytical techniques which make these flying procedures possible.

Data collected using this flying method can clearly be additionally processed before or after the use of either the above processes. Any of the existing techniques for noise reduction or data combination can be utilised. Topographical effects can be largely eliminated by use of high quality DEMs, indeed any of the processes associated with more conventional data processing can be used.

The aeroplane can be any geophysical survey craft which is suited to low flying. In practice, the plane will be fitted with as much equipment as it can carry commensurate with safety and reasonable flying time. The equipment chosen will measure signals associated with the geology and signals associated with the plane motion. The latter signals are used, particularly in gravity related surveys, to correct the measurements taken by the instrumentation chosen to measure geological properties.

In flat areas, all flying strategies converge but it is in hilly or genuinely mountainous areas that the proposed method generates significant advantages. In such areas, each part of the survey area is analysed relatively independently and a strategy is developed as to how to fly the area.

In such situations, the area for survey is analysed for its topography and a flight survey strategy is developed based upon the geological targets of interest and how it is possible to fly the area to best identify those geological targets. If, as is usually the case, it is necessary to get very close to the ground at all points of the survey, then the survey lines will become disconnected with relatively little connectivity. The plane can only climb or descend at a given rate and adequate clearance is needed to get safely over ridges. Lines flown going "down" or "up" valleys will not intersect lines "crossing" the valleys or will intersect at many fewer crossovers. Crossovers are usually considered important to get rid of low frequency noise (often called drift) and this is another consideration on designing the survey.

Equipment desirable for a standard survey include * Gravity gradiometer * Gravimeter * Inertial measurement unit (ll4U) * GPS * LIDAR * Magnetic gradiorneter * Magnetometer * Altimeter * Radar height meter * Barometer * Spectrometer * 1-lyperspectral imager * EM * TDEM * Navigation S/W * Data collection Once the flight strategy is developed, the actual survey is flown and the data is collected. The analysis uses the real data from the potential field measurement system, this data is corrected for any known artefacts associated with the measurement (e.g. tilt of plane, accelerations of the plane, pressure affects and so on), and then the data is analysed as indicated above. The reason for flying low is because of the rapid height dependence of the measured signal, particularly those having short wavelengths, i.e. large wavenumber.

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.

Claims (24)

1. CLAIMS: 1. A method of processing measured potential field data from a

   potential field survey of the earth to determine map data for mapping said field, the method comprising: inputting said measured potential field data, said measured potential field data comprising data defining a plurality of potential field measurements and associated positions, each said position defining a position of a said potential field measurement in three dimensions; determining a plurality relationships between said potential field measurements and said positions, each said relationship relating a said potential field measurement to a function of a said associated position in three dimensions multiplied by a field mapping parameter; and determining a substantially self-consistent set of said field mapping parameters for said plurality of relationships to thereby determine said map data.
2. 2. A method as claimed in claim I wherein said map data comprises said set of

   field mapping parameters.
3. 3. A method as claimed in claim 1 wherein said map data comprises data mapping said field on a surface, the method further comprising determining said field on said

   surface from said set of field mapping parameters.
4. 4. A method as claimed in claim 1, 2 or 3 wherein said relationships comprise a set of simultaneous equations and wherein said field mapping parameter determining comprises solving said simultaneous equations.
5. 5. A method as claimed in claim 4 wherein said simultaneous equations are over- constrained, and wherein said solving comprises solving by a noisereduction procedure.
6. 6. A method as claimed in any preceding claim wherein said position function substantially satisfies Laplace's Equation.
7. 7. A method as claimed in any one of claims I to 6 wherein said field is represented by a plurality of surface elements each having an associated source strength, and wherein each said field mapping parameter comprises a said source strength.
8. 8. A method as claimed in any one of claims 1 to 6 wherein said field is represented by a harmonic expansion having a plurality of expansion coefficients, and wherein each said field mapping parameter comprises a said expansion coefficient.
9. 9. A method as claimed in claim 8 wherein in said representation said coefficients are substantially independent of height above the earth.
10. 10. A method as claimed in any preceding claim wherein said measured potential field data comprises one or more of gravimeter data and gravity gradiometer data, and

    wherein said field represents a gravity field.
11. 11. A carrier carrying processor control code to implement the method of any preceding claim.
12. 12. A data processing system for processing measured potential field data from a potential field survey of the earth to determine map data for mapping said field, the system comprising: data memory for storing said measured potential field data, said measured potential field data comprising data defining a plurality of potential field measurements and associated positions, each said position defining a position of a said potential field measurement in three dimensions; program memory storing processor control code; and a processor coupled to said data memory and to said program memory, to load and implement said control code, said code comprising code for controlling the processor to:

    input said measured potential field data;

    determine a plurality relationships between said potential field measurements and said positions, each said relationship relating a said potential field measurement to a function of a said associated position in three dimensions multiplied by a field mapping parameter; and determine a substantially self-consistent set of said field mapping parameters for said plurality of relationships to thereby determine said map data.
13. 13. A method of processing measurement data from an airborne gravity survey to provide data for a gravity map, the measurement data comprising a plurality of gravitational potential field measurements each with an associated measurement position in three dimensions, the method comprising using said measurement data to estimate coefficients multiplying a function of said three-dimensional position in a harmonic expansion to thereby determine a representation of a gravity field in an area of said survey for said gravity map.
14. 14, A method as claimed in claim 13 wherein said harmonic expansion comprises a Fourier expansion and wherein said function of threedimensional position comprises a function having a form equivalent to where x, y and z define said position in said three-dimensions, km and k are positive, and Ik11, = +
15. 15. A method as claimed in claim 13 or 14 wherein said coefficients are formally independent of height.
16. 16. A method of processing measurement data from an airborne gravity survey to provide data for a gravity map, the measurement data comprising a plurality of gravitational potential field measurements each with an associated measurement position in three dimensions, the method comprising using said measurement data to estimate a plurality of mass elements each multiplying a function of said three- dimensional position in an equivalent source representation of said gravity field to thereby provide said data for said gravity map.
17. 17. A method of mapping a gravity field, the method comprising: inputting measurement data from an airborne gravity survey, the measurement data comprising a plurality of gravitational potential field measurements each with an associated measurement position in three dimensions; determining a plurality relationships between said potential field measurements and said positions, each said relationship relating a said potential field measurement to a function of a said associated position in three dimensions multiplied by a field mapping parameter; and mapping said gravity field by determining a substantially self-consistent set of said field mapping parameters for said plurality of relationships.
18. 18. A method as claimed in any one of claims Ito 10 and 13 to 17, further comprising surveying said field using an aircraft to capture said measured data for said determining of said plurality of relationships, and wherein said surveying comprises flying a set of paths which have one or more of the following features: two paths cross at heights differing by at least 50 metres; in a region of said survey paths in the same general direction are non-parallel by greater than 5 degrees; said paths include curved paths; said paths of said set of paths, taken together, do not substantially lie in a surface; said paths of said set of paths, taken together, define a surface wherein at least one of said paths defines one of two orthogonal directions in said surface such that said surface has a rate of change of height with distance in the other orthogonal direction of greater than 5 percent.
19. I 9. A carrier carrying processor control code to implement the method of any one of claims 13 to 18.
20. 20. A method of flying an airborne potential field survey, the method comprising flying an aircraft along a set of paths and measuring potential field data at points on said paths, wherein said set of paths has one or more of the following features: two paths cross at heights differing by at least 50 metres; in a region of said survey paths in the same general direction are non-parallel by greater than 5 degrees; said paths include curved paths; said paths of said set of paths, taken together, do not substantially lie in a surface; said paths of said set of paths, taken together, define a surface wherein at least one of said paths defines one of two orthogonal directions in said surface such that said surface has a rate of change of height with distance in the other orthogonal direction of greater than 5 percent.
21. 21. A carrier canying navigational data defining a set of airborne potential field flight survey paths as claimed in claim 20,
22. 22. A method of generating airborne flight survey data, in particular for the method of claim 20, the method comprising: determining a flight survey pattern by draping a surface over the terrain underlying the survey in accordance with an aircraft maximum climb/descent rate condition; and generating flight survey data for a plurality of paths substantially on said surface; and wherein the method further comprises adjusting one or both of said surface and one or more of said paths to reduce the height of one or more of said paths such that in at least one direction in said surface orthogonal to a said path said maximum climb/descent rate condition is exceeded.
23. 23. A carrier carrying processor control code to implement the method of claim 22.
24. 24. A data carrier carrying the flight survey data of claim 22.