RU2486549C2 - Processing of gravity survey data - Google Patents

Abstract

FIELD: measurement equipment.

SUBSTANCE: method of processing of measured data of a potential field in aviation or marine surveys of the potential field includes identification of an array of field mapping parameters for field mapping and introduction of potential field data into the processor. Measured data of the potential field contains data that identify multiple measurements of the potential field, each with the appropriate position of measurement and time of measurement. The array of field mapping parameters is identified using a model containing a combination of a spatial part displaying spatial variation of the specified potential field, and a time area displaying time noise in the specified measured data of the potential field. Besides, the specified identification includes adaptation of measured data of the potential field both to spatial and time parts of the specified model.

EFFECT: problem of low-frequency noise in process of board survey of a potential field is solved.

27 cl, 5 dwg

Description

Technical field

The present invention relates to improved methods for processing potential field measurement data during aerial surveys, for example during gravimetric surveys, and to methods, apparatus and computer program code for such techniques.

State of the art

The present description relates to aerial surveys and, more specifically, to gravimetric surveys. However, the described techniques are not limited to these shooting examples and can be applied to other surveys of a potential field, including but not limited to magnetic field surveys, for example magnetotelluric surveys, electromagnetic surveys, etc.

A potential field is taken by measuring potential field data, which, for gravimetric surveys, may contain one or more gravimetric data (gravitational field measurement) or gravitational gradiometer data (gravitational field gradient measurement), vector magnetometer data, reliable magnetic gradiometer data, and other types data well known to those skilled in the art. The overall goal of a geophysical survey of a potential field is to find characteristic features that potentially indicate valuable mineral deposits.

Traditionally, aerial surveys of potential fields, such as gravimetric surveys, are carried out along a grid. The grid is defined by orthogonal rows of parallel lines (flight paths) on a two-dimensional surface that envelops the earth's surface. However, the enveloped surface is limited to the approximation to the Earth permitted for the aircraft and the maximum speed of the aircraft ascending / lowering. Some improved techniques for aerial surveys of potential fields that facilitate the collection of data at a short distance from the Earth are described in the PCT pending patent application "Processing of gravimetric survey data", PCT / GB2006 / 050211, thereby fully incorporated into this description by reference.

After the data of the potential field is collected, but before the interpretation of the data, the correction of the terrain is usually applied, compensating for surface variations in height. The surface profile can be obtained in the form of digital elevation data or determined from (D) GPS ((Differential) Global Positioning System) and / or from on-board techniques such as LIDAR (Laser Detection, Imaging and Distance Sensing System) and SAR (radar with synthesized aperture). Aircraft acceleration, orientation, angular velocity, and angular acceleration data can also be used to correct the output of potential field measurement equipment. The authors describe some improved techniques for correcting terrain during geophysical surveys in their pending United Kingdom Patent Application "Terrain Correction Systems", No. 0601482.3, filed January 25, 2006 and also incorporated herein by reference in its entirety.

However, despite improved techniques, the authors previously noted that low-frequency noise, or drift, in on-board measurements during shooting remains a problem. (Although the term drift sometimes refers to random monotonous variation, in the present description it is used to denote any kind of correlated noise with a characteristic frequency much lower than the frequency of the main signals of interest.)

The term “leveling” is used in the art as a general term referring to conventional techniques for matching data. These techniques include removing the low-frequency drift, matching the low-frequency contribution of adjacent lines, and assigning data to a plane of fixed height. For example, the intersection points of a standard coordinate survey can be used for cross-leveling, when the data on the survey lines are adjusted to minimize differences at these points.

However, there is a need for improved data processing and, in particular, for improved noise processing.

SUMMARY OF THE INVENTION

Thus, in accordance with a first aspect of the invention, there is provided a method of processing measured potential field data from an aerial or marine survey of a potential field to determine an array of cartographic field parameters for mapping a field, comprising: introducing said measured potential field data, said measured potential field data comprising a plurality of measurements defining potential field data, each with a corresponding measurement position and measurement time; and determining said array of cartographic field parameters using a model comprising a combination of a spatial part representing a spatial variation of said potential field and a temporal part representing temporal noise in said measured potential field data, said definition comprising fitting said measured potential field data to both of said , spatial and temporal, parts of the mentioned model.

In some preferred embodiments, the spatial and temporal parts of the model comprise the spatial and temporal terms of the model equation, which are defined together to determine the cartographic parameters of the field. Preferably, the model is linear; this simplifies the inversion. In particular, the model may have the form

,

,

where ρ and λ are the vectors of the spatial and temporal model parameters, respectively, A and B are the matrices, and f contains the vector of direct calculation of measurements estimated from the model.

Preferably, the model has a sufficiently high temporal resolution for the intervals between measurements, less than the expected correlation time of the device or the measurement noise. Thus, it is preferable that the time between the nodal points is selected depending on the measured, or expected, characteristic or correlation time of the noise in said measured data; this can be obtained from the specifications of the manufacturer of the measuring device.

In some preferred embodiments, the spatial portion of the model satisfies the Laplace equation. Preferably, the spatial portion of the model comprises an equivalent source model.

Preferably, the temporal part of the model contains a line defined by discrete points in time (nodal points). In embodiments, the line comprises sets of piecewise linear sections connecting these nodal points; in other embodiments, higher order interpolation functions may be used. Preferably, the nodal points are in time intervals of at least 60, 120, 180, 240, 300, 360 or 420 seconds. Preferably, the noise model displays a measured drift (which should not be monotonous).

Preferably, the determination of the cartographic parameters of the field comprises a solution to a system of joint equations. Since the equations may not be correct, regularization is preferably used, but preferably the spatial and temporal parts of the model are regularized separately to facilitate the difference between drift and true spatial variations. Thus, in some preferred embodiments, gradient regularization is used for the spatial part of the model, and another type of regularization, for example Tikhonov regularization, is used for the time part of the model. This helps, for example, to distinguish between true spatial variation and instrumental drift, which may, for example, be in different spatial directions in adjacent survey lines carried out in opposite directions. It will be clear that, in embodiments of the techniques, cartographic field parameters are described that are determined by fitting to the model using data from potential field measurements from many different survey lines at the same time (in the sense that fitting the measured data to a model defines a global fit for multiple survey lines, in embodiments for all processed measurement data).

Thus, in embodiments, the fit comprises jointly tailored potential field measurement data from a plurality of different survey lines for (aerial) surveying a potential field extending along the lengths of the lines. These survey lines may be substantially parallel or substantially not parallel, for example substantially orthogonal. In embodiments, the survey lines contain more than two, and preferably substantially all of the flight paths when shooting, in other words, in preferred embodiments, substantially all of the data collected from the survey is provided for a model that models both spatial and temporal variation for full shooting. (In this latter case, the survey preferably comprises a plurality of flight lines positioned so as to encompass the surveyed area, which may be of any reasonable physical length.)

The measured potential field data preferably contains gravitational data, for example gravimetric data (gravitational field measurement), and / or gravitational gradiometric data (gravitational field gradient measurement), and may include multiple components. For example, the measured potential field data may contain a component G _{zz of the} gravitational gradient. Thus, in embodiments, the mapped field contains a gravitational field, although this can be displayed using gravitational gradient data. However, in other embodiments, the measured potential field data may comprise vector magnetometric data, true magnetic gradiometric data, or other potential field data (e.g., scalar gravimetric data).

The cartographic field parameters preferably contain parameters or coefficients that allow characterizing the potential field and, in particular, which can be used to create similar data of the potential field from direct calculation in the embodiments.

Thus, in embodiments, the parameter values determine a potential field, for example, a gravitational field or a gravitational gradient field, which displays the best approximation to the measured potential field data. This best approximation preferably minimizes the remainder, which depends on the difference between the direct calculation of the potential field from the model and the measured data of the potential field, for example, the mean square error, the value of the difference modulus, or some other remainder.

In embodiments, the cartographic field parameters are extracted from the spatial part of the model and may contain values for equivalent mass source elements or, in general, equivalent source elements with corresponding source strength when the cartographic field parameters can contain these source forces.

In embodiments, the method further comprises defining a map using an array of cartographic field parameters, and it can display either potential field data or data extracted from potential field data, for example, after further processing. Such a map can display either two-dimensional or three-dimensional data, and in embodiments, the map is constructed by performing direct calculation to determine a potential field, for example, on a surface, from an array of cartographic field parameters. The cartographic parameters of the field can be used in direct calculation, for example, to allow the determination and mapping of a potential field, for example, a gravitational field, or a gravitational gradient field, or the corresponding scalar potential.

One skilled in the art will appreciate that in the general case, when measuring potential field data, gravity and / or gravity gradient are actually measured (although other quantities derived from the spatial derivatives of the potential field can be measured additionally or alternatively).

In some preferred embodiments, a method that uses an equivalent source mapping, surface mass values (or density) for a plurality of surface mass elements, preferably on a flat surface (but not necessarily, since the surface can be chosen to approximate the surface of the earth in the shooting area), used in direct calculation to map one or more field components on a flat cartographic plane or some other convenient surface.

Another related aspect of the invention provides an equivalent source method for processing measured potential field data for determining field cartographic data in which said potential field measured data is modeled as a combination of a spatial signal and a temporal noise component, the method comprising determining the overall best approximation of said potential field measured data to said simulated combination.

The invention further provides a processor control code for implementing the above methods, in particular on a storage medium such as a disk, CD or DVD-ROM, programmed memory, such as read-only memory (Firmware), or on a storage medium such as for example, optical or electrical signal carrier. The code (and / or data) for implementing the embodiments of the invention may comprise a source, object or executable code in a conventional programming language (interpreted or compiled), for example C, or assembler code, ASIC (specialized integrated circuit) installation or control code, or FPGA (Field Programmable Gate Array) or code for a hardware description language, such as Verilog (Trademark) or VHDL (Ultra High Speed Integrated Circuit Description Language). One skilled in the art will appreciate that such code and / or data can be distributed among a plurality of related components in communication with each other, for example, distributed over a network.

The invention further provides a data processing system configured to implement embodiments of the above methods.

Thus, in an additional aspect of the invention, there is provided a data processing system for processing measured potential field data from a potential field survey to determine an array of cartographic field parameters when mapping a field, the system comprising: a data storage device for said potential field measured data, said potential data being measured the fields contain data defining a plurality of measurements of a potential field, each with a corresponding measurement position niya and time of measurement; and a program storage device storing processor control code; and a processor associated with said data storage device and said program storage device for downloading and executing said control code, said code comprising a code for controlling a processor for: inputting said measured potential field data; and determining said array of cartographic field parameters using a model comprising a combination of a spatial part representing a spatial variation of said potential field and a time part representing temporal noise in said measured data of a potential field, said code for determining said array of cartographic field parameters being configured to fit the mentioned measured data of the potential field to both mentioned, spatial and temporal, frequent m mentioned model.

The invention additionally provides a carrier transferring the data structure, the data structure containing data defining values of at least A and ρ, wherein A and ρ are determined so that for the array of measurements of the potential field m (x, y, z, t) L (f, m) is minimized, where L (f, m) represents the measure of the difference between f and m, where

,

,

where ρ and λ are the vectors of the spatial and temporal model parameters, respectively, A and B are the matrices, and moreover, the data mentioned carried by the data structure are compared with the mentioned array of the mentioned measurements m (x, y, z, t) or include them .

The invention further provides a carrier transferring a data structure defining a map obtained from data defining values of at least A and ρ, wherein A and ρ are determined such that for the array of measurements of the potential field m (x, y, z, t) the value of L (f, m) is minimized, where L (f, m) represents the measure of the difference between f and m, where

,

,

where ρ and λ are the vectors of the spatial and temporal model parameters, respectively, A and B are the matrices, and moreover, the data mentioned carried by the data structure are compared with the mentioned array of the mentioned measurements m (x, y, z, t) or include them .

The invention further provides a storage medium, such as described above, for transferring an array of cartographic field parameters determined by a previously described method.

Brief Description of the Drawings

These and other objects of the invention are described below by way of example only in connection with the accompanying drawings, in which:

figure 1 shows an example of a temporal interpolation function for a joint model of an equivalent source and temporal drift in accordance with an embodiment of the invention;

figure 2 shows a synthetic model and a simulated picture of the flight;

3a-3c show a comparison of the predicted G _zz from inversions of the equivalent source using synthetic measurements: a) without drift simulation (spatial equivalent source only), b) while drift modeling (augmented model) in accordance with an embodiment of the invention, and c) without synthetic noise measurement added to the model (true signal - no noise);

figure 4 shows an aircraft with aerial survey data and an example of a data processing system configured to implement an embodiment of the method in accordance with the invention; and

5 is a flowchart of a procedure for processing measured potential field data for implementing an embodiment of a method in accordance with the invention.

Detailed Description of Preferred Embodiments

The following describes techniques that, in embodiments, use joint simulation of an equivalent source and time drift for measured potential field data.

Where it is a field, in particular a gravitational field, it is not limited to a vector field, but includes a scalar field and a tensor field, a potential field and any derivatives obtained from a potential field.

Potential field data includes, but is not limited to, gravimetric data, gravitational gradiometric data, vector magnetometric data, and true magnetic gradiometric data. Elements and representations of a potential field can be obtained from a scalar quantity.

For gravity, the corresponding potential is the gravitational scalar potential Ф (r), defined as

,

,

where r, ρ (r '), G - respectively, the measurement position of the gravitational field, mass density at the location r' and the gravitational constant. Gravitational acceleration, which is a measure of the gravitational field, is the spatial derivative of the scalar potential. Gravity acceleration is a vector that sets its direction. It is represented by three components with respect to any chosen Cartesian coordinate system:

.

.

Each of these three components changes in each of the three directions, and nine quantities thus form the tensor of the gravitational gradient:

.

.

The fundamental equations and relations for potential fields follow from an analysis of the properties of a scalar potential function, its derivatives, its Fourier transforms, and other mathematical quantities.

It follows from one of Green's theorems that if any spatial derivatives of the scalar potential (including the scalar potential itself) over a closed surface are known, then the value of the spatial derivative is known at all points within the volume bounded by this surface. The consequence of this is that as soon as this quantity is known at all points, by differentiation and integration all other derivatives of the scalar potential can be obtained, including the scalar potential itself. Thus, the scalar potential and all its derivatives are actually known at all points within the volume if only one of its derivatives is known over the surface enclosing this volume. This means that a complete 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. It also follows from this that, at least in theory, it does not matter which quantity is measured, in this case the apparatus that measures the desired value with the highest signal-to-noise ratio is simply selected.

Differentiation of the above gravitational scalar potential with the need gives:

,

,

which, regardless of areas, goes over to the Laplace equation - an important fundamental relation in the theory of gravity:

.

.

Harmonic functions satisfy the Laplace equation, and they have many properties that can be used to analyze the data collected during surveys of a potential field.

Data can be analyzed and processed using a number of techniques that work with data collected during the survey as a starting point, but which then change the data and / or their format so that all values correlated with the measured values appear on the correct 2D grid, which is located on a horizontal, with a fixed height of the analytical plane (leveling and gridding).

In gridding, generally speaking, the surveyed area is divided into rectangular cells, the sides of which are preferably aligned with the main flight directions for shooting, and then the actual measurement data is replaced by data (coordinate-referenced data), which are "equivalent" to the measured data, but which are now marked values at the points in the middle of each cell. The size of each cell can be selected based on the average separation of the flight lines in two orthogonal directions. Once the data is placed in this coordinate format, it is much easier to mathematically process it. Data can be processed as a set of numbers, for example, by statistical or other methods, providing the best estimate of the potential field on the horizontal analytical plane.

The authors tentatively described an improved equivalent source method (see PCT / GB2006 / 050211, incorporated herein by reference).

Equivalent Source Methods

To understand the embodiments of the invention, it is useful to describe some aspects of modeling methods for an equivalent source.

It is known that the gravitational field outside the body can be modeled as coming from a substance located entirely in an infinitely thin layer at the surface of the body and which exactly repeats the shape of the body surface. Such a layer defines a two-dimensional equivalent source - that is, a source of gravity that produces essentially (theoretically accurate) the same gravitational features as the body itself. There are many options for determining equivalent sources: they can be as described above, or they can be strictly horizontal, they can be completely or partially above or below the surface of the earth; they can be two- or three-dimensional. However, they have the same property that they usually form essentially the same gravitational field as that formed by the removed plot of land. For further information, reference is made to the book by R.J. Blakely, "Potential Theory in Gravity and Magnetic Applications", Cambridge University Press, 1995.

When shooting, you can measure G _zz as a function of position, r _measure , using a gravitational gradiometer and working with it without the need to create other elements of the gravitational gradient tensor. This can be used to create a display of the main mass distribution. The aforementioned Green theorem shows that the main mass distribution can be obtained, in principle, from G _zz , although the later Green's theorem, which is not required to be used explicitly, is described in the embodiments of the techniques.

In an example of an equivalent source method, the surface of the survey zone is divided into small parts, usually with a side of the order of 50 m, which can be considered plates or mass elements. The resolution of the model can be chosen to roughly match the resolution of the measurements. It is easy to directly calculate the gravity of each plate (see, for example, R.J. Blakely, "Potential Theory in Gravity and Magnetic Applications", Cambridge University Press, 1995), the mass of which is adjusted until the best complete fit to the measured data is obtained. This mass definition can use the standard least squares fit procedure. A fit can be obtained by matching the data in the true measurement positions with the gravitational field created by the proposed equivalent source in identical true measurement positions. This process is mathematically rigorous and does not require any artificial adjustments to the data to fit the horizontal rectangular survey.

Once fitting is obtained, the model is considered the primary data array. All subsequent analysis to determine the geological structure may then (but not necessarily) then consist in comparing and minimizing the difference between the gravitational field that any given geological structure can create and the gravitational field of an equivalent source. One significant advantage of this technique is that the best fit arises from the mass distribution, albeit synthetic, and therefore the best fit solution will automatically satisfy the Laplace equation. This represents an improvement in the method that gives a numerical best fit, but which does not impose an additional restriction on the satisfaction of the data with the Laplace equation, that is, it can result from a real mass distribution.

The equivalent source method does not use a surface corresponding to topography, it can use sources that cover any surface that can be constant in height or not, above or below the true surface of the earth, which can cut the true surface of the earth or not, etc. Additionally, the equivalent source method does not use a surface, but can use a three-dimensional model. Similarly, the plates can be of any size or shape in two or three dimensions (or even can be point-like), and they do not even have to correspond in size or shape; indeed, this contributes to the effectiveness of the analysis if their size and geometry can vary depending on how quickly topography and / or geology can vary in each area.

An equivalent source following the topography or even basic geology, as far as it can be reasonably marked, is likely to produce less variation in the mass of individual plates, but the final result, in principle, will not be significantly affected by any reasonable choice. Thus, for example, if a specific geological layer or, for example, a geological anomaly, for example a kimberlite pipe, is suspected, an equivalent source can be determined with this in mind. Likewise, more sources can be added in areas where there is a more detailed or faster change in geological structure. Thus, in some respects, the choice of forms of equivalent sources can be conceived as similar to the discretization of a technical product in finite element analysis.

In part, the mathematical complexity of the process is determined by the number of sources used and how many of these sources are used in the analysis in each position when shooting. The advantage of the technique is that for some components of gravity or the gravitational gradient, it is possible to use direct sources in the region of the measuring point, which can significantly reduce the complexity of the analysis. The threshold distance is usually several kilometers, for example between 1 and 10 kilometers, although this depends to some extent on geography (for example, it may be necessary to extend the distance to cover a large nearby mountain). For example, the source element on one side of the survey area is unlikely to make a significant contribution to the measured field on the other side of the survey region, and where the influence of the source element is much less than the expected or actual noise, such contributions can be ignored (set to zero). This is useful since such a matrix may contain, for example, 100 K per 500 K matrix elements, and processing difficulties are significantly reduced if they are mostly zeros.

An equation is defined below in which the contributions of the source elements to the measured signal depend on the matrix values of the source influence, and where (as is usually the case in practice) the influence of many sources can be ignored, this matrix becomes a sparse matrix, allowing more efficient numerical processing.

As soon as an equivalent source is created, it becomes possible to predict any derivative of the scalar gravitational potential on the surface through direct direct calculation. This process is useful for both analysis and perspective visualization. Thus, as soon as an array of mass elements is found, direct calculation (that is, summing up the influence of mass elements) allows one to obtain other components of the scalar potential and, by differentiation, other components of G. The obtained values of G can be compared with a geological model (denoted as “interpretation”) to determine basic geological structure.

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

.

.

In the above equation, the function F is called the Green function (Blakely, ibid., P. 185, incorporated by reference), and r _mass-element determines the location of the mass element (for example, the center of gravity or some other specific point).

Functions F are standard functions well known to those skilled in the art (and their applications are not limited to gravimetric surveys). This, in essence, is the influence of the source (mass element) of a unit mass or density and a certain shape, available at a certain (measuring) point. The source may be a point source, a sphere, or an ellipsoid, but in practical applications of the technique, the source is often a prism, which may be incorrect. Many textbooks list Green's functions for simple forms; functions for more complex source geometries can be found in the literature. In addition, the influence of the source is superimposed so that if a complex form can be discretized into many simpler forms, then the Green's functions for discrete forms can be summed. In principle, this allows one to determine the numerical values of the Green's function for any arbitrary form, although relatively simple forms are usually preferable. As an example, the Green's function F for a rectangular prism is defined in the same place in Blakely, p. 187, and in this case it is included by reference; there are 8 members, each of which corresponds to the top of the prism.

The equivalent source inversion method can combine multiple measurements of the potential field distribution (the observed values of the potential field are actually spatial derivatives) into a single model. The method includes finding the parameters of the source model (for example, density or geometry) in such a way that the data calculated from the model are adjusted to this array of measurements. The inverted equivalent source model can then reconstruct the signal during measurements and, within certain limits, can be used to re-design the data for different locations. In general, the number of adjustable model parameters is less than the number of independent measurements, and thus, in accordance with the law of large numbers, the signals recalculated from the model will have greater signal-to-noise ratios than in the initial measurements.

The design of an equivalent source is larger than the measurement of the averaging method, since it displays a convenient and real way of compiling an array of spatial basis functions that support the functional form of the main measuring signal. Basic functions may allow the union of various types of potential field measurements, for example; gravitational measurements and measurements of the gravitational gradient tensor can all be inverted to a single model. The basic functions can also contribute to the greater stability of the inverted model relative to some false variations in the measurements, since if they do not appear as true variations of the potential field, then they are not supported by the functional form of the model. An equivalent source model can further focus its functional space by eliminating solutions that are known as physically invalid. For example, an equivalent source model placed near or below ground level will naturally filter out non-physical high-frequency variations in measurements (although the model can, in principle, be placed anywhere). Otherwise, such variations can create large errors in a less limited model.

An increase in the ratio of the number of measurements to the number of degrees of freedom of the model improves the inversion accuracy, making the resulting model of an equivalent source less sensitive to measuring noise and interference. However, in many cases, due to practical difficulties associated with collecting a complete sample of variations of the potential field, this ratio often becomes small, and the inversion of the equivalent source can become difficult.

Low-frequency noise (drift) in aerial survey measurements worsens the inversion process, since these phenomena may resemble a true signal (although the term drift usually refers to random monotonic variation, in this description it is used to display any kind of correlated noise with a characteristic frequency much lower than the dominant frequency and signals of interest). Such a change can therefore easily be accounted for by the equivalent source model, especially in places where it is not defined. The consequences of successive direct calculations of the equivalent source then show false signs that correlate with the measurement pattern. In addition, the coherent noise components are not averaged during the inversion process and therefore lead to a worse estimate of the model parameter.

An addition to the equivalent source technique is described below, in which a separate model is introduced to compensate for the correlated temporal noise. By using a separate model that works simultaneously with the equivalent source model, it is possible to make the estimate of the source parameter less susceptible to this type of noise during measurements. After the inversion, the augmented model can separately predict the desired spatial distribution of the potential field and the long-wave time noise accompanying the measurements.

Joint simulation of equivalent source and time drift

The measured value m during geophysical surveys, corresponding to the location (x, y, z) and time t, can be divided into a number of terms

,

,

where S (x, y, z) represents the signal, I (t) is the source of interference, C (t) is the correlated noise, and n is purely random noise. The main goal of post-acquisition processing is to determine with the best accuracy the true signal S (x, y, z). Interference sources I (t) are corrected by suitable measurements of the main disturbances and corresponding transfer functions with an error. In theory, after these corrections, the only remaining unwanted terms are the correlated noise C (t) and white noise n. White noise is fundamental, and its influence can only be reduced by the law of large numbers in multiple measurements. Most of the power of the correlated noise can, however, be modeled by an interpolation function with a time period that is similar to or shorter than the characteristic time period of the noise.

Figure 1 shows a simple temporal interpolation function formed by stitching piecewise-linear sections between them at equal intervals. More specifically, FIG. 1 shows an exponentially correlated noise with a characteristic time of 500 seconds modeled by a piecewise linear interpolator connecting nodal points (squares) every 400 seconds.

The interpolator is continuous within certain limits, but has discontinuous derivatives at nodal points. The model is fully defined by means of the values at the nodal points and therefore makes the number of degrees of freedom equal to the number of nodal points. In this case, the nodal values become model fitting parameters during optimization (see below). If more complex interpolation functions are used, for example, high-order polynomials or splines, then the number of degrees of freedom can be large. As can be seen from figure 1, the long-wavelength components of the correlated noise can be quite modeled by the main interpolator with a time period slightly less than that for noise. The described simultaneous model of the equivalent source and time drift allows the low order drift model to absorb most of the temporal correlated noise within the measurement array, leaving the model of the equivalent source only to establish correspondence with the spatially correlated signal. The remaining unmodeled portion of the measurements (residual fit) is then largely white noise, indicating that the inversion statistics are compatible with a good estimate of the model parameter.

In the case of a linear inversion model, the augmented model can be written as

,

,

where f is the model vector of the direct calculation of measurements, ρ is the vector of model parameters of the equivalent source (for example, density or discrete modules-mass sources), A is the matrix that connects the responses of the source elements to the measurement locations, λ is the vector of parameters of the drift model (for example, nodal values in a piecewise linear model), and B is the matrix that relates the interpolation of the drift model to the moments of measurement. In this form, with a given array of measurements m (x, y, z, t), the model parameters ρ and λ can be determined using any standard optimization technique that minimizes the remainder when fitting

,

,

where L represents the measure of the remainder; the L ₂ norm, for example, performs least squares optimization. Minimizing the squared error, the norm of L ₂ , has the advantage of being easy to calculate, and there are many commercially available math / library components to implement the least squares fit. Another option is to use a module instead (norm L ₁ ). This can lead to a possible omission of data point outliers, which makes optimization more time consuming; in this case, suitable mathematical components of packages / libraries are commercially available, but the computational cost may be too large (since a more complex mathematical apparatus is required).

After complete inversion, a member of the equivalent source is used in isolation or to evaluate the initial signal underlying the measurements,

,

,

or to predict different values at possible different locations. Similarly, the second term in (2) can be used in isolation to estimate the drift in measurements

The model of equation (2) can be supplemented for joint modeling of signals from multiple data sources of the measured potential field. In this case, the general model of the equivalent source can be used for different data sources of the measured potential field, but preferably members of the individual noise (drift) are included in the model for each different source of measured data. Different sources of measured potential field data can be, for example, different devices for measuring various components of G (G _xx , G _yy , G _zz , G _xy , G _xz , G _yz - thus, six different devices of this type are possible) and / or different measuring devices g (g _x , g _y , g _z - thus, three devices of this type are possible). For example, one or more data components of the gravitational and gravitational gradient potential fields can be jointly measured in this way with multiple measuring devices (although they can be combined in a single measuring device). Adjustable minimization can be performed in the same way as noted above (and discussed further below), although with a large number of determined drift parameters.

At this stage, it is useful to define in more detail the matrices A and B and the vectors ρ and λ. Continuing to consider the spatial part of the model, the initial stage is usually to build a model that follows the landscape in the shooting area. The model is discretized into a series of finite elements, each of which determines a mass element, the density of which is determined by the inverse solution. The matrix element A _ij displays a model, directly calculated signal, giving a contribution from the mass element j in the measurement position i, when the unit density is attributed to the mass element. The type of direct calculation used corresponds to the type of measurement signal, for example; if there are n measurements and m were taken in the same position, but n measurements are measurements of the gravitational gradient G _zz , and m measurements are measurements of the gravitational gradient G _xx , then the corresponding direct calculations from the mass element j, A _{n, j} and A _{m , j} will correspond to G _zz and G _xx, respectively.

The vector ρ defines an array of actual densities of the mass elements that determine the model, and after optimization contains a part of the solution of the equivalent source. Thus, the vector ρ actually provides an array of parameters representing the field. As mentioned above, the product Aρ is used to predict the actual signal during measurements.

To facilitate optimization (2), approximations can be introduced when constructing matrix A. This can take the form of truncation when negligible small matrix elements are replaced by zeros), or averaging, when the effects of several individual elements are combined into a single element that displays the average effect. Both of the above approximations become valid when the distance between the measurement point and the source element becomes large compared to the size of the source. Using such approximations, the spatial matrix A becomes sparse, allowing more efficient numerical matrix calculations.

The elements of matrix B determine the equations that relate the predicted drift in the measurements to the parameters λ of the drift model. For the case of a piecewise linear interpolation function (referring to FIG. 1), if the ith measurement was performed at time t _i , which falls between the nodal parameters λ _j and λ _{j + 1 of the} drift determined at the (nodal point) time points T _j and T _{j + 1} , then the linear interpolation prediction for the drift in the measurements is specified by

Expressed in accordance with the matrix equation (2), the above equation defines the matrix elements B _{i, j} and B _{i, j + 1} as follows:

Times T _{j of} nodal drift parameters determine the resolution of the drift model. When there is a gap in the array of measurement data, for example, at the end of the survey line, or when the recording of data is suspended for some reason, it is preferable to restart the drift model by forcing the drift parameter at the point in time at which measurements are resumed.

In the case of multiple measurement types in the same survey (for example, when there are several instrument output channels) it is convenient to divide matrix equation (2) so that each measurement type occupies a block of rows in matrices A and B. Thus, the first group of drift parameters in λ displays the simulated drift for the first type of measurements, the second array displays the second type of measurements, and so on. It is important to note that although each type of measurement has its own drift model, there is only one equivalent source model.

If the remainder spectrum (f-m) is not sufficiently white in the time domain, then this indicates that the resolution of the drift model may be inadequate, and a shorter time period or a higher order interpolator may be beneficial. Similarly, if the remainder is represented in the spatial domain, this may indicate that the equivalent source model is too limited. If at the same time, the remainder still shows a noticeable correlation, then a more fundamental problem is possible, for example, excessive interference in measurements or a poorly designed model of an equivalent source, for example.

Since the nature of the signal and the drift is different, one is correlated in the spatial domain and the other is correlated in the time domain, the separation of these quantities into two terms on the right-hand side of (2) seems natural during optimization. To illustrate this, consider two adjacent aerial survey lines. Although the drift is strongly correlated along individual lines, when viewed in the spatial domain there is no correlation between the lines - adjacent lines could be obtained with an interval of weeks, or even with different instruments. In optimization (3), since a variation of this type is not well approximated by an equivalent source model, it can be correctly described by a member of the drift model. Conversely, a true signal usually appears outside the range of the drift model and is therefore preferably described by a member of the equivalent source model.

At the same time, there are many situations, especially in the indefinite case, when it is difficult to distinguish between the low-frequency noise of the device and the true long-wave signal. If it happens that drifts along adjacent lines have the same general tendency, the noise may seem spatially correlated and during the inversion can enter and distort a member of the equivalent source. Conversely, a long wavelength of a geological signal can easily occur as a drift in time, therefore, self-discovery refers to the wrong term in (2). The result is a loss of signal intensity in predictions (4) of the equivalent source associated with undesired noise trends.

To add a control level on this problem, one can introduce a separate regularization for two parts of the extended model (2). Gradient regularization is preferable for the density distribution ρ of the equivalent source, because, while the geometry of the model resembles something physically plausible, real density distributions tend to have smooth variations. It is preferable to choose smoother density solutions that therefore tend to improve the inversion integrity. For drift model parameters, λ, Tikhonov standard regularization is the preferred choice. With regularization, expression (3) takes the form

,

,

where a is the coefficient of full regularization, and b is the coefficient that controls the relative regularization between the two models (in practice, the terms of the regularization are normalized by model matrices). Increasing b makes the equivalent source model more meaningful by attributing ambiguous variation in measurements, preferably density variation rather than drift. There is no definite way to get the right values for the regularization coefficients, but L-curve analysis and other empirical methods can give direction to determine the best regularization trade-off between fitting measurements and solution behavior. Ultimately, a sequence of parameter tests is required to determine the best values for a and b.

Thus, generally speaking, in some preferred embodiments, a regularization is introduced, which in the spatial part of the model is designed to find the smoothest distribution fit for measurements (the expected large changes in the underground density are preferably taken into account, where possible, by constructing the initial model of the source element). In a broad sense, the value of a determines to what extent a smooth solution with the model should be achieved, and in practice this parameter can vary in accordance with the expected basic geology and to achieve what turns out to be a good result. In addition, the selection of a large value for b cancels the selection of members of the drift model (since the smaller members of the drift model are then necessary to completely minimize equation (8)).

However, the most effective way to resolve the ambiguity between drift and signal is to ensure that there are enough measurements when shooting, which are close in space but well separated in time. Consider, for example, points in a survey data array where two lines (often perpendicular) intersect. At these locations, the signal S (x, y, z) is common to two measurements on each line. (A pair of measurements can be considered a correct intersection point if their separation is sufficiently small compared to the correlation width of the signal commensurate with the height of the survey.) The difference between the measurements taken at time t ₁ and t ₂ should therefore be due to the correlated noise, I (t ₂ ) -I (t ₁ ) plus some random noise. In matrix construction (2), these pairs of points lead to duplicate rows in the equivalent source term, but independently for rows in the drift model term. Therefore, they together provide equations that only substantially affect the drift solution, and therefore can resolve the ambiguity between the long wavelength signal and the drift. By taking into account the spatial variation of the signal, simultaneously solving the equivalent source model and the drift model is more accurate than leveling intersection, since the drift model is essentially derived from measurements in the survey at multiple positions on different survey lines, potentially all measurements, and not just those for intersection points.

A synthetic example is presented to demonstrate the characteristics of the combined equivalent source model and drift model. The synthetic model contains a selection of 5 kimberlite pipes with various sizes taken by an onboard gravity gradiometer in accordance with the regular aerial survey, as shown in FIG. 2, which illustrates the synthetic model and simulated flight picture.

In modeling the component of the gravitational gradient, G _zz was directly calculated for the location of the survey picture to display the signal measured by the gravitational gradiometer. An exponentially correlated temporal noise was then added to the synthetic signal to display drift and other low-frequency noise that would be present in a real measurement system. Simulated measurements were used to invert the model of the equivalent source, first, without and then with the addition of the drift model. In both cases, the inverted model was used to directly calculate G _zz for an array of coordinate points having the same height as in the original survey. Figure 3 shows a comparison of the predicted value of G _zz from inversions of the equivalent source using synthetic measurements: a) without, b) while simulating drift, c) the perfect response achieved without synthetic measurement noise added to the model.

In the above example, to illustrate the current method, no additional data processing procedures for noise reduction were performed, although in practice conventional noise reduction techniques, such as filtering / limiting, can (and often preferably should) also be used. The results from the extended equivalent source model (b) clearly show the level of reliability and protection relative to the correlated noise in the measurements. Without this extension, results (a) show significant artifacts correlated with the survey picture and an almost complete masking of the true signal.

Thus, the addition of a time-drift model of the equivalent source inversion method and the subsequent solution for both models at the same time provides a method for reliable processing of potential field data arrays that contain significant low-frequency noise. Although low-frequency noise can be reduced by conventional leveling procedures, these methods do not work, for example, when shooting has insufficiently intersecting points. The advantage of associating a drift model with an equivalent source model is that the drift behavior can be inferred from a larger number or essentially from all measurements when shooting.

Figure 4 shows an example of an aircraft 10 for capturing a potential field to obtain data for processing in accordance with the method described above. The aircraft 10 contains an inertial platform 12 on which a gravity gradiometer 14 (and / or vector magnetometer) is installed, which provides potential field data to a data collection system 16. The inertial platform 12 is equipped with an inertial measurement module (IMU) 18, which also provides data to the data acquisition system 16, typically containing spatial position data (e.g., pitch angle, roll and course deviation), angular velocity and angular acceleration data and aircraft acceleration data. The aircraft is also equipped with a differential GPS system 20 and 22 LIDAR or the like to provide data on the height of the aircraft above the main terrain. Position and time data is preferably obtained from (D) GPS, optionally in combination with IMU for accuracy.

The aircraft 10 may also be equipped with other equipment 24, for example a magnetometer, a TDEM (Temporary Electromagnetic System) system and / or hyperspectral imaging system, which also feed the data acquisition system. The data acquisition system 16 also has input from the general equipment 26 of the aircraft, which may include, for example, an altimeter, data of airspeed and / or speed relative to the ground, etc. The data acquisition system 16 may provide some initial data preprocessing, for example, to correct LIDARa data for aircraft movement and / or to combine data from IMU 18 and DGPS 20. The data acquisition system 16 may be provided with a communication line 16a and / or non-volatile device 16b of a memory so as to be able to store the collected potential field data and position data for subsequent processing. A network interface (not shown) can also be provided.

Data processing for creating cartographic data for surveying a potential field is generally (but not necessarily) autonomous, sometimes in another country where the survey data were collected. It can be seen that the data processing system 50 comprises a processor 52 connected to a code and data storage device 54 and an input / output system 56 (for example, containing interfaces for a network and / or storage media and / or other communication) with a user interface 58 , for example, containing a keyboard and / or mouse. The code and / or data stored in the memory device 54 may be provided on a removable storage medium 60. In operation, the data includes data collected from a potential field survey, and the code contains a code for processing this data to generate map data in the embodiments in accordance with the procedure shown in FIG. 5, described below.

FIG. 5 shows an example of a procedure performed by a data processor, which in embodiments may comprise a universal computer system for processing data from a flight survey in accordance with previously described techniques. Thus, in step S200, the procedure enters the measured data of the potential field, for example, from a gravitational gradiometer and the corresponding 3D position data. If necessary, in step S200a, some pre-processing may be applied, for example, to remove anomalies and / or to reduce (or increase) or select the data to be processed.

In step S202, the procedure fills in the matrices A and B as described above, and then creates a sparse matrix equation related to the measured potential field data, for example G _zz , the measured 3D position and temporal data, spatial cartographic parameters (ρ) of the field and temporal drift parameters (λ). An equivalent source terrain model (S202a) may be specially constructed, or a standard regular grid may be used. The procedure then solves the matrix equation (S206) to determine the array of cartographic parameters of the field, more specifically the array of forces of elements of equivalent sources, and they can provide output from the procedure as cartographic data. Alternatively, the map can be directly calculated as part of the procedure (S208) by performing direct calculation using these parameters. For example, a two-dimensional map of the type shown in FIG. 3 or a three-dimensional map (a more accurate geological model) of the type shown in FIG. 2 can be created.

Without a doubt, many other effective alternatives will be apparent to those skilled in the art. For example, the described technique can be carried out in Fourier representation, for example, along the lines described in the considered PCT patent application in the same place (PCT / GB2006 / 050211, incorporated by reference).

Although a technique is described herein using a preferred example of aerial survey of a potential field, embodiments may also be used for marine surveys of a potential field performed on a ship and, in general, for surveys of a potential field performed on other moving platforms or vehicles.

The techniques are not limited to the processing of gravitational data, but can also be used, for example, when processing magnetic field data. The measured potential field data can thus be obtained by measuring the magnetic field and / or the flux density vector and / or its magnitude, for example, using a measurement made with a magnetic gradiometer. Elements of an equivalent source may then have, for example, surface current density or field strength.

It should be understood that the invention is not limited to the described embodiments and encompasses modifications obvious to those skilled in the art that fall within the spirit and scope of the appended claims.

Claims (27)

1. A computer-implemented method for processing measured potential field data obtained from an aerial or marine survey of a potential field to determine an array of cartographic field parameters for a model for mapping a field, the method comprising:
introducing into the processor the measured potential field data from an aerial or marine survey of the potential field, said measured potential field data containing data defining a plurality of potential field measurements, each with a corresponding measurement position and measurement time; and
determining, using a processor, said array of cartographic field parameters by fitting said model to said measured potential field data, the model comprising a combination of a spatial part representing a spatial variation of said potential field and a temporal part displaying temporal noise in said potential field measured data. 2. The method according to claim 1, wherein said spatial and temporal parts of said model comprise spatial and temporal terms of a model equation, and said definition comprising jointly determined estimates of said spatial and temporal terms for determining said cartographic field parameters. 3. The method according to claim 1, in which said model has the form:
$$f(x,y,z,t)=A\rho (x,y,z)+B\lambda (t),$$

where ρ and λ are the vectors of the spatial and temporal model parameters, respectively, A and B are matrices, and f contains the vector of direct calculation of the measurements estimated from the model. 4. The method according to claim 1, wherein said spatial portion of said model comprises an equivalent source model. 5. The method according to claim 1, wherein said time portion of said model comprises an interpolation function for interpolating between nodal points in time intervals greater than 60 s. 6. The method according to claim 1, wherein said temporary part of said model comprises a piecewise linear model. 7. The method according to claim 1, wherein said noise comprises a measured drift. 8. The method according to claim 1, wherein said definition comprises minimizing a measure of difference between said measured potential field data and data predicted by said array of cartographic field parameters to determine said array of cartographic field parameters. 9. The method according to claim 1, wherein said definition comprises regularizing said spatial part and said time part of said model. 10. The method according to claim 1, wherein said definition comprises regularizing said spatial part of said model with a first regularization and regularizing said temporal part of said model with a second, different regularization. 11. The method according to claim 9, containing the use of gradient regularization for said spatial part of said model. 12. The method according to claim 9, containing the use of Tikhonov regularization for said time part of said model. 13. The method according to claim 1, wherein said fitting comprises co-fitting data from potential field measurements from a plurality of different survey lines of said potential field survey. 14. The method according to claim 1, in which the said measured data of the potential field contain one or more gravimetric data and gravitational gradiometric data. 15. The method according to claim 1, in which the aforementioned cartographic field parameters contain values for elements of a magnetic equivalent source. 16. The method according to claim 1, wherein said cartographic field parameters contain values for elements of an equivalent mass source. 17. The method according to claim 1, wherein said potential field measured data contains data from a plurality of potential field measurement sources, said spatial portion of said model being common to said plurality of measurement sources, and wherein a separate said temporary part of said model is provided for each said source of potential field measurement. 18. The method according to claim 1, additionally containing a map definition using said cartographic field parameters. 19. A medium carrying a code controlling a processor for executing when starting the method according to any one of the preceding paragraphs. 20. A data processing system for processing the measured data of a potential field by surveying a potential field to determine an array of cartographic parameters of a field for mapping a field, the system comprising:
a data memory for said potential field data being measured, said potential field data being measured comprising data defining a plurality of potential field measurements, each with a corresponding measurement position and a measurement time; and
program memory storing processor control code; and
a processor associated with said data memory and said program memory for downloading and executing said control code, said code comprising a code to control a processor for: inputting said measured potential field data; and
determining said array of cartographic field parameters using a model comprising a combination of a spatial part representing a spatial variation of said potential field and a time part representing temporal noise in said measured data of a potential field, said code for determining said array of cartographic field parameters being configured to fit said measured data of the potential field and to the aforementioned spatial parts, and to the mentioned temporal m part of the said model. 21. A computer-implemented method of an equivalent source for processing measured potential field data for determining cartographic field data for a model in which a potential field is modeled as a combination of a spatial signal and a temporary noise component, the method comprising
measurement of potential field data using aerial or marine surveys of potential fields;
determining with a processor the joint best approximation of the measured potential data of the potential field to the simulated combination of the spatial signal and the temporal noise component to determine the cartographic field data for the model. 22. The method according to item 21, further comprising including one or both spatial regularization and temporal regularization in said simulated combination of said spatial signal and said temporary noise. 23. The method according to item 21, in which the said measured data of the potential field contain gravitational data. 24. A medium that transfers code that controls the processor, which, when launched, implements the method according to claims 21, 22, or 23. 25. A data processing system including a medium according to claim 24. 26. A medium that transfers an array of cartographic field parameters determined by the method according to any one of claims 1-17. 27. A medium that transfers cartographic data determined by the method according to any one of paragraphs.21-23.

Images