"Method and apparatus for detection using magnetic gradient tensor"

Cross-Reference to Related Applications

The present application claims priority from Australian Provisional Patent Application No 2007903249 filed on 18 June 2007, the content of which is incorporated herein by reference.

Technical Field

The present invention relates to detection of magnetic bodies, and in particular relates to defining the location and magnetic moment vector of a detected magnetic body by analysing the complete magnetic gradient tensor.

Background of the Invention

Airborne, ground, marine and downhole magnetic measurements have played an important role in detection of magnetic targets such as land mines, naval mines, submarines, shipwrecks, unexploded ordnance (UXO), archaeological artefacts and structures, buried drums containing, for example, toxic waste, ferrous particles in food and hospital laundry, and many others. Magnetic surveys also have major applications to mapping geology, in exploration for minerals and energy resources, and in environmental and archaeological surveys. Total magnetic intensity (TMI) sensors may be used in such surveys. However, limiting factors when using TMI sensors include complicating effects of the orientation of the geomagnetic field, which produces asymmetric anomalies that are displaced from source locations, the sensitivity of the shape and strength of TMI anomalies to direction of source magnetisation, relatively low resolution compared to magnetic field gradients, frequent need for post- processing techniques, and the fact that only a limited, distorted and fuzzy view of geology or terrain is obtained due to cost constraints on sampling density. Moreover, it is not possible to locate targets _, , uniquely between survey lines or to delineate adversely oriented structures in low latitudes near the Earth's equator.

Magnetic surveys using vector magnetometers such as fluxgate-based or SQUID-based detectors have also been considered, because of the extra information that is in principle obtainable from measuring three components of the field, rather than a single

TMI component. Vector surveys, where the direct measurement of vector components has been attempted, have met with mixed success. The accuracy of direct measurement of the field vector is largely governed by orientation errors, which are generally so large that the theoretical derivation of the components from sufficiently densely sampled TMI is actually preferable. Vector magnetometry is so sensitive to orientation errors because the anomalies that need to be detected are usually very weak compared to the background geomagnetic field. For example a 1° change in orientation produces changes in vector components of the geomagnetic field of up to approximately 1000 nT, whereas anomalies are often of the order of a few nT, or even sub nanoTesla.

An important application of magnetic sensors is the detection, location and classification (DLC) of magnetic objects. Compact magnetic bodies can be well represented by a point dipole source, except very close to the body. A number of methods have been proposed for locating dipole targets from magnetic gradient tensor data. Some of these methods are based on point-by-point analysis of the eigenvectors of the tensor, and are adversely affected by noise in individual measurements of the gradient tensor elements. Calculated eigenvectors are inherently sensitive to noise in gradient tensor elements and eigenvector-based methods may not always be robust. If two eigenvalues are almost equal, for example, small perturbations of the eigenvalues can produce approximately 90° jumps in the corresponding eigenvector orientations. Methods that rely on recalculation of tensor elements with respect to principal axes of the tensor (which are defined by the eigenvectors) also suffer from this problem. Moreover, as the analysis is essentially point-by-point, the solutions are unreliable for individually noisy tensor measurements.

Furthermore, there is an inherent four-fold ambiguity in obtaining solutions for dipole location and orientation of its moment from point-by-point analysis of gradient tensors, which must be resolved by comparing solutions from different sensor locations, rejecting those that are not consistent (the so-called "ghost" solutions) and retaining the solutions that exhibit the best clustering. Existing methods of dipole tracking are also not robust to the contamination of the measured signal by variable background

gradients, interfering anomalies, instrument drift or departures of the target from a pure dipole source.

Proposed methods for resolving ambiguity of gradient tensor measurements by eliminating ghost solutions rely on numerical differentiation of gradient tensor time- series to obtain gradient rate tensors, or by using arrays of triaxial magnetometers that constitute two or more tensor gradiometers with small spatial separation. Both approaches effectively amount to incorporating second order gradients of magnetic components into the analysis, which are more subject to noise than the first order gradient tensor. Furthermore, the incorporation of higher order gradient information greatly increases the mathematical complexity and computational demands of the solutions, making near real-time tracking difficult. Moreover, numerical differentiation amplifies noise in the original data, and the gradient rate tensor has a 1/r ⁵ fall-off, which makes the signal very weak, except very close to the source.

Any discussion of documents, acts, materials, devices, articles or the like which has been included in the present specification is solely for the purpose of providing a context for the present invention. It is not to be taken as an admission that any or all of these matters form part of the prior art base or were common general knowledge in the field relevant to the present invention as it existed before the priority date of each claim of this application.

Throughout this specification the word "comprise", or variations such as "comprises" or "comprising", will be understood to imply the inclusion of a stated element, integer or step, or group of elements, integers or steps, but not the exclusion of any other element, integer or step, or group of elements, integers or steps.

Summary of the Invention According to a first aspect the present invention provides a method of locating and characterising a magnetised body, the method comprising:

during relative movement between a magnetic gradient tensor sensor and the magnetised body, along a profile, obtaining a plurality of magnetic gradient tensor measurements corresponding to points along the profile, at different relative orientations between the sensor and magnetised body; and determining at least one characteristic of the magnetised body using a rotational invariant calculated from the eigenvalues of the magnetic gradient tensor measurements.

According to a second aspect the present invention provides a device for locating and characterising a magnetised body, the device comprising: a sensor for obtaining a plurality of magnetic gradient tensor measurements at different relative orientations between the sensor and a magnetised body, during relative movement between the sensor and the magnetised body along a profile; and a processor for determining at least one characteristic of the magnetised body using a rotational invariant calculated from the eigenvalues of the magnetic gradient tensor measurements.

The present invention recognises that by measuring the magnetic gradient tensor rather than total magnetic intensity, it is possible to calculate from each such measurement a limited number of solutions for characteristics of the magnetised body. Such characteristics might be the range and direction of the magnetised body relative to the sensor and/or the magnitude and orientation of the magnetic moment of the magnetised body. Moreover, by obtaining measurements at each of multiple relative orientations between sensor and body, it is possible to eliminate spurious solutions and identify a single solution for range and direction information and for the dipole moment vector from such measurements. Thus, the present invention provides for analysis of a number of gradient tensor measurements along a profile to obtain a unique solution, which also averages out noise in the individual measurements, producing more robust and precise estimates of the target parameters. The present invention has the flexibility to use relatively few measurements or a large number, and to analyse a portion of the

gradient tensor anomaly or the complete profile, depending on the signal-to-noise ratio, the desired precision and computational demands.

Additionally or alternatively, the determined characteristics of the magnetised body may in some embodiments of the present invention include an orientation of a magnetic dipole representing the magnetic moment of the magnetised body. For longitudinally magnetised bodies, the dipole alignment can enable identification of the alignment of the body itself.

The gradient of a field in three dimensions is a second rank tensor, and accordingly has nine components. In the case of static magnetic fields, for regions free of electric currents Maxwell's equations impose restrictions which imply that only five of the gradient components are independent. To determine the gradient tensor of a magnetic field it is therefore necessary to measure all five of these components. In preferred embodiments of the invention, each magnetic gradient tensor measurement may be obtained by use of a device of the type set out in International Patent Application No. WO 2004/015435, the content of which is incorporated herein by reference. However, any form of sensor that produces full gradient tensor measurements with sufficient sensitivity to detect targets of interest and an adequate bandwidth to define the entire frequency spectrum of the response due to relative motion of the sensor and magnetic source, can be used in this invention.

Embodiments of the present invention may comprise a single tensor gradiometer that moves relative to the target. Such a sensor may comprise an array of at least three non- collinear triaxial magnetometers, a set of three intrinsic single axis gradiometers rotating about differently oriented axes, arrays of intrinsic axial and/or planar gradiometers with referencing magnetometer(s), or other devices.

The present invention further recognises that obtaining measurements of the gradient tensor is advantageous, because measuring field components directly imposes heavy constraints on sensor orientation, sensor orientation stability, and accuracy of

measurement of sensor orientation, if an anomalous field associated with an arbitrarily oriented magnetised body is to be determined in the presence of the very much larger geomagnetic field.

The present invention also recognises that the gradient tensor of a point dipole has the three mutually orthogonal eigenvectors, one of which must, from symmetry considerations, coincide with the normal to the plane that contains the dipole moment vector m and the sensor-to-dipole displacement vector r. This eigenvector is parallel to the vector, or cross, product m x r, and corresponds to the eigenvalue with the smallest absolute value. The other two eigenvectors lie within the plane containing m and r.

Furthermore, measuring the magnetic gradient tensor allows rotationally invariant parameters to be worked with, reducing the burden of sensor orientation. The eigenvalues of the gradient tensor, and combinations of them, are rotational invariants, which are much less sensitive to the presence of noise in the measured tensor elements than are the eigenvectors.

Preferred embodiments of the present invention further involve a two-step analysis of the gradient tensor, whereby scalar parameters of the source including the magnetic moment magnitude and the distance to the source at each measurement point are estimated from a profile segment of rotational invariants, and these solutions are then used to extract the vector parameters of the source, such as the direction to the source and the orientation of its magnetic moment. This invention further recognizes that specific rotational invariants that are independent of dipole orientation and spherically symmetric about the dipole location (the scaled moment and a parameter proportional to its square root), can be derived simply from the eigenvalues of the tensor and have particular utility for separating the scalar parameters of the DLC problem from the vector parameters, greatly simplifying its solution. This two-step approach also lends itself to iterative improvement of estimates of source parameters, by estimating and removing background interference terms from the measured rotational invariants and tensor elements, then repeating the analysis on the corrected profile data.

In preferred embodiments of the invention, one or more characteristics of the magnetised body are determined by deconvolving the scaled moment of a point dipole representation of the magnetised body. Alternatively, the square root of the scaled moment may be deconvolved for this purpose. Such embodiments recognise that the far field of a magnetised body is dominated by the dipole component and that working with the point dipole representation may thus provide sufficiently accurate measurement results while simplifying data processing. Such embodiments further recognise that the scaled moment of a dipole, and its square root, are useful rotational invariants suitable for such deconvolution, even in the presence of interference from other sources that produce background gradients.

In preferred embodiments of the invention the plurality of measurements are obtained as the sensor is moved along a substantially linear path.

In a first group of embodiments of the invention, deconvolution of the scaled moment profile of such measurements preferably comprises: modelling the scaled moment and background interference terms determined from the measurements; obtaining sufficient measurements to obtain an over-determined system of linear equations using the model; solving the system of linear equations for unknown parameters of the model and for background gradients of various orders; modelling the individual gradient tensor elements and background interference terms, incorporating scalar parameters that have been estimated in the previous steps; and solving the resulting system of linear equations for the remaining unknown parameters of the model and for background gradients of various orders.

In such embodiments the systems of linear equations may be solved directly by a least- squares best-fit pseudo inverse solution, but is preferably solved by a numerically robust method such as QR decomposition or singular value decomposition.

The scaled moment of a dipole source is defined as:

3m

where m is the magnitude of the dipole moment (in appropriate units) and r is the distance from the point dipole to the sensor. In the absence of interference from other sources, the scaled moment can be expressed as the following combination of eigenvalues: μ = _λ ]- λ; - λ,λ ₃ , where the eigenvalues of the gradient tensor in descending order are λι > λ ₂ > λ ₃ .

In the first group of embodiments of the invention the model of the scaled moment, as calculated from the eigenvalues of the gradient tensor, preferably comprises:

where xo is the co-ordinate of the point of closest approach along the measurement profile, S = V(Y ² + h ² ) is the slant distance from the point of closest approach to the dipole, Y is the horizontal offset of the target dipole from the profile, h is the height of the profile above the target level, a is the base level, and b, c are coefficients of linear and quadratic terms that represent interference from other anomalies. In such embodiments, where the number of measurements obtained n is greater than 11, an over-determined system of linear equations arises, allowing solution of the system for the unknown parameters. Preferably, parameters for which more than one estimate arises are used as a check on the solution, or are solved by a parameter-specific system of over-determined equations. Such parameters, once calculated, are preferably used in turn to allow calculation of a solution for remaining parameters, until the moment m itself may be calculated by using the calculated parameters in the model being used.

In embodiments utilising such a model an iterative solution may be determined by first identifying the point of closest approach x ₀ and then repeating all calculations with this point being used as origin.

In embodiments using the above model, at each sensor location the range from the sensor to the magnetised body may be determined by

where the scaled moment μi has been corrected for background fields.

Preferably, complete information for source moment and location is derived from the individual magnetic gradient tensor elements. Such embodiments preferably comprise: modelling the five independent measured gradient tensor elements in a second model based on known parameters of the model of the scaled moment; solving the deconvolution problem of the second model for parameters giving at least one of: interference terms; magnetic moment source relative to the sensor; and sign and magnitude of the lateral displacement of the source from the point of closest approach.

In a second group of embodiments of the invention, deconvolution of the square root of the scaled moment is undertaken. Such embodiments may then solve for parameters of a model of the square root of the scaled moment and may solve for parameters of a model of the individual components of the gradient tensor, in a similar manner as for those embodiments of the first group described above. However, embodiments of the second group may be particularly advantageous in making the deconvolution problem mathematically equivalent to a Werner deconvolution of the TMI anomaly of a thin sheet. Furthermore, numerical stability of the inversion is greater for the second group, because the denominator of the square root of the scaled moment is only the second power of the distance, rather than the fourth power. Such embodiments may thus

require a number of measurements n to be greater than 5 in order to provide for an over-determined system of linear equations for deconvolution.

Brief Description of the Drawings An example of the invention will now be described with reference to the accompanying drawings, in which:

Figure 1 illustrates a method of detecting and locating a magnetised body in accordance with the present invention;

Figures 2a to 2d are a flowchart illustrating the algorithm for deconvolution of the scaled moment data; and

Figure 3 is a flowchart illustrating the algorithm for deconvolution of the square root of (the scaled moment)/3.

Description of the Preferred Embodiments Figure 1 illustrates a method of detecting and locating a subterranean magnetised body 100 in accordance with the present invention. A sensor 1 10 for obtaining measurements of the complete magnetic gradient tensor is mounted in an aeroplane 120 which travels along a substantially straight and horizontal path 130. The path 130 is at an elevation h above the magnetised body 100. The path 130 is not directly above the magnetised body, but has a point of closest approach which is laterally spaced from the magnetised body by a horizontal distance Y.

While illustrated with respect to an airborne sensor in Figure 1, it is to be appreciated that the present invention may be applied with respect to a sensor which is conveyed on land, underground, above water or underwater. Alternative applications may comprise a stationary tensor gradiometer which is used to track a magnetic object that moves along a quasi-linear path past the sensor at approximately constant known speed.

The present embodiment addresses the need for deconvolution along a profile of gradient tensor data arising from a point dipole source, in the possible presence of background gradients. The magnetic gradient tensor arising from the point dipole

source is referred to as an anomaly. The point dipole model is applicable to any sufficiently compact magnetised body. Outside a spherical surface that completely encloses an arbitrary magnetised body, the field can be expressed as a sum of multipole contributions, with a 1/r ³ fall-off for the dipole term, 1/r ⁴ fall-off for the quadrupole, 1/r ⁵ fall-off for the octupole term, and so on. Thus for any compact body with a nonzero dipole moment the far field is always dominated by the dipole term. More specifically, the external fields and gradients due to a uniformly magnetised sphere are identical to that of a point dipole, because the quadrupole and all higher multipole moments vanish identically.

For a uniformly magnetised body, the dipole moment m = JV, where V is the volume, J = (J _x , J _y , J _z ) = J(L, M, N) is the magnetisation vector, J is the magnetisation intensity, and L,M,N are the magnetisation vector's direction cosines with respect to the x,y,z axes. By Poisson's relation the scalar potential U due to a point pole p and the scalar potential V due to a point dipole m are related by:

m.VU

V = - = -m.V

v ^{m r}

where r is the displacement vector from the dipole to the observation point. In this and all subsequent equations, it is assumed that the magnetic moment incorporates a numerical factor, dependent on the system of magnetic units used, that relates calculated magnetic fields and gradients to the magnetization and volume of the source.

For example, if the magnetization is given in AJm and distances are measured in metres, the moment in Am ² should be multiplied by 100 nTm/A (= 100 nH/m) to yield calculated fields in nT and gradients in nT/m.

The gradient tensor due to a point dipole is the gradient of the field vector, i.e.

B _1J = VB = -VVV,

which gives:

3m,r _J +3m _J r, + 3(ITIr)O _y ^(m.r)^

^B « = (i,j = 1,2,3).

For static magnetic fields and for the common case where conduction currents are negligible, the magnetic gradient tensor is symmetric and traceless. Then, at an observation point (x, y, 0) due to a point dipole at (x _o ,yo,h), explicit expressions for the individual gradient tensor elements are:

3m f

B _xx =-^{r ² [3L(x-x ₀ ) ₊ M(y-y ₀ )-Nh)-5(x-x ₀ ) ² [L(x-x ₀ ) + M(y-y ₀ )-Nh]}, r

3m f , B _xy =^{r[L(y-y _o ) + M(x-x _o )]-5(x-x _o )(y-y _o )[L(x-x _o ) + M(y-y _o )-Nh]}, r

B _xz =^{r ^: [-Lh + N(x-x _o )] + 5(x-x _o )h[L(x-x _o ) + M(y-y _o )-Nh]},

B _yy -^{r ² [L(x-x _o ) + 3M(y-y _o )-Nh]-5(y-y _o ) ² [L(x-x _o ) + M(y-y _o )-Nh]},

B _yz =^{r ² [-Mh + N(y-y _o )] + 5(y-y _o )h[L(x-x _o ) + M(y-y _o )-Nh]},

B _2z =^ ^m -{r ² [L(x-x _o ) + M(y-y _o )-3Nh]-5h ² [L(x-x _o ) + M(y-y _o )-Nh]}.

Note that the diagonal components are not independent because B _xx +Byy+Bzz = 0 (the tracelessness property) and that the remaining off-diagonal elements are determined by symmetry (B _yx = B _xy ; B ₂₇ = B _yz ). Step 1 of Figure 2a represents acquisition of at least all five independent gradient tensor elements, from which all nine components of the full tensor can be determined from symmetry and tracelessness, at a series of points along a traverse or at a series of times as a magnetic source passes by.

Eigenvectors of B satisfy by definition the relationship Bv = λv, for some scalar λ, which is the eigenvalue corresponding to eigenvector v. The eigenvalues are found by solving the characteristic equation det(B - λl) = 0, where I is the identity matrix. The characteristic equation is a cubic equation in λ, which therefore has three roots which we label λ|, λ ₂ , λ ₃ . Because the gradient tensor is symmetric, eigenvectors corresponding to distinct eigenvalues are mutually orthogonal and, in the case of two degenerate eigenvalues, three mutually orthogonal eigenvectors can always be found. Therefore three mutually orthogonal principal axes of the tensor, defined by unit eigenvectors, can always be found. Due to the tracelessness property and the rotational invariance of eigenvalues, the sum of the three eigenvalues is zero, so at least one eigenvalue is positive and at least one eigenvalue is negative. The intermediate eigenvalue may have either sign.

Step 2 of Figure 2a represents calculation of the eigenvalues for each of the gradient tensors acquired in step 1 , using standard matrix methods. For each measurement point, the numerical order of the eigenvalues is determined. The eigenvalue with the smallest absolute value is denoted λ ₂ for reasons that will become evident from the following.

The symmetry of the source-sensor geometry is lower for a dipole source than for a point pole. For a pole the source strength is a scalar and the geometry reflects only the source-sensor displacement vector. For a dipole, on the other hand, the source strength is a vector, the dipole moment, that has arbitrary orientation with respect to the displacement vector. The dipole moment and displacement vectors together define a plane, with normal parallel to their cross-product.

Since the gradient tensor has three orthogonal two-fold axes of symmetry, it is clear that one of these axes must coincide with the normal to the plane that contains the moment and displacement vectors. Thus one eigenvector is parallel to m x r and the other two eigenvectors lie within the plane containing m and r. It has been noted that the eigenvector μ x r has eigenvalue μ. r = μ cosφ, where φ is the angle between μ and

r. To verify this, choose a co-ordinate system (xi, x ₂ , X ₃ ) with x _\ directed along r, so that X| = r, X ₂ = X ₃ = O, and with X ₃ parallel to μ x r . Then mi = m cosφ and m ₂ = m sinφ, where φ is positive from X| towards X ₂ . The gradient tensor elements with respect to these co-ordinates are:

_, -6m.x? -6m.r -6m.f _ _λ .

B ₁ , = -J- ^L = — - = — = -2μ.r = -2μcosφ,

3m ₂ x, 3m, 3m m.

B ₁₂ = -7 f = — _ ₄ 4 " = μsmφ,

B ₁₃ =O,

R ^3m ι ^x ι 3m,

D^ = ; — = μcosφ,

B,,=0,

_3m,x, 3m.r 3m.f

D '3 ₁ 3, — ; = μ.r = μcosφ.

Thus the tensor has the form:

which has characteristic equation:

= (B _n -λXB ₂₂ -λXB ₃₃ -λ)-B, ₂ (B ₃₃ -λ)B ₁₂ =O,

. ^• .(B ₃₃ -λ)[(B _ll -λXB ₂₂ -λ)-Bf ₂ ] = O,

which gives λ = B ₃₃ = μ.f = μcosφ

or

λ ² -(B ₁₁ +B ₂₂ )λ + B,,B ₂₂ -B ² ,=0.

As expected from symmetry, an eigenvector corresponding to the eigenvalue B ₃₃ is (0,0, l) ^τ , i.e. it is directed perpendicular to the (m,r) plane:

We label this eigenvalue λ?, as it turns out that this eigenvalue is always intermediate in numerical value.

The other eigenvalues are obtained by solving the quadratic equation given above:

_ B ₁ , +B ₂₂ +V(B _n - B ₂₂ ) ² +4B ² ₂ _ -B ₃₃ +VB ² ₃ +4(B ² ₂ -B _n B ₂₂ ) λ, =

2 2

_ B ₁ , +B ₂₂ -V(B ₁₁ -B ₂₂ ) ² +4B ² ₂ _ -B ₃₃ -VB ² ₃ +4(B ² ₂ -B ₁₁ B ₂₂ ) λ ₃ = ^■ 2 2

Note that:

B ₁ , B ₁₂ 0 B ₁₂ B ₁₁ B I ₂ + B _l2 (λ-B,,)

B ₁₂ B ₂₂ 0 λ-B,, = B? ₂ + B ₂₂ (λ-B,,)

0 0 B,, 0 0

where the x ₂ component of the RHS has been simplified by rearranging the quadratic equation for λ.

Therefore eigenvectors u, corresponding to λ, (i = 1,3) are given by:

B ₁₂ ^u , = λ, -B ₁₁ 0

B I ₂ u, = λ ₃ -B _n 0

The x? components of these two eigenvectors have opposite signs, as:

_- (B ₁ ,-B ₂₂ )+ V(B ₁ , -B ₂₂ ) ² + 4Bf _{2 ^} λ, -B ₁₁ >0, 2

and

-(B _n -B ₂₂ J-V(B ₁₁ -B ₂₂ ) ² +4Bf ₂ λ ₃ -B ₁ , - <0,

with equality only if B ₁₂ = 0, and Bn > B ₂₂ (for the first case) or Bn < B22 (for the second case).

In terms of the scaled moment μ and the angle φ, the three eigenvalues in descending order are:

λ, = μcosφ,

λ ₃ = -(- cosφ - _λ /4 + 5cos ² φ]

where λ| > λ ₂ > λ ₃ and λ ₂ is always the eigenvalue with the smallest absolute value.

A commonly used rotational invariant is the contraction of the tensor on itself:

B..B = ∑∑B _u B _JI + BL .

' J

With respect to principal axes therefore

B..B = λi + λ; + λ ²

The contraction can be written in terms of the scaled moment and φ as

B..B = 2μ ^: (l + 2cos ² φ)= -^ ^m -(l + 2cos ² φ)

Thus for a given range the tensor contraction is maximum when the moment vector is aligned with the source-sensor axis (the so-called Gauss A position), and minimum when φ = 90°, i.e. when the sensor is broadside or in the Gauss B position with respect to the moment. Although the contraction measured along a track is very sensitive to the range and peaks in the vicinity of the point of closest approach, which for example

makes it a useful parameter for homing onto a target by a tensor gradiometer system mounted on a manoeuverable platform, in general the maximum does not strictly coincide with the point of closest approach, due to the dependence on φ. The surfaces of constant B..B around a point dipole are rotationally symmetric about the dipole axis. Their long axes are aligned with the dipole moment and are 1.147 times larger than the equatorial short axes.

A particularly useful rotational invariant is the scaled moment, which can be expressed as the following combination of eigenvalues:

3 3mm

- — - y- λ,- - λ,λ, ₃ .

Thus, in the absence of interfering anomalies, the scaled moment of a dipole source at any measurement point can be calculated directly from the eigenvalues of the gradient tensor at that point. The present invention recognises that a major virtue of the scaled moment is that it depends only on the range and is independent of the orientation of the dipole moment and the source-sensor axis, unlike the aforementioned invariants, which depend on φ. The surfaces of constant μ around a point dipole are spherical.

The present embodiment thus recognises that a useful rotational invariant for inversion to obtain certain source parameters is the scaled moment, which can be expressed as the following combination of eigenvalues:

3 ^mm

4~ - -\/~ λ, - λ _| λ _j r Step 3 of Figure 2a represents calculation of this initial estimate of the scaled moment from the eigenvalues. For a pure dipole signature, this estimate will correspond to the actual scaled moment of the source. If the estimated scaled moment exceeds a threshold level, indicating detection of a signal from a magnetic source above the background noise level, logic switch 4 of Figure 2a retains this value for the next stage of processing, otherwise data continues to be acquired until the scaled moment threshold is exceeded. When sufficient (> 1 1) successive values of the estimated scaled

moment have been acquired, logic switch 5 of Figure 2a initiates deconvolution of the scaled moment profile (step 6 of Figure 2a), to extract the desired characteristics of the magnetised body. At successive points x = Xj ( i = 1 , 2,..., n) along a straight and level path, defined by y - y ₀ = Y, the scaled moment determined from the measured data can be modelled as:

3m , ^μ i ⁼ 7 ; — λ T + ^a + ^bχ i + ^cχ i ^" >

((X ₁ - X ₀ ) ² + S ² ) where S = V(Y ² + h ² ) is the slant distance from the point of closest approach to the dipole, a is the base level, and b, c are linear and quadratic terms that represent interference from other anomalies and departures from a pure dipole signature. It should be noted that the method can be generalized, in a relatively straightforward fashion, to include higher order background terms in cases of large and complex interference from other anomalies, at the possible expense of some robustness of the solution. Correction of the observed gradient tensor profiles for interference terms and departures from a pure dipole signature is a unique aspect of this method and suppresses systematic errors that can afflict other solutions to the DLC problem.

The deconvolution problem is to solve for the unknown parameters X ₀ , S, m, a, b, c. Multiplying out, rearranging and writing in matrix form gives Mu = v, where:

.

3m + (Xo +S ^: ) ² a

-4x _o (x^+S ² )a + (x ² +S ² ) ² b

(6x ² +2S ² )a-4x _o (x ² +S ² )b + (x ² +S ² ) ² c

-4x ₀ a + (6x ² +2S ² )b-4x ₀ (x ² +S ² )c a-4x ₀ b + (6x ² +2S ² )c u = b-4x ₀ c c

-(x ² ₊ S ² ) ²

4x _o (x ² +S ² )

-6x ² -2S ²

4x _λ

μι ^χ ι v = μ ₂ ^χ .

M is an n x 11 matrix, u is an ll x l column vector and v is n x 1. For n > 11 this matrix equation defines an over-determined system of linear equations in combinations of the unknown parameters. The least-squares best-fit solution to this system of equations can be written formally as:

U,

U = = (M ^T M) ^~ 'M ^T V.

MO

In practice, because the square matrix M ^T M is often ill-conditioned, the matrix equation Mu = v should be solved by a numerically robust method such as QR decomposition or singular value decomposition. The unknown parameters of interest are given by:

X ₀ =u,,/4,

S = V(-u _l0 -3uf,/8)/2, or

C = U 7> b = u ₆ +u,,u ₇ , a = u ₅ +u _M (u ₆ +u _ll u ₇ ) ₇ +u ₁₀ u ₇ , or a = -[u ₂ +u ₈ (u ₆ +u,,u ₇ )]/u ₉ , m = (u, + au _g )/3.

Three independent estimates of the slant distance S and two estimates of the base level term, a, occur and can be used as a check on the solution. A more rigorous approach to estimating XQ and S is to solve the over-determined system of equations:

which has the least squares solution:

Then X ₀ and S can be determined by taking square roots, noting that S > 0 and that x ₀ has the same sign as uπ.

Once the origin and slant distance are estimated the interference parameters a, b, c can also be least-squares estimated by solving the over-determined system:

-4x _o (Xo +S-) (X ² +S ² ) ² 0

6x ² +2S ² -4x _o (x ² +S ² ) (x ² +S ² ) ²

-4x ₀ 6x ² +2S ² -(X ² +S 2 ² \)2

1 -4x ₀ 6 xx ²² ++22SS ²

0 1 -4x _λ

0 0

Finally, the moment is calculated as: m=[u,-(x ² +S ² ) ² a]/3.

An iterative procedure may be useful, where the point of closest approach x = X ₀ is identified and the calculations repeated with this value subtracted from all profile co- ordinates, i.e. with the closest approach point as origin (step 7 of Figure 2b). Then:

(M ^T M) ^" 'M ^T V«U.

The parameters of interest are directly determined, with the magnitude of the u ₂ and u ₄ terms providing an indication of error. The process can be repeated with the estimated interference removed from the original data (step 8 of Figure 2b). Logic switch 9 decides whether the fit of the model to the measured data is acceptable. If not, another data point is included in the analysis (step 10 of Figure 2b). The background-corrected scaled moment peaks at the point of closest approach, where it is equal to m/S . When the modelled scaled moment is well fitted by the data, the refined estimate of m is used in subsequent analysis (step 11 of Figure 2c).

At each measurement point the distance to the dipole is given by (step 12 of Figure 2c):

r, = V( ^x . - ^χ o) ² + S ² , where the scaled moment in this and all subsequent expressions is assumed to have been corrected for background fields.

The minor eigenvalue at each point is:

,. . , (X ₁ - X ₀ ) Y h (X ₁ - X ₀ ) S

(λ ₂ ), = μ.cosφ, = m _x ' + m _y — - m ₂ — = m _x ¹ ^ ₅ + m _z , -y ,

where z' denotes the rotated co-ordinate axis from the point of closest approach to the dipole. Therefore a linear regression of r, ⁵ (λ ₂ ), on x, can be used to estimate the along- profile magnetic moment m _λ from the slope of the linear regression and the "slant" magnetic moment m _z - from the intercept. That is,

.

Complete information about the source moment and location is unobtainable from the values of these two scalars alone. The extra information to determine these quantities must be extracted from the individual tensor elements.

Once the origin of x and slant distance are determined and the scaled moment, μ,, and distance to source, r,, at successive points are known, the five independent measured gradient tensor elements can be modelled by:

B<'>= -^2Lx? +4M'x ² -3LS ² X ₁ -M'S ³ ]+a _xx +b _xx x,+c _xx x ² ,

B^ = ^-[MX ³ -4LYX; -(4MY ² -Mh ² -5NYh)x, +LYS ² ]+a _xy +b _xy x, + c _xy x ² ,

B^ = ^-[NX ³ + 4Lhx ² +(NY ² -4Nh ² + 5MYh)x, -LhS ² ]+a _xz + b _xz x, + C _x^ x ² ,

B^ = -^-[LX ³ +(3MY-Nh)X ² -L(4Y ² -h ² )x, +MY(3h ² -2Y ² ) + Nh(4Y ² -h ² )]

^{+ a} yy ^+b yy ^x . ^+c >X' B^ =ii- [(NY-Mh)xf+ 5LYhx, +Mh(4Y ² -h ² ) + NY(Y ² -4h ² )]+a _yz +b _yz x,+c _yz x ² . where the distances along the profile, x,, are now with respect to an origin at the point of closest approach. M' is direction cosine of the slant component of magnetization and a quadratic interference term is assumed for each component. M' obeys the following relations:

M' = ^, m

M'S = MY-Nh.

The deconvolution problem is to solve for the unknown parameters L, M, N, Y, h and the interference terms a, _p b _ψ c, _y Multiplying out and rearranging gives (step 13 of Figure 2c):

Mu = v,

.-.U=(M ^T M) ^~ 'M ^T V, where M is an n x 7 matrix, u is 7 x 5 and v is an n x 5 column vector, which are given by:

x ^X3 l /μ, r ³ x /μ> ^r i ^3χ i /μ, x; r, ³ x. /μ. F ₂ ³ X* /μ ₂

M = r ₂ ³ /μ ₂

< r _n ³ /μ _n r ³ x _n /μ _n /μ _n

M 'S ³ LYS^ -LhS ² MY(3h ² -2Y ² )-Nh(h ² -4Y ² ) (NY-Mh )(S ² -5hY)

3LS ² 5NYh-4MY ² +Mh ² NY ² -4Nh ² +5MYh L(h ² -4Y ² ) 5LYh

-4 M 'S -4LY 4Lh 3MY-Nh NY-Mh u = -2L M N L 0 ^a xx ^a xy ^a xz ^a yy ^a yz ^b xx ^b xy ^b xz ^b yy ^b yz ^c xx ^c xy ^c xz Cyy ^c yz

_B (l)

/μ, r, ³ «y /μ, « /μ I B ⁽ ;; /μ, T BS /μ 1 r ³ _B (2)

B ⁽ ;>

V = /μ ₂ r ³ xy /μ ₂ /μ 2 r ³ BS ^> /μ ₂ r ₂ ³ B? /μ 2 _β (n)

/ ³ B ⁽ ;: /μ» r ³ xy /μ _n T ³ B ¹ ^ /μ n r ³ B£ ^> /μ _n r ³ B ⁽ ; _Z > /μ n

Step 14 of Figure 2c represents solving this overdetermined matrix equation by standard methods. The interference parameters are obtained directly from the 5 generalized inverse matrix u: ^a xs = ^U 5l ' ^a xy ^{= U} 52> ^a xz ^{= U} 53> ^a yy = ^U 54 > ^ = U ₅₅ , b _x , = u _6l , b _xy = u ₆₂ , b _x2 = U ₆₃ , b _yy = u ₆₄ , b _yz = u ₆₅ ,

^C xx = ^u _7l >C _xy ^U^.C^ =U ₇₃ ,C _yy =U ₇₄ ,C _yz =U ₇₅ .

At this stage it may be optimal to remove the interference terms from the measured tensor elements and recalculate the eigenvalues and scaled moment (step 15 of Figure 10 2c). Using the new estimates of X ₀ , S, m the deconvolution of the tensor elements can be repeated. If the process is convergent, as expected, the revised interference terms should be small (as tested by logic switch 16 of Figure 2d) and the source parameters should be more precisely determined.

The other parameters can be obtained from the elements of u in a sequential manner (step 17 of Figure 2d). First we have: M = u.,

N = u

Because the other parameters occur in more than one element, each subsequent step can be a least-squares estimate. The least-squares estimate for L is: u,, - 2u<

L = * 44 Ml

Now that the direction cosines of the magnetic moment (L, M, N) are determined, the parameters Y and h can be estimated by inverting the overdetermined system:

At this stage all the parameters of the dipole model have been estimated. However, a further check on the consistency of the solutions can be made by inverting the over- determined system:

The values of Y ² , h ² and Yh determined from the generalized inverse of this matrix equation should be consistent with the values of Y and h estimated from the previous step, and the sum Y ² + h ² should be close to the value of S ² obtained from the analysis of the scaled moment. When the solution of the source and background parameters is complete for the given data window, the data window can be shifted and the entire process repeated (step 18 of Figure 2d).

While the embodiment described in the preceding models the scaled moment, other parameters determined from the measured data may be suitable for modeling in alternative embodiments of the invention. One such alternative embodiment is described in the following, which provides a simplified deconvolution of the scaled moment profile.

Step 19 of Figure 3a represents acquisition of all five independent gradient tensor elements, from which all nine components of the full tensor can be determined from symmetry and tracelessness, at a series of points along a traverse or at a series of times as a magnetic source passes by. Step 20 of Figure 3a represents calculation of the eigenvalues for each of the gradient tensors acquired in step 1, using standard matrix methods. For each measurement point, the numerical order of the eigenvalues is determined. In this embodiment, a new parameter v is defined at each point by: μ Vm

V = . ' -

3 (x - x _o ) ² + S ^{2 '}

At successive points x = Xj ( i = 1, 2,..., n) along a straight and level path, defined by y - yo = Y, this new parameter, as determined from the measured data (step 21 of Figure 3a), can be modelled by:

where S = V(Y ² + h ² ) is the slant distance from the point of closest approach to the dipole, a is the base level, and b, c are linear and quadratic terms that represent interference from other anomalies.

If the estimated value of v exceeds a threshold level, indicating detection of a signal from a magnetic source above the background noise level, logic switch 22 of Figure 3a retains this value for the next stage of processing, otherwise data continue to be acquired until the scaled moment threshold is exceeded. When sufficient (> 7) successive values of the estimated v have been acquired, logic switch 23 of Figure 3a initiates deconvolution of the scaled moment profile (step 24 of Figure 3a), to extract the desired characteristics of the magnetised body.

The deconvolution problem in this embodiment is to solve for the unknown parameters X ₀ , S, m, a, b, c. Advantageously, this is equivalent to conventional Werner deconvolution of the TMI anomaly of a thin sheet. Multiplying out, rearranging and writing in matrix form gives Mu = v, where:

1 ^x , x. x; x; v, V ₁ X ₁

1 X ₂ V ₂ X ,

M = x. X ³ : x.

_n

V ₁ X ^"

V =

V _n X ²

Vm +a(xjj +S ² ) b(x ² +S ² )-2x ₀ a a-2x ₀ b + c(X _g + S ² )

U = b-2x ₀ c

2x

M is an n x 7 matrix, u is a 7 x 1 column vector and v is n x 1. For n > 7 this matrix equation defines an over-determined system of linear equations in combinations of the unknown parameters. The least-squares best-fit solution to this system of equations can be written formally as the pseudo inverse: u I u 2 u 3 u = u 4 = (M ^T M) ^" 'M ^T v.

U 5

U 6 u 7

In practice, because the square matrix M ^T M is often ill-conditioned, the matrix equation Mu = v should be solved by a numerically robust method such as QR decomposition or singular value decomposition.

The unknown parameters of interest are given by:

X ₀ = u ₇ /2,

C = U 5' b = u ₄ +u ₅ u ₇ , a = u ₃ + u ₇ b + u ₅ u ₆ , or a = -(u, + u ₆ b)/u ₇ , m = (u,+au ₆ ) ² .

As can be seen two estimates of the base level term, a, occur and can be used as a check on the solution. A more rigorous approach to estimating a, b, c is to solve the over- determined system of equations:

0

1 - ^U 7 - u

0 1 - u

0 0 1

An iterative procedure may be useful, where the point of closest approach x = x ₀ is identified and the calculations repeated with this value subtracted from all profile coordinates, i.e. with the closest approach point as origin (step 25 of Figure 3b). Then:

(M ^T M) ^~ ' M ^T V * U . The parameters of interest are directly determined, with the magnitude of the U ₂ and U ₄ terms providing an indication of error. The process can be repeated with the estimated interference removed from the original data (step 26 of Figure 3b). Logic switch 27 of Figure 3b decides whether the fit of the model to the measured data is acceptable. If not, another data point is included in the analysis (step 28 of Figure 3b). The background-corrected scaled moment peaks at the point of closest approach, where it is equal to Vm/S ² . When the modelled v is well fitted by the data, the refined estimate of m is used in subsequent analysis (step 29 of Figure 3c).

At each measurement point the distance to the dipole is given by (step 30 of Figure 3c)

A number of practical considerations apply in implementing the present invention.

Solutions may be numerically unstable if the algorithm is applied to a short, quasilinear segment of data, particularly on the flanks of the profile caused by the anomaly. It is therefore desirable to fit a substantial portion of the anomaly, either by including many

data points in the fit or, if this is too demanding computationally, fitting every n ^th point over a segment, then moving the data window by one point and repeating the inversion. Successive solutions that satisfy reliability criteria can be then averaged, after they have been reduced to a common origin.

Solutions are more robust if the interference term is simplified. For example if the interference is assumed to be a constant base level shift over the fitted segment, the matrix equations simplify to:

1 x , V ₁ X , x: V ₂ X ,

M =

1 x. ^V n ^X n

V ₁ X ^" v, x ² , v =

V _n X

M is an n x 5 matrix, u is an 5 x 1 column vector and v is n x 1. For n > 5 this matrix equation defines an over-determined system of linear equations in combinations of the unknown parameters. The least-squares best-fit solution to this system of equations can be written formally as: u, u, u = = (M ^T M) ^~ ' M ^T V.

The unknown parameters of interest are given by: x ₀ = u ₅ /2,

S = -^(-U ₄ - Xo ), a = (u, - u ₂ /u ₅ )/2, m = (u, + au ₄ ) ² .

As the data window is shifted, the background level will change slowly. A linear or quadratic fit to the estimated values of a as a function of x can then be used to subtract the background from the whole data set and a more robust fit to the source parameters obtained from the background-corrected data.

An alternative approach to determining the unknown parameters of the target is to undertake a non-linear least-squares inversion, using well known techniques. In this approach we model measured values of the parameter v, sampled at x = X|, X ₂ ,..., x _n along a profile by:

In this case the optimisation problem is to find estimates of the parameters x ₀ , S, m, a, b, c that minimise the sum of shares of the residuals δv,, which are given by:

^δv ' ^{= v} ' ^~ l τ ^m , + a + bx, + cx, ² , (i = l, 2 n).

The Gauss-Newton algorithm employs a linearised estimate of the changes δp = (xo, S, m, a, b, c) ^τ that are needed to reduce the residuals to zero. In matrix form, this yields:

δv = = Jδp,

where

A solution for the six parameters requires n > 6. With n > 6 the system of equations is over-determined, with a least-squares solution obtained by solving the normal equations:

J ^τ Jδp = J ^τ δv.

The updated set of parameters p + δp is then used to calculate a new set of residuals and the process is repeated until the fit is acceptable.

Reasonable initial guesses for the parameters are the approximate location of the maximum value of v for xo, the approximate half width at half maximum (HWHM) for S, the square of the HWHM value of S, multiplied by v _max for Vm, and zero for a, b, c, which is equivalent to ignoring the background. Alternatively, the average value of v well away from the anomaly peak could be used as the initial estimate of a.

Dividing v by v _max and normalising all values of horizontal co-ordinate to S makes all quantities dimensionless and aids numerical stability.

The Levenberg-Marquardt algorithm modifies the inversion by adding a term to the LHS of the normal equations:

(J ^τ J + λI)δp = J ^τ δv, where λ > 0 is an adjustable damping parameter that controls the rate of convergence and the stability and I is the identity matrix. This algorithm is a compromise between the Gauss-Newton method, which converges rapidly but may require good initial estimates, and the method of steepest descent, which is robust but may converge slowly. Assuming that the Jacobian matrix has been normalised appropriately, so that diagonal elements are of the order unity, a suitable initial choice of the damping parameter might be λ = 10-20, which corresponds to an initial steepest descent approach. After each reduction in the squared residual, λ can be reduced by a factor of 10, until no further improvements in the fit occur. For λ « 1, the method reduces essentially to the Gauss-Newton algorithm, ensuring rapid convergence.

At each measurement point the distance to the dipole is given by:

where the parameter v in this and all subsequent expressions is assumed to have been corrected for background fields.

Once the origin of x and slant distance are determined and the background-corrected scaled moment, μι,= 3v; ² , and distance to source, η, at successive points are known, the five independent measured gradient tensor elements can be analysed as described above to estimate the remaining parameters of the dipole: Y, h, L, M, N and the interference terms for each of the gradient tensor elements. In this second embodiment, these parameters are calculated using steps 29-33 of Figure 3c and steps 34-36 of Figure 3d, in the same fashion to steps 1 1-18 of Figure 2.

The present invention thus recognises that gradient measurements are relatively insensitive to orientation. This is because gradients arise largely from anomalous sources, and the background gradient is low, in contrast to the field vector, which is dominated by the background field arising from the Earth's core. Gradient measurements are therefore most appropriate for mobile sensor applications. Another advantage of gradient measurements is that they obviate the need for base stations through common-mode rejection of time-varying geomagnetic fields during the survey. They also greatly reduce the need for regional corrections, which are required by TMI surveys due to fields arising from sources that are deeper or broader than the features of interest, or to the normal (quasi-) latitudinal intensity variation of the global field. However, it is notable that variable background gradients of geological origin or from neighbouring sources do present problems for conventional methods of dipole location/tracking, and the present invention is designed to minimise these problems by estimating and removing these effects. Gradient measurements also provide significantly higher resolution of closely spaced sources than TMI or vector field measurements.

Gradient measurements also provide valuable additional information, compared to conventional total field measurements, when the field is undersampled. Undersampling is common perpendicular to flight lines in airborne surveys, is usual in ground surveys, and always pertains in down-hole surveys. Although TMI gradients can be useful in the above-mentioned applications, the present invention recognises that full magnetic gradient tensor contains significantly more information. In particular, measurements of

the full magnetic gradient tensor can be used to uniquely locate and characterise a compact magnetic source from measurements at only a few nearby locations. Furthermore, rotational invariant quantities, which are independent of sensor orientation, have simple mathematical form and are relatively easy to interpret, and can be calculated from the gradient tensor.

Embodiments of the invention may thus provide improved imaging and mapping of compact source locations and their magnetic moments, even in the presence of instrumental drifts, relatively high levels of noise or substantial interference from deeper and/or broader anomaly sources, while suppressing effects of time-varying magnetic fields, regional geological gradients and errors in orientation of the sensor. In mapping applications, embodiments enabling identification of source location from a single nearby pass further may allow flight line spacing to be doubled, for example, and/or may permit movement of the magnetised body during a time period from one set of measurements to the next to be detected and determined. Moreover, the ability to determine location from a single pass means that flight lines do not need to pass directly over, or on both sides of, the magnetised body, but may suffice if passing only to one side of the magnetised body.

Embodiments of the present invention may thus be applicable to data acquired by a single tensor gradiometer of any type, acquiring data along an approximately linear profile segment. Embodiments of the invention use specific rotational invariants, calculated directly and simply from the eigenvalues of the tensor, which have the virtues of perfect isotropy around the dipole source, eliminating asphericity errors. These invariants are the scaled moment and a quantity proportional to the square root of the scaled moment. Preferred embodiments do not require calculation of eigenvectors or calculation of the tensor with respect to principal axes, thereby reducing computational requirements and avoiding numerical instability of solutions.

Preferred embodiments of the invention are relatively simple mathematically and rely on readily implemented numerical algorithms for extraction of dipole source

parameters, and further improve signal-to-noise ratio by incorporating sufficient measurements from a time series to average out effects of individually noisy measurements. Preferred embodiments simplify inversion of the segment of gradient tensor data, to obtain the dipole location and moment vector, by separating the scalar and vector aspects of the inversion. Preferred embodiments further provide for estimating and removing background interference terms, drifts and non-dipole terms, which enables the final solution to be based on an almost pure dipole signature, minimising or avoiding systematic errors.

It will be appreciated by persons skilled in the art that numerous variations and/or modifications may be made to the invention as shown in the specific embodiments without departing from the scope of the invention as broadly described. For example, the case of a stationary sensor and a moving magnetic body may arise in embodiments of the present invention, such as systems for detecting metal in items such as food or laundry on a conveyor belt, security systems at doors or portals, or the like. Where the magnetic body moves past the sensor at an unknown substantially constant speed, a single sensor may yield a solution which is not unique for: (a) range at the point of closest approach; and (b) magnetic moment magnitude. In such cases a unique solution may be obtained by use of additional sensor(s) at a suitable known location nearby, and/or by incorporating additional information regarding the point of closest approach or other constraints. The present embodiments are, therefore, to be considered in all respects as illustrative and not restrictive.