AU681492B2 - Method of interpreting potential field data - Google Patents

Description

WO 95/20174 PCT/AU95/00039 Method of Interpreting Potential Field Data Field of the Invention 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 trialand-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 P=cosa and Q=cosp where a is M's angle from the (horizontal) x-axis and 3 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 zo above the Earth's surface.

Suppose also that the direction of the total magnetic field at measuring point (x',Z 0 is a unit vector T whose direction cosines are p=cosa' and q=cos6'. Then the total magnetic field measured by a magnetometer at point (x',Z 0 can be expressed as [1,2] SUBSTITUTE SHEET (Rule 26) WO 95/20174 PCT/A'U95/00039 2 =TVG z 0 The subscript a denotes that the f ield intensity mentioned here is the net anomaly field of the source body and are the coordinates of a point within the source body. The gradient operator V is def ined as a- V=i+-J where i and j are unit vectors in the directions of the x and z axes, and G(x',ZOxz) is the so-called miagnetic scalar potential which is defined as G x, z) -f V~dxdz where r= 2 2 Since M.?V=Af(P 0 where M is the magnitude of vec IrfM, we have x, M [P(x +Q(z z,,)JIdz (2) Similarly, because WO 95/20174 PCT/AU95/00039 3 t V=p-+qequation becomes b x,z)Mdxdz (3) where K is the data kernel defined as K(x',x 3[P(x Q(z z)] x (Z zo) 2 5 2 P+qQ (4) (z 2 ]3/ 2 Here, P,Qpandq are considered as constants. Modern magnetometers, such as three-component magnetometers and differential vector magnetometers, can provide information about these paraneters. Thus, the data kernel can be considered to be a known function.

The interpretation task is to solve equation 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 WO 95/20174 PCT/AU95/0039 4 in an area of interest with a new distribution function defined as ,\Mif(x,z)eS w(x,z) M if(xz) )0 otherwise Therefore the governing formula for the interpretation of the magnetic anomaly profile becomes b.(x',z,)=fsK(x',,x,z)w(x,z)drdz (6) Equation is a Fredholm integral equation of the first kind, and is an ill-posed problem in the sense of Hadamard's standard 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 0 rather than continuously for all values of one method of solving equation is to convert it into a system of linear equations of the form Ay b (7) where A is an mxn matrix and y and b are nxl and mxl vectors, respectively. There are many numerical quadrature schemes available for approximating an integral equation as a system of linear equations However, mid-point quadrature has been found 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

I

WO 95/20174 PCT/AU95/00039 One such method is Green's strategy To use Green's strategy for converting equation into a linear system of equations, suppose that the area of interest S' is a rectangle given by the Cartesian product S' (8) where is the length of S' and is the width.

Divide S' into JK pixels, each of length Ax and width Az, where A andAz=(d-c) K.

Then the coordinates of the centre of pixel are xi and Zk (k-1/2)Az; k For simplicity, assume that the direction cosines of M and I are P=p=0 and Q=q=1. The mid-point quadrature then defines the discrete kernel function as 2(z, 2 (x',Zo x, x)2 (zk )2 2 With a finite number of measuring points, equation (6) becomes J K ZjAw,; (9) j=1 k=1 1- WO 95/20174 PCT/AU95/00039 6 where N is the number of measuring points, Ak is given by 2(2k 0 2

X,)

2 A, 1AxAz [xi X)+(zk 0 2 2 and w.i is a constant within a pixel but may change from pixel to pixel. In terms of matrices equation can be rewritten as Aw =b where A is a (NxJK) matrix and w and b. are JK and N vectors, respectively.

Now the governing formula has been converted into the system of linear equations 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 bare 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 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 WO 95/20174 PCT/AU9500039 7 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 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 8; 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 8, where the values of 8 are evenly distributed within the range of n] 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 ak, is computed according to a, K(x= Iox, zwv'dxdz I WO 95/20174 PCT/AU95/00039 8 where i denotes a location index, k denotes a strip index, S. is the kth strip of an area of interest, K(x.z,x,z) is a known kernel function, x, is the i 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 inteormediate system of linear equations into a second intermEidiate system of linear equations of the form A*y,=4 where A is a matrix of coefficients obtained from the first intermediate system of linear equations, y, is the i th row of the unknown matrix Y, and b, 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 Ib 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; I o WO 95/20174 PCT/AU95/00039 9 Figure 2 illustrates how an area of interest and a subterranean object is partitioned into strips perpendicular to a direction specified by an angle 08; 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 f in a region S' can be uniquely reconstructed

I

WO 95/20174 PCT/AU95/00039 from line integrals of f 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 x and y, or its polar coordinates r and 0. For an arbitrary straight line in the plane, let 0 denote the angle between the x-axis and the direction normal to the line, and let r denote the signed perpendicular distance from the origin to the line. The line can then be described by the equation xcosO+ysin9= r (11) LetF(r,0) be the Radon transform line integral, or projection of Then f in polar coordinates can be reconstructed uniquely from F(r,8) by the inverse transformation F25d (12) f (77, 0) J, d (12) where is the partial derivative of F(r,9) with respect to r.

In practice, F(T,O) is only known at discrete points along lines, for r=r k =k for nonnegative integer values of k, and 0=0, Here, each value of 6 defines one of N views taken around a semicircle -I WO 95/20174 PCT/AU95/00039 11 surrounding the object, and 1 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 f4, depending on given by F(k(rk, Ojxcos, +ycosO-r,). (13) 2 J=0 k=1 Here, A0 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=1 Image reconstructing techniques based on equation (13) are known to as back projection methods.

Now apply the principles of -cmputerised 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 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 has a magnetic dipole moment density M whose direction is defined by its direction cosines P=cosa and Q=cosp where a is d's angle from the x-axis and fl 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: above the Earth's surface giving an anomaly as illustrated by curve 14 in Fig. lb.

I WO 95/20174 PCT/AU95/00039 12 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 9 with respect to the x-axis, draw parallel lines 16 with the same interval q between two adjacent lines, the normal to all lines making an angle 9 with the x-axis such that they can all be described by equation Let the area between two adjacent lines define a strip. Then q is the width of the strips, and is chosen (to ensure that a desired spatial resolution is obtained in the reconstructed image such that 77min(&,&), where is the size of a pixel to be produced in the reconstructed.

Let SI denote the area of the kth strip, and let Ak denote the intersection of the source body S with S.

Ak=snlS). Then t, J Z, x, z)o(x,z)d i N (14) where t, represents the contribution of the source body from the kth strip to the measured field intensity at the ith measuring point, and N is the total number of measuring points. It follows that ba(i)= ~tik; k=1 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).

WO 95/20174 PCT/AU95/00039 13 Let a denotes the strip integral ai, (16) where w' is a constant and w'=M Since a, and t are numbers, it will always be true that 1& =Yjk% =ykak (17) where ye is an unknown real number and Substituting equation (17) into equation results in

N

ba Z y1kai; i k=1 (18) This can be written in the matrix form Ya=b (19) where Y2 y, yn 2 a, b2 Y= i a= Y. b and, for consistency, the field intensity is denoted by b We now have a system of linear equations (19) in which the unknown is the matrix Y. Unlike equation i WO 95/20174 PCT/AU95/00039 14 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 ith row of Y is computed from equation (19) then, according to the contribution of the source body to the ith measuring point from the kth strip is simply t,=ya. It is then possible to determine the contribution of the source body from all of these strips.

Defining an image function as eS' (21) it is evident that y=ak f(x,z).dz (22) That is, that ya, is the strip integral of the image function 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 9, the set {yeaklk=1,2,...,N} can be considered to be a projection of an image function f(x,z) according to the definition used in computerised tomography as described above.

For each value of 9, there is a set of strips and a corresponding projection of the image function as defined in equation Construct and solve a system of linear equations for each value of 9 given by 9,=jrl/N, to produce N projections of the image function f(x,z) along N directions. With these N projections the image f(x,z) is readily reconstructed through the back projection method of equation Once the image f(x,z) is reconstructed, it is possible to find a WO 9520174 PCT/AU95/00039 the distribution function since the data kernel K(x,,zo,x,z) is known.

Recall from its definition in equation that the data kernel is considered to be a function of the argument while P,Q,p and q are all considered as constants. For a given source body, P and O 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 p 2 (23) the data kernel K can be considered to be a function of the extended argument That is 3[P(x Q(z zo[p(x F z-z) z P, x, z, p) For a given direction given define vector a' as aj a a] where the elements a are computed using WO 95/20174 PCT/AU95/00039 16 a f= -x(xk ,x,.p,)w'dxdz for k (26) For a given 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 [b bbi, b] (27) Then one has Ya' bi; j= (28) Defining the (nxn) matrices A and B by juxtaposing vectors as follows a 2 b 2 (29) enables equation (28) to be written as YA=B. Equation (30) has a unique solution if A is a nonsingular or, equivalently, if the vector set {ajj is a linearly independent set. From the definition of a i it is clear that the linearity of this vector set is determined by the linearity of the data kernel K(x,,z,x,z,pj) with respect to the variable p,.

However, is a nonlinear function of p because, from q= Therefore, are WO 95/20174 PCT/AU95/00039 17 linearly independent vectors, A is non-singular, and has a unique solution given by Y= BA-' (31) where A denotes the inverse of A.

It was indicated before that only one row of the solution matrix is required. If one transposes equation 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'Ye, B'e (33) which gives J= (34) where y, is the ith row of the unknown matrix Y and b, is the ith row of the matrix B.

If we rewrite equation (30) explicitly as 1 Y12 Y' y l a' a) b:l b b S Y2 2 2 y: :22 Y2 a a" b: b, b3"5 yi y. 2 yn, a. b' b2 b then we can see that the elements of the ith row of the matrix B, are the components in the direction p 1 of the total field intensity measured at the same point.

That is WO 95/20174 PCT/AU95/00039 18 b (36) The problem that arises here is how to obtain the vector sets {aJi= and {bij= in practice. The magnetic field intensity measured by most modern threecomponent 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 Y 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, osa' (38) where a' is the angle between C and the x-axis. Let bi=b, (39) where b, is the magnitude of C. We have p= a+(j-)Aa where Aa is the difference between angles pi and and (j-1)Aa is the biasing angle of the pi from It follows that -I U- WO 95/20174 PCT/AU95/00039 19 b, =b,'co(ji-1)Aa] (41) Therefore, b, can be obtained from the magnitude of the total field intensity and direction cosines cos(i-1)Aa] at the measuring point By substituting {pjj=1,2,...,nr into equation the vector set 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 However, the vector set 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 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 socalled damped least squares method in the Marquardt-Levenberg's formulation A'b (42) where A is a damping factor or regularization parameter, and I is the identity matrix. An important issue in the use of this method is the choice of an optimal A, and the generalised cross validation (GCV) method is commonly used for this purpose Quite often, the GCV method fails to find an optimal A, and there now exists a more powerful L-curve method for choosing an optimal A.

WO 95/20174 PCT/AU95/00039 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 surrounding area of interest S' refer to strip numbers within The dimensions of the pixels were set at r=Sz=0.5nxOSm. 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 0 and the total number of views is set to be 33. For the chosen object, with strips in a direction 9054.5°, only strips 28 to 40 intersect with the object, as shown in Fig. 4.

Measurements were simulated at an elevation above the surface of the Earth using equation 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 i in equation (39) was chosen to be i=32. Random errors were added to this perfect data to give b =b 2

+S

32 A corresponding projection of the object-related image function as defined in (21) is shown as line 18 in Fig. 6, with a typical A 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 farntion f(x,z) defined in equation 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 I WO 95/20174 PCTAU95/00039 21 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: it can find a physically acceptable solution to a Fredholm integral equation of the first kind without any a priori knowledge about the solution, its computational cost is low, and 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.

Claims (7)

1. A method for obtaining magnetic dipole density image data concerning a subterranean object S contained in an area of interest 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 clject 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 8; 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 generalised projection; repeating the above three steps for different angles 8, where the values of 8 are evenly distributed within -the range of performing an inverse Radon transformation on the generalised projections to obtain the magnetic dipole moment density image data concerning the object; and displaying a representation of the magnetic dipole moment density image data.

2. The method of claim 1 wherein the strip integral aik is computed according to ak i K(x, z w'dxdz I WO 95/20174 PCT/AU95100039 23 where i denotes a location index, k denotes a strip index, S is the kth strip of an area of interest, K(xz,,x,z) is a known kernel function, x is the i 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.

3. The method of claim 2 wherein 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,=b where A is a matrix of coefficients obtained from the first intermediate system of linear equations, y, is the i th row of the unknown matrix Y, and b is a vector formed from the total magnetic field.

4. The method of claim 3 wherein the step of solving the second intermediate system of linear equations uses a regularization method.

5. The method according to any one of claims i, further comprising the step of displaying a representation of the magnetic dipole moment density image data on display means.

6. The method of claim 1 wherein 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 where Y is an unknown matrix, a is a vector obtained form the strip integrals, and b is I I WO 95/20174 PCT/AU95/00039 24 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 where A is a matrix of coefficients obtained from the first intermediate system of linear equations, y, is the i th row of the unknown matrix Y, and b, is a vector formed from the total magnetic field.

7. The method of claim 6 wherein the step of solving the second intermediate system of linear equations uses a regularization method.