This invention relates to the analysis of subterranean objects and, more particularly, to the determination of properties of subterranean objects by using principles of computerised tomography to interpret magnetic anomaly data associated with those objects. Background of the Invention

Magnetic surveying is a well-established method of conducting geophysical exploration based on detecting anomalies in the Earth's magnetic field. A limitation of this method at present is that it is difficult to use with three dimensional modelling in order to obtain vertical geological sections or profiles from measured magnetic anomaly field data. In the past, efforts have been made to derive automatic methods for interpretation of the magnetic or potential anomaly profiles, but the most widely used methods at present are still various trial- and-error methods which require good a priori knowledge concerning the object being analyzed.

Suppose that there is a buried two-dimensional magnetic source body of area S with constant magnetic dipole moment per unit area M whose direction is determined by its direction cosines =cosα and Q = cosβ where α is M's angle from the (horizontal) x-axis and β is M ' s angle from the (vertical) z-axis. The Earth's magnetic field is measured along a line parallel to the x-axis at elevation z ₀ above the Earth's surface. Suppose also that the direction of the total magnetic field at measuring point (x',∑ ₀ ) is a unit vector t whose direction cosines are p = cos ' and q = cosβ' . Then the total magnetic field measured by a magnetometer at point (x',∑ ₀ ) can be expressed as [l,2]

b _a (x',z ₀ ) = -t - VG(x',∑ ₀ ,x,z) (1)

The subscript a denotes that the field intensity mentioned here is the net anomaly field of the source body and (x,z) are the coordinates of a point within the source body. The gradient operator V is defined as

V=— i +— j dx dz

where / ^' and j are unit vectors in the directions of the x and z axes, and G(x',z ₀ ,x,z) is the so-called magnetic scalar potential which is defined as

where

Since

M-V = M\ P— + 0— )

¹ ά ~ dz)

where M is the magnitude of vector M , we have

Similarly, because

δ δ ϊ - = p — + q — δx δz

equation (1) becomes

where K is the data kernel defined as

Here, P,Q,p and q are considered as constants. Modern magnetometers, such as three-component magnetometers and differential vector magnetometers, can provide information about these parameters. Thus, the data kernel can be considered to be a known function.

The interpretation task is to solve equation (3) to find M and to determine the boundary of S from the measured field data. From the point of view of image processing, this interpretation task is to reconstruct, from the measured signal, a magnetic dipole density distribution function which describes the source body.

If there is a priori knowledge about the shape of the boundary, the location, and the value of M for the source body, then the interpretation can be accomplished in a straightforward manner using least squares curve fitting. This is the principle of various trial-and-error strategies where the a priori knowledge may either come from well-log data or the interpreter's guess work. Without good a priori knowledge, an hypothesis that is often adopted is that the boundary of S is contained

in an area of interest S' (/.e., S cz S') with a new distribution function defined as [3].

Therefore the governing formula for the interpretation of the magnetic anomaly profile becomes

(2.7) b _a (x',z ₀ ) = Jl _sι K(x',z ₀ ,x,z)w(x,z)ώcdz (6)

Equation (6) is a Fredholm integral equation of the first kind, and is an ill-posed problem in the sense of Hadamard's standard [4]. If the distribution function w(x,z) could be reconstructed by solving the integral equation then the interpretation of magnetic anomaly profile is accomplished or, in other words, the image of the source body is reconstructed from the measured signal. Since the real data are measured at discrete points

(x',z ₀ ) rather than continuously for all values of x', one method of solving equation (6) is to convert it into a system of linear equations of the form

Ay = b (7)

where A is an m x n matrix and y and b are n x \ and m x \ vectors, respectively. There are many numerical quadrature schemes available for approximating an integral equation as a system of linear equations [5,6]. However, mid-point quadrature has been found [9] "to be a good scheme for ill-posed Fredholm integrals of the first kind. Examples which adopt such a quadrature scheme in geophysical inversion problems are discussed in references

[7,8]. One such method is Green's strategy [7]. To use Green's strategy for converting equation (6) into a linear system of equations, suppose that the area of interest S' is a rectangle given by the Cartesian product

S' = (b -a) x(d-c) (8)

where (b-a) is the length of S' and (d-c) is the width. Divide S' into JK pixels, each of length Δx and width Δz, where

Δx = (b-a)/ J, md Δz = (d-c)/K.

Then the coordinates of the centre of pixel (j,k) are

x _J = -l/2)Δx; = 1,2,..., J and

z _k =(*-l/2)Δz; k = \,2,...,K

For simplicity, assume that the direction cosines of M and t are P = p = 0 and Q = q = \ . The mid-point quadrature then defines the discrete kernel function as

With a finite number of measuring points, equation (6) becomes

where N is the number of measuring points, A, _jk is given by

and w _jk is a constant within a pixel but may change from pixel to pixel. In terms of matrices equation (9) can be rewritten as

Aw = b _a (10)

where A is a (Nx JK) matrix and w and b _a are JK and N vectors, respectively. Now the governing formula (6) has been converted into the system of linear equations (10). This linear system, from the point of view of inversion theory, serves as the forward model in discrete form.

Unfortunately, it can be shown that matrix A is ill-conditioned. Therefore, it is difficult to obtain an accurate estimate of w even if very accurate measurements b _a are available.

Further, if the number of pixels is increased in an attempt to obtain better accuracy, A becomes even more ill-conditioned. Yet further, A is a large dense matrix, so computing w from (10) is computationally intense, and thus the strategy for solving equation (6) by dividing the area of interest into pixels is not optimal. Therefore, in the following, we describe a method based on computerised tomography that reduces some of these problems of finding the boundary of S .

Summary of the Invention

The present invention consists in a method for obtaining magnetic dipole density image data concerning a subterranean object S contained in an area of interest S' , the object having properties that cause a total magnetic field comprising a vector sum of the Earth's magnetic field and an anomaly, the method comprising the steps of: measuring the total magnetic field at a plurality of locations in a vicinity of the object to obtain a plurality of magnetic field observations; dividing the area of interest S' into a set of narrow parallel strips, all of which are perpendicular to a direction which is specified by an angle θ ; relative to a primary axis of a reference coordinate system; computing a strip integral of the plurality of magnetic field observations for each one of the plurality of locations; forming and solving a system of linear equations using the strip integrals to obtain a generalised projection; repeating the above three steps for different angles θ, where the values of θ are evenly distributed within the range of [0, π] performing an inverse Radon transformation on the generalised projections to obtain the magnetic dipole moment density image data concerning the object.

A representation of the magnetic dipole moment density image data may then be displayed in any suitable format. Preferably, the strip integral _Q is computed according to

^a * ⁼ \ _k ' ^κ(χ * ' ^z ° ^,x ' ^ ^,dxdz

where / denotes a location index, k denotes a strip index, S _k ' is the kth strip of an area of interest, K(x _κ z ₀ ,x,z) is a known kernel function, x _i is the / th location, z ₀ is the height of the plurality of locations, and w' is a constant substantially equal to a magnetic dipole moment per unit area of the object.

Preferably, the step of forming and solving the set of linear equations further comprises the steps of: forming a first intermediate system of linear equations of the form Ya = b where Y is an unknown matrix, a is a vector obtained form the strip integrals, and b is a vector formed from the plurality of magnetic field observations; and converting the first intermediate system of linear equations into a second intermediate system of linear equations of the form A *y _i = b _i where A is a matrix of coefficients obtained from the first intermediate system of linear equations, y _t is the / th row of the unknown matrix Y , and b _t is a vector formed from the total magnetic field.

Preferably, the step of solving the second intermediate system of linear equations uses a regularization method. Brief Description of the Drawings The invention will now be described with reference to the accompanying drawings in which:-

Figure la is a schematic cross-sectional representation of a subterranean object S contained in an area of interest S' ; Figure lb is a graph representing the measured magnetic field in the vicinity of the subterranean object, the measured field being caused by the object's distortion of the Earth's magnetic field;

Figure 2 illustrates how an area of interest and a subterranean object is partitioned into strips perpendicular to a direction specified by an angle θ; Figure 3 illustrates a vector summation of the Earth's magnetic field and an anomaly field to produce a total field;

Figure 4 is a schematic cross-sectional representation of a subterranean object in an area of interest, illustrating how the area of interest in an example is divided into strips;

Figure 5 is a graph of simulated magnetic field density for the example of Figure 4;

Figure 6 is a generalised projection computed in accordance with the present invention; and Figure 7 is the reconstructed image of Figure 4 using projections in 33 directions in accordance with the present invention.

Detailed Description of the Invention

As described above, the task at hand is to reconstruct a magnetic dipole density image of the cross section of S' from which the boundary of S can be determined. The following describes a method for accomplishing this task using techniques related to computerised tomography. Computerised tomography (CT) and tomographic image reconstruction techniques provide a well-defined automatic method for three-dimensional modelling using anomaly data. Such techniques have been extensively studied and developed for applications in medicine and in some forms of non-destructive testing. They have also been used for seismic data processing in the petroleum industry, and for studying hidden or buried objects.

Briefly, computerised tomography relates to the reconstruction of an unknown function from line integrals. In particular, it relies upon the property that an unknown function / in a region S' can be uniquely reconstructed

from line integrals of / taken along all directions through S'.

Generally speaking, modern CT methods are applicable to imaging situations where line integrals or strip integrals of a parameter, such as X-ray attenuation in medical imaging or wave slowness in seismic tomography, are available as collected data.

To define the operation of computerised tomography more precisely, consider a plane in two dimensions in which any point can be represented in terms of its

Cartesian coordinates JC and y , or its polar coordinates τ and φ . For an arbitrary straight line in the plane, let φ denote the angle between the x-axis and the direction normal to the line, and let τ denote the signed perpendicular distance from the origin to the line. The line can then be described by the equation

xcos#+.ysinr9=r (11)

Let (r,6? ) be the Radon transform [11], line integral, or projection of f(x,y) . Then / in polar coordinates can be reconstructed uniquely from F(τ,θ) by the inverse transformation

_f ⁽ _η ,φ) = - π ^∞ ' dθdτ (12)

where F{τ,θ) is the partial derivative of F{τ,θ) with respect to τ .

In practice, F(τ,θ) is only known at discrete points along lines, i.e., for τ = τ _k = kη for nonnegative integer values of k , and Θ = Θ _J = jπ/ N;j =0,1,2...,N- 1. Here, each value of θ defines one of N views taken around a semicircle

surrounding the object, and η is the separation between two adjacent parallel lines in each view.

Hence, in practical implementations, the integrals in equation (12) are replaced by summations to obtain an approximate reconstruction formula for f _ξ , depending on ξ , given by

f _ξ (x,y) = ₊ ycosθ _J - τ _k ) . (13)

Here, Aθ is the angle between two adjacent views, ξ is a filter function introduced by regularisation related to approximation theory and L is the number of line integrals within a view. It can be shown [13] that L and

N are related optimally by the expression L = N/ . Image reconstructing techniques based on equation (13) are known to as back projection methods.

Now apply the principles of computerised tomography to the anomaly interpretation task at hand. Fig. la illustrates a schematic cross-section of object 10 beneath Earth's surface 12. Object 10 is contained within area of interest S' . As shown in Fig. la, an x-axis is defined to be oriented parallel to the Earth's surface, and a z-axis is defined to be oriented vertically. Object 10 has a magnetic dipole moment density M whose direction is defined by its direction cosines =cosαr and Q = cosβ where a is M ' s angle from the x-axis and β is M ' s angle from the z-axis.

The Earth's magnetic field is measured along a line parallel to the x-axis at elevation z _o above the Earth's surface giving an anomaly as illustrated by curve 14 in Fig. lb.

To apply the principles of computerised tomography, instead of dividing the area of interest S' directly into a number of pixels as is conventionally done, divide it into a number of narrow strips as shown in Fig. 2. For a given angle θ with respect to the x-axis, draw parallel lines 16 with the same interval η between two adjacent lines, the normal to all lines making an angle θ with the x-axis such that they can all be described by equation (11). Let the area between two adjacent lines define a strip. Then η is the width of the strips, and is chosen (to ensure that a desired spatial resolution is obtained in the reconstructed image [10]) such that 77<min(<St,&), where δxδz is the size of a pixel to be produced in the reconstructed. Let S _k ' denote the area of the A:th strip, and let A _k denote the intersection of the source body S with S _k '

'* =JJ _Δέ K(x„z ₀ ,x,z)ω(x,z)dxdz ; i=1,2,..., N (14)

where t _ik represents the contribution of the source body from the k h strip to the measured field intensity at the / ^' th measuring point, and N is the total number of measuring points. It follows that

b _a (ή = ∑ty, = ,..., N (15) t=ι

where N is the total number of the strips. (For the sake of simplicity we set the number of strips equal to the number of measuring points) .

Let a _Λ denotes the strip integral

= _s K(x„z ₀ ,x,z)w'dxd∑ ι=U N (16)

where w' is a constant and w'« M

Since a _Λ and t _Λ are numbers, it will always be true that

where y _Λ is an unknown real number and a _k = a _a .

Substituting equation (17) into equation (15) results in

N i}) = y _Λ ay, i = ,2,...,N (18) k=\

This can be written in the matrix form

Ya = b ( 19 )

where

and, for consistency, the field intensity b _a (J) is denoted by b, .

We now have a system of linear equations (19) in which the unknown is the matrix Y . Unlike equation (10),

this system theoretically has no approximation error (ignoring round off errors) because it is not a discrete form of a Fredholm integral of the first kind.

If the / ^' th row of Y is computed from equation (19) then, according to (18), the contribution of the source body to the / ^' th measuring point from the k strip is simply t _Λ = y _Λ a _k . It is then possible to determine the contribution of the source body from all of these strips. Defining an image function as

f(x,∑) = K(x„z ₀ ,x,z)w(x,z); (X,Z) €S' (21)

it is evident that

y, _k a _k = jj _s (x,z)dxdz . (22)

That is, that y _tk a _k is the strip integral of the image function f(x,z), and we refer to it as the generalised Radon transform. Since all of these strips are divided in parallel in a direction specified by the angle θ, the set |.y _Λ α _t |&= l,2,...,N> can be considered to be a projection of an image function /(x,z) according to the definition used in computerised tomography as described above.

For each value of θ, there is a set of strips and a corresponding projection of the image function as defined in equation (21). Construct and solve a system of linear equations for each value of θ given by θ _j =jπ/N, j = 0,\,...,N- \ to produce N projections of the image function /(x,z) along N directions. With these N projections the image /(x,z) is readily reconstructed through the back projection method of equation (13). Once the image /(x,z) is reconstructed, it is possible to find

the distribution function w(x,z) , since the data kernel

K{x _j ,∑ ₀ ,x,z) is known.

Recall from its definition in equation (4) that the data kernel is considered to be a function of the argument (x',z ₀ ,x,z), while P,Q,p and q are all considered as constants. For a given source body, P and Q are , indeed, constants. However, p and q , the directional cosines of the unit vector t , are variable; when t is parallel to the total magnetic field intensity, the measured magnetic field data is the magnitude of the total magnetic field intensity, otherwise the data is the component of the total field in the direction of t . Since p and q satisfy

pyq' = \ , (23)

the data kernel K can be considered to be a function of the extended argument (x',z ₀ ,x,z,p) . That is

For a given direction t , , (i.e., given p _i ) define vector α' as

α ^J = [α{, ..., α _k ^J , ..., α„ ^J ] (25)

where the elements αi are computed using

al = JK(x _k ,z ₀ ,x,z,p _} )w'dxdz for * = 1,2,...,N . ( 26 )

For a given p _} , denote the measured anomaly field data (a component of the total field in the direction of tj ) by the vector b ¹ which is defined as

Then one has

Ya ^J = b'; y=l,2,...,«. (28)

Defining the (nxn) matrices A and B by juxtaposing vectors as follows

A = (a ^l a ² ... a"), B = (b ^λ b ² ... b"), (29)

enables equation (28) to be written as

YA = B . (30)

Equation (30) has a unique solution if A is a non- singular or, equivalently, if the vector set (α-'| =l,2,...,«jis a linearly independent set. From the definition of a ¹ , it is clear that the linearity of this vector set is determined by the linearity of the data kernel K\x _k ,∑ ₀ ,x,z,p with respect to the variable P _j . However, K\ x _k ,∑ ₀ ,x,∑,p is a nonlinear function of p

because, from (23), q = jl -p ² . Therefore, a a ² ,...,a" are

17

linearly independent vectors, A is non-singular, and (30)

Vhiaase aa πun ni i rqriuioe esnoll nuit-i ϊ ΛonΠ rgri-i wvoen

Y = BA- (31)

where A ^'1 denotes the inverse of A.

It was indicated before that only one row of the solution matrix is required. If one transposes equation (30) and multiplies the transposed matrix equation by the vector e, where e, is an n-dimensional vector with a one in its ith position and zeros elsewhere, then one has

AΥe. = B'e. (33)

which gives

y, = (ATb, (34)

where y _t is the /th row of the unknown matrix Y and b _t is the / th row of the matrix B .

If we rewrite equation (30) explicitly as

(35)

then we can see that the elements of b _t , the th row of the matrix B , are the components in the direction p _} of the total field intensity measured at the same point. That is

b, = [b:,b ² ,...,b,y...,b:] ( 36 )

The problem that arises here is how to obtain the vector sets {α ^J |/=l,2,..., and =1,2,...,«j in practice. The magnetic field intensity measured by most modern three- component magnetometers is that produced by vector summation of the Earth's main field and the anomaly field introduced by the magnetic body. This vector summation is shown schematically in Fig. 3 where H represents the Earth's main field, C represents the anomaly field, T represents the total field measured by the magnetometer, and γ is the angle between H and T . The Earth's main field at a surveying point can reasonably be considered to be a known constant. Thus, the anomaly field C can be obtained simply from

C = T~H (37)

Once C is computed, one can find

p,=cos ' (38)

where a' is the angle between C and the x-axis. Let

b = b, (39)

where b _t is the magnitude of C . We have

_Pj = a + (j- l)Aa (40)

where Δα is the difference between angles p _} and p _J+] and (/-l)Δα is the biasing angle of the /?. from /?, . It follows that

Z>, = £ _. ¹ cos[(/ ^' - l)Δα] . ( 41 )

Therefore, b _t can be obtained from the magnitude of the total field intensity and direction cosines cos[(/-l)Δαl at the measuring point (x _l z ₀ ) • By substituting _y |/=1,2,...,«J into equation (24), the vector set α|/=1,2...,«J defined in equation (25) is readily obtained by equation (26).

The problem has now been reduced to that of solving an ordinary system of linear equations (34). However, the vector set |σ-'|y ^' =l,2...,nj may be far from orthogonal, even though it is linearly independent. Therefore, A may be ill-conditioned or rank deficient. Thus, standard methods of linear algebra for solving (34), such as LU, Cholesky, or QR factorization, may not be reliable if one uses them in a straightforward manner to compute the solution to the problem.

The most common methods used to solve systems of linear equations involving ill-conditioned matrices are regularization methods. The simplest of these is the so- called damped least squares method in the Marquardt-Levenberg's formulation

y = (A'A + M)-' A ^* b (42)

where A is a damping factor or regularization parameter, and / is the identity matrix. An important issue in the use of this method is the choice of an optimal λ, and the generalised cross validation (GCV) method is commonly used for this purpose [12]. Quite often, the GCV method fails to find an optimal λ , and there now exists a more powerful L-curve method for choosing an optimal λ .

The following describes a numerical example that illustrates the use of the method for interpreting magnetic anomaly profiles. The object considered here is a square object S within area of interest S' as shown in Fig. 4. The labels 1,2, ...,28,...,66 surrounding area of interest S' refer to strip numbers within S' . The dimensions of the pixels were set at δxδz=0.5/wx 05m . The total number, N , of pixels within area of interest S' is 61x21=1281. This requires S' to be divided into 66 strips for a given angle θ and the total number of views is set to be 33. For the chosen object, with strips in a direction 0«54.5°, only strips 28 to 40 intersect with the object, as shown in Fig. 4.

Measurements were simulated at an elevation z _o =0_5m above the surface of the Earth using equation (3). This simulated anomaly field data is shown as line 16 in Fig. 5, where the simulated magnetic field measurements are plotted against the strip number. The subscript / in equation (39) was chosen to be ι=32. Random errors δ _t were added to this perfect data to give b = b ₃₂ + δ _i2 .

A corresponding projection of the object-related image function as defined in (21) is shown as line 18 in Fig. 6, with a typical λ chosen using the GCV method. In Fig. 6, the projection is a plot of strip integrals of the image function against the strip number, and it can be seen that the strip integrals are significantly different from zero at least for strips 30 to 37.

Figure 7 shows the reconstructed image function

/(x,z) defined in equation (21). This is accomplished by using a back projection algorithm with projections obtained by solving 33 systems of linear equations.

Although in the conventional application of CT methods, the image of an object is obtained from projections directly measured at different angles around

the object, the present invention illustrates how the principles of CT can be used to obtain the image of an object using only the potential field measured on a plane located some distance away from the object. Some advantages of the method presented here are that: (1) it can find a physically acceptable solution to a Fredholm integral equation of the first kind without any a priori knowledge about the solution, (2) its computational cost is low, and (3) it can be extended easily to three-dimensional problems.

Thus, there has been described a method for interpreting magnetic anomaly profiles. Those skilled in the art will appreciate that various modifications can be made without departing from the scope of the invention as broadly described, and the embodiment described is presented for purposes of illustration and not limitation.