Abstract:
The invention relates to a method of determining the spatial field distribution and/or the spatial position of a useful field source, producing said field distribution, by means of a multichannel field measuring device which comprises a plurality of spatially distributed sensors which produce measured values which contain on the one hand useful measured values stemming from the useful field source and on the other hand noise measured values stemming from at least one noise field source, correction values being formed and superposed on the measured values so that compensated measured values are obtained which are compared with similarly compensated, mathematically derived reference values which are produced by a fictitious reference field source by means of the field sensors, the reference field source being determined for which the pattern of compensated reference values corresponds best to the pattern of compensated measured values.

Description:
BACKGROUND OF THE INVENTION 
     The invention relates to a method of determining the spatial field distribution and/or the spatial position of a useful field source, producing the field distribution by means of a multi-channel field measuring device which comprises a plurality of spatially distributed sensors which generate measured values which contain on the one hand useful measured values stemming from the useful field source and on the other hand noise-measured values stemming from at least one noise field source, correction values being formed and superposed on the measured values so that compensated measured values are obtained which are compared with similarly compensated, mathematically derived reference values which are produced by a fictitious reference field source by means of the field sensors, the method thus determining the reference field source for which the pattern of compensated reference values corresponds best to the pattern of compensated measured values. 
     Using suitable known sensors, such a method enables the measurement of electrical or magnetic field distributions or the localization of associated field sources. For example, magnetoresistive elements constitute suitable magnetic sensors. Sensors utilizing the rotation of the magnetization in ferromagnetic films are also highly sensitive. SQUIDs are the sensors customarily used for measurements in the medical field. 
     The method of the invention is preferably used for the measurement of biomagnetic fields which are produced by currents flowing in the body, notably in the head or the heart. 
     A method of the kind set forth is known from DE-A-43 04 516. 
     Because on the useful fields to be measured there are usually superposed much stronger noise fields from noise field sources which are remote in comparison with the useful field source to be examined, the measured values contain comparatively large noise-measured values in addition to the desired useful measured values. Therefore, steps must be taken so as to suppress the noise measured values. In the known method noise signal vectors are derived from a pattern of measured values determined during a given period of time, said noise signal vectors being used to correct the measured values. It is assumed that the noise signal only slightly, resembles the useful signal. The comparison of the measured values with calculated field values is based on a corrected measured signal vector and a corrected source signal vector. 
     If the use of complex, heavy and expensive shielding devices is to be avoided, the steps taken for noise field compensation must be particularly effective. This holds notably if gradiometer devices are not or cannot be used for medical SQUID measurements, as is the case for example for SQUIDs ceramic high-temperature superconductors. This is because the manufacture of the superconducting connections required for conventional gradiometer devices is problematic. 
     SUMMARY OF THE INVENTION 
     It is an object of the invention to provide a method of the kind set forth whereby the effect of noise fields can be substantially eliminated and whereby not only a spatially uniform noise field is to be detected. 
     This object is achieved in that the pattern of measured values is approximated by a correction function which contains a number of terms of a function series which are essential to the spatial variation of the noise-measured values, in the correction function there being absent at least one term of the function series which is essential to the spatial variation of the useful measured values but not essential to the spatial variation of the noise measured values, that compensated measured values are formed as the difference between the measured values and the transformed measured values occuring at the areas of the associated field sensors, on the basis of the correction function and that the reference values are transformed into compensated reference values in the same way. 
     The method in accordance with the invention enables noise field suppression factors of the order of magnitude of 10 5  to be achieved. The noise field suppression can be performed at any instant, independently of the temporal variation of the noise signals. 
     Converging mathematical function series are suitable for the approximation of the measured pattern of measured values by a correction function. 
     Particularly suitable are power series, but in special cases other function series may also be attractive, depending on the spatial arrangement of the sensors and also on the characteristic of the noise fields, for example Tschebyscheff polynomials, trigonometric functions or spherical functions, which are particularly suitable for helm systems for measurement of cerebral fields. 
     The number of sensors used to form the correction function should be at least &#34;1&#34; larger than the total number of function coefficients to be determined. A redundant set of equations is then obtained and correction function coefficients offering optimum adaptation of the measured values are determined (linear least squares fit). The larger the number of measured values taken into account from different sensors, the less the intrinsic sensor noise will affect the result. 
     Evidently, the spatial position of the sensors taken into account should be such that they are capable of supplying linearly independent information for the mathematical determination of the function coefficients. With a view to achieving a high accuracy, the information should be significantly independent. This means in practice that the sensors used should be distributed as uniformly as possible across the space available for the measuring device and should be situated as far as possible from one another. If field components or field vectors pointing in a plurality of coordinate directions are to be determined, a corresponding number of sensors must be assigned to each coordinate direction. 
     In a simple version of the method in accordance with the invention the field sensors are arranged in a single measuring plane so that field components extending perpendicularly to the measuring plane are determined. 
     Further field components can be measured by means of additional field sensors which are arranged in measuring planes which are inclined preferably orthogonal to the first measuring plane. 
     Because field sensors cannnot be arranged in a measuring plane with absolute precision due manufacturing tolerances, field components from directions orthogonal to the desired measuring direction could also be detected in the case of angular deviations. Correction is then possible in that with the field sensors arranged in a measuring plane there are associated correction sensors for measuring field components extending parallel to the measuring plane, and that via correction circuits the measured values from field sensors unintentionally inclined relative to the measuring plane are corrected by means of the measured values obtained by the correction sensors. 
     Evidently, the method of the invention can be used analogously when measured values are present in the form of differences between field strengths detected at a distance from one another, for example by means of gradiometer devices. Field sensors are then arranged in a plurality of mutually parallel planes in order to measure field gradients. 
     Because of the large noise component in comparison with the actual useful signal, the sensor signals have a very high level. For better use of the dynamic range of subsequent amplifiers, notably A/D converters, it is particularly attractive to choose the measured values to be the differences between the values measured by means of the field sensors and by means of a coaxial level sensor, and to reduce the reference values also by the value of the correction function at the location of the level sensor. The differences on the output of the sensors should be formed by analog circuits, notably by superposing the output voltage of the level sensor on the output voltages of the measuring sensors. 
     The object of the invention is also achieved by means of a measuring device for determining the spatial field distribution and/or the spatial position of a useful field source producing the field distribution, which device consists of: 
     a multichannel field measuring device which comprises a plurality of spatially distributed sensors whereby measured values are produced which contain on the one hand useful measured values stemming from the useful field source and on the other hand noise measured values stemming from at least one noise field source, 
     a correction circuit for forming correction values which are superposed on the measured values, so that compensated measured values are formed as the difference between the measured values and the values occurring at the areas of the associated field sensors because of a correction function, 
     an arithmetic circuit for determining reference values produced by a fictitious reference field source by means of the field sensors, and for compensating these reference values in the same way as in the correction circuit, thus producing compensated reference values, 
     a comparison circuit for comparing the compensated measured values with the compensated reference values, the reference field source thus being determined for which the pattern of compensated reference values corresponds best to the pattern of compensated measured values, the pattern of measured values being approximated by a correction function which contains a number of terms of a function series which are essential to the spatial variation of the noise measured values, in the correction function there being absent at least one term of the function series which is essential to the spatial variation of the useful measured values but not essential to the spatial variation of the noise-measured values. 
     These and other aspects of the invention will be apparent from and elucidated with reference to the embodiments described hereinafter. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     In the drawings: 
     FIG. 1 shows a flow chart illustrating the method of the invention, 
     FIG. 2 shows diagrammatically an array of sensors arranged at the same angle in a measuring plane, 
     FIG. 3 shows sensors arranged in three mutually parallel planes, and 
     FIG. 4 shows a circuit diagram for the analog subtraction of the signal of a level sensor from the signals of the other sensors. 
    
    
     DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 
     The overall execution of the method of the invention will be described in detail hereinafter with reference to FIG. 1. 
     The field distribution of a useful field N, on which a large noise field S is superposed, is to be measured in a two-dimensional or three-dimensional zone by means of a measuring device 1. The measuring device 1 is constructed, for example as described in Appl. Phys. Lett. 63, 1963, pp. 2700-2702 and comprises a plurality of sensors in the form of spatially distributed DC SQUID magnetometers and also comprises the necessary auxiliary devices such as amplifiers and the like. 
     Biomagnetic sources which are situated at a distance of from approximately 5 cm to 10 cm from the measuring device produce flux densities of approximately 100 fT at the areas of the sensors. Noise fields stem from noise sources which are much more remote and cause flux densities of approximately 100 pT. 
     The individual sensors supply individual measured values Mi, the totality of which is denoted by the letter M. This concerns a pattern of individual measured values determined in a two-dimensional or three-dimensional zone. 
     Because of the substantially different distances between on the one hand the useful field source and the measuring device and between the measuring device and the noise field source on the other hand, the spatial characteristics of the useful measured values and the noise-measured values significantly deviate. When the useful measured values and the noise measured values are adapted by way of a respective, same mathematical function series (notably a power series), it appears that for some terms the associated function coefficients are pronounced upon adaptation of the useful measured values, whereas for the adaptation of the noise measured values the corresponding coefficients of these terms can be practically ignored. 
     The measured values M are adapted by a correction function M&#39; in which at least one term of the selected function series is absent, the coefficient of said term being of significance only for the detection of the useful measured values but insignificant for the detection of the noise measured values. M&#39; is then obtained, using an arithmetic circuit 2, as the sum of a constant value Mo and a sum of the products of selected function terms Fm and the associated coefficients Mm. 
     An arithmetic circuit 3 forms compensated measured values Mk as differences between the original measured values M and the values M&#39; calculated from the correction function. The values Mk are applied to a comparison circuit 4 whereby, the useful field distribution and/or the spatial position of the associated useful field source is calculated on the basis of the values Mk. To this end, the arithmetic circuit 5 calculates reference values V of various fictitious useful field sources, which reference values were measured in the respective sensors. The pattern of these fictitious reference values V is adapted, in the same way as the pattern of measured values M, by a correction function V&#39; as the sum of a constant value Vo and the sum of the product of the term Fm chosen for the measured values M and the associated coefficients Vm (arithmetic circuit 6). In the arithmetic circuit 7 compensated reference values Vk are formed from the difference between the calculated reference values V and the transformed values V&#39; calculated for the same measuring locations. The comparison circuit 4 compares the values Mk and Vk. In the case of deviations, the procedure is repeated, using new reference values V from another, fictitious useful field source, until optimum correspondence is obtained. Subsequently, the data of the optimum reference source can be applied to an output device 8 as the data of the useful field source. 
     The quantities mentioned with reference to FIG. 1 may be vectors or scalar components. 
     Referring to FIG. 2, 19 z-magnetometers 9 are coaxially arranged in a measuring plane 10. 19 different inductions Bz can thus be measured in the z direction, extending perpendicularly to the measuring plane 10. In addition, correction magnetometers 11 and 12 are arranged at the area of the central z magnetometer 13, said correction magnetometers enabling measurement of inductions Bx and By in the x direction and the y direction, respectively, parallel to the measuring plane. 
     The z magnetometers (numbering n+1=19) measure inductions Bi=Bo to Bn at the locations having the coordinates x i  and y i  in the measuring plane 10. 
     For the correction function, a Taylor power series terminated after the quadratic term was found to be suitable. Its expansion coefficients G k  =G 0  to G 5  can be calculated by way of the following set of equations: ##EQU1## in vector notation this set of equations becomes: ##EQU2## 
     The coefficients G k  are determined by solution of the set of equations. For the method of the invention, n is preferably substantially larger than 6. The system is then redundant and the coefficients G k  are determined which offer optimum adaptation the measured values B i  (linear least squares fit). The larger the number of magnetometers, the less the intrinsic sensor noise will affect the result. 
     
         G=B*X.sup.T (X*X.sup.T).sup.-1 
    
     Subsequently, the noise signal calculated by means of this approximation at the location of the i th  magnetometer is subtracted from the signal measured at the location i: 
     
         B*.sub.i =B.sub.i -(G.sub.0 +G.sub.1 x.sub.i +G.sub.2 y.sub.i +G.sub.3 x.sub.i y.sub.i +G.sub.4 x.sub.i x.sub.i +G.sub.5 y.sub.i y.sub.i) 
    
     
         B=B-G*X 
    
     It is thus achieved that not only the mean magnetic field background and the mean gradient of the first order are subtracted, but also that a quadratic curvature of the magnetic field is compensated. Because of the spatial proximity of the sources to the magnetometers, the biomagnetic signal exhibits a strong curvature (terms of third and higher order) in the variation of the magnetic field which is extracted as the actual measured signal. 
     In order to ensure that the dynamic range of the electric circuits succeeding the magnetometers is not exceeded, the signal B o  of the central z magnetometer 13 is subtracted from all other signals B 1  to B n . This can be attractively performed by means of a small coupling-in coil connected to each other sensor (FIG. 4). These coils are usually present any way because they serve to drive the magnetometers as zero detectors in conjunction with a fed back control circuit. An analog subtraction circuit preceding each A/D converter 15 is also feasible. There may also be reasons for the signal of the central z-magnetometer to be generated as an overall external counter field by a large coil system, so that the same effect is obtained. 
     In order to compensate the inclination error of the magnetometers, this error is first measured exactly. This measurement can be by means of a system of Helmholtz coils whereby uniform field are applied successively in the x direction, the y direction and the z direction. The three measured signals in [V], divided by the magnetic field strength in (T), can be considered to be components of a normal vector of the magnetometer, the length of the vector indicating the sensitivity: Q ix , Q iy , C iz , (Q denotes &#34;transverse sensitivity&#34;). The sensitivity of the x magnetometer and the y magnetometer in the x direction and the y direction is assumed to be C ox  and C oy , respectively. 
     The vector of the magnetic field in this location is measured by means of the central x, y and z magnetometers in the proposed device: B ox  =S ox  /C ox , B oy  =S oy  /C oy , B oz  =S oz  /C oz . The signals of all other magnetometers S i  [V] can then be corrected as follows (the analog subtraction of the signal from the central z-compensation magnetometer S oz  /C oz  as described above was also taken into account) herein: 
     
         B.sub.irci =S.sub.i /C.sub.iz -(S.sub.0x /C.sub.0x)*(Q.sub.ix /C.sub.iz)-(S.sub.0y /C.sub.0y)* (Q.sub.iy /C.sub.iz)-S.sub.0z /C.sub.0z 
    
     It is assumed herein that the noise field at the location of the individual magnetometers approximately equals the field at the location of the central three-axis magnetometer. This is, of course, a coarse approximation only. The coefficients Q ix  and Q iy , however, are only 1/100 of C iz  when the device is accurately constructed, so that they represent only a correction quantity. Therefore, this approximation can be used in this case. 
     When a high-precesion construction with extremely small inclination errors of the z magnetometer is achieved, the correction by the x and y magnetometers can also be dispensed with. 
     If a more accurate correction of the inclination errors is required, more accurate correction functions for the x and y field components can be formed by means of a plurality of x and y magnetometers per measuring plane so as to represent their spatial variation more accurately. 
     Evidently, all manipulations performed on the measured signal must be taken into account in the reconstruction algorithm. This is comparatively simply possible; because an analytic solution of the inverse problem does not exist, always the &#34;forward problem&#34; is solved and the calculated field distribution is adapted to the measured field distribution. In the &#34;forward problem&#34;, however, all corrections described above can be exactly implemented. 
     In a preferred embodiment of the device as shown in FIG. 3 the derivative of the magnetic fields in the z direction is also taken into account. This is of importance if not all biomagnetometers are arranged in one plane. The line 16 in FIG. 3 is intended to indicate that the z magnetometers 17 are distributed across a curved surface. 
     A plurality of magnetometers are provided for the z direction, and a respective correction magnetometer 18, or 19 is provided for the x direction and the y direction, respectively. An arbitrary number of biomagnetometers 17 in the z-direction are enclosed by compensation magnetometers 20, 21 and 22 which are arranged in different planes which are situated at a distance from one another in the z- direction. These z-compensation magnetometers are arranged at the corners of three large isosceles triangles situated at a distance of each time from 50 mm to 100 mm. In each plane z- compensation gradiometers are arranged on an isosceles triangle. A tenth z-compensation magnetometer 23 is arranged, together with the x and y correction magnetometers, at the centre of the triangle in the central plane. 
     The set of equations for determining the coefficients of the transformed function is now: ##EQU3## or in the simpler vector notation: ##EQU4## This set of equations can be solved exactly. The individual measured values of the biomagnetometers are then corrected as follows: ##EQU5## 
     The above correction of the inclination angle of the magnetometers by the magnetometers in the x-direction and the y-direction and the improved use of the dynamic range by analog subtraction of a central z-magnetometer signal could also take place. 
     For the described arrangement or for other feasible arrangements of magnetometers it is only important that the coefficient matrix X is not singular. Moreover, the information contained in this matrix should be substantially linearly independent so that the coefficients G i  can be calculated in a stable fashion. This is case notably for the described arrangements. 
     In an environment which is less disturbed, the Taylor series can be terminated already after the linear term. 
     It is not very decisive whether the magnetometers used for calculating the correction coefficients themselves have a large biomagnetic signal or not. It is to be ensured that the correction does not lead to excessive loss of biomagnetic signal. As a measure in this respect use can be made of the mean signal amplitude of a current dipole (typical biomagnetic source): ##EQU6## 
     The mean amplitude for near current dipoles should not drop to values smaller than approximately 10% after the noise field compensation. This is the case in the described devices. 
     An analogous comparison of the calculated amplitude of a disturbing magnetic dipole at a distance of a few meters before and after the compensation has revealed that a power series of up to the second order yields noise field suppression factors of 10 5 . 
     Other functions can also be used as polynomials for a series expansion. It is only important that the termination of the series expansion after one of the first terms leads to spatial high-pass filtering. 
     The correction of noise signals by a power series as described above can also be used for devices comprising conventional gradiometers. 
     The above correction formules also apply when the output signals of the gradiometers are inserted for B. The transverse sensitivity of the gradiometers can be corrected accordingly, by means of x and y magnetometers. 
     Analogously to the described method, magnetic field components other than the z-component described above can also be determined. For example, in many locations each time three orthogonal magnetometers may be provided. In that case the z as well as the x and y components will be expanded in a power series. The measured values are again appropriately corrected by subtraction of the linear and/or quadratic term. The locations for the three-axis magnetometers must be situated so that none of the coefficient matrices is singular. Such arrangements are important for systems in which the biomagnetic signals are to be recorded around the entire head (helm system) or around the torso, the magnetometers being uniformly distributed around the head or the torso. When the x, y and z components of the magnetic field are treated completely independently, for the same noise field suppression factors, the number of magnetometers is tripled in comparison with the measurement in only one coordinate direction. 
     Three-dimensional systems of this kind can be further simplified by strictly taking into account the relations between the components of a flux density vector field which follow from the Maxwell equations. 
     When the correction function is formed up to the quadratic terms of a Taylor series, it can be demonstrated that 15 linearly independent equations can be indicated for the Taylor coefficients to be determined in the correction algorithm, so that the number of the magnetometers required to achieve for the same noise field suppression can be reduced by 15. 
     The linearly independent equations for the Taylor coefficients are obtained as folows. 
     Formulation for the Taylor expansion: 
     
         B.sub.x (x,y,z)=A.sub.0 +A.sub.1 x+A.sub.2 y+A.sub.3 z+A.sub.4 x.sup.2 +A.sub.5 y.sup.2 +A.sub.6 z.sup.2 +A.sub.7 xy+A.sub.8 yz+A.sub.9 xz 
    
     
         B.sub.y (x,y,z)=B.sub.0 +B.sub.1 x+B.sub.2 y+B.sub.3 z+B.sub.4 x.sup.2 +B.sub.5 y.sup.2 +B.sub.6 z.sup.2 +B.sub.7 xy+B.sub.8 yz+B.sub.9 xz 
    
     
         B.sub.z (x,y,z)=C.sub.0 +C.sub.1 x+C.sub.2 y+C.sub.3 z+C.sub.4 x.sup.2 +C.sub.5 y.sup.2 +C.sub.6 z.sup.2 +C.sub.7 xy+C.sub.8 yz+C.sub.9 xz 
    
     Calculation of the derivatives: ##EQU7## Analogous equations hold for B y  and B z . It follows from the Maxwell equation div B=0 that: 
     
         A.sub.1 +2A.sub.4 x+A.sub.7 y+A.sub.9 z+B.sub.2 +2B.sub.5 y+B.sub.7 x+B.sub.8 z+C.sub.3 +2C.sub.6 z+C.sub.8 y+C.sub.9 x=0 ∀x,y,z 
    
     The following linearly independent equations result therefrom: 
     
         A.sub.1 +B.sub.2 +C.sub.3 =0 
    
     
         2A.sub.4 +B.sub.7 +C.sub.9 =0 
    
     
         A.sub.7 +2B.sub.5 +C.sub.8 =0 
    
     
         A.sub.9 +B.sub.3 +2C.sub.6 =0 
    
     From the Maxwell equation rot B=0 (no currents flow outside the body!) it also follows that: 
     
         (B.sub.3 +2B.sub.6 z+B.sub.8 y+B.sub.9 x)-(C.sub.2 +2C.sub.5 y+C.sub.7 x+C.sub.8 z)=0 
    
     
         (C.sub.1 +2C.sub.4 x+C.sub.7 y+C.sub.9 z)-(A.sub.3 +2A.sub.6 z+A.sub.8 y+A.sub.9 x)=0 
    
     
         (A.sub.2 +2A.sub.5 y+A.sub.7 x+A.sub.8 z)-(B.sub.1 +2B.sub.4 x+B.sub.7 y+B.sub.9 z)=0 ∀x,y,z 
    
     Further linearly independent equations follow therefrom: 
     
         B.sub.3 -C.sub.2 =0 C.sub.1 -A.sub.3 =0 A.sub.2 -B.sub.1 =0 
    
     
         B.sub.9 -C.sub.7 =0 2C.sub.4 -A.sub.9 =0 A.sub.7 -2B.sub.4 =0 
    
     
         2B.sub.6 -C.sub.8 =0 C.sub.9 -2A.sub.6 =0 A.sub.8 -B.sub.9 =0 
    
     Of the three equations denoted by *, only two are linearly independent, so that the third is identically solved. All other equations are linearly independent, so that the number of linearly independent equations is 15. 
     If the measured magnetic field is to be expanded in a Taylor series up to the second order as described above, the following set of equations must be solved: ##EQU8## 
     For example, if B x , B y  and B z  are measured in at least five suitable points, the coefficients a . . . p can be calculated. It is again advantageous to perform measurements in substantially more than five points and to determine the coefficients a . . . p which provide optimum adaptation of the measured values. To this end, use can be made of the following method: ##EQU9## n j  are the normal vectors of the magnetometers: 
     
         n.sub.i =(n.sub.xi,n.sub.yi,n.sub.zi).sup.T and X.sub.i =(1,x.sub.i,y.sub.i,z.sub.i,x.sub.i.sup.2 . . . ).sup.T 
    
     The derivatives ##EQU10## then provide the desired coefficients a . . . p of the matrix U. All measured values can then be corrected as follows: 
     
         B.sub.i.sup.corr =B.sub.i.sup.ori -U*X.sub.i 
    
     Another attractive embodiment (for example, for the so-called &#34;helm system&#34;) comprises measuring sensors which are uniformly distributed across a strongly curved surface and which determine the components of the magnetic field perpendicularly to the surface. It is to be noted that the measuring sensors provide measured values which contain contributions from the field components in the x and/or y and/or z direction. This necessitates more work calculation and complicates the set of equations for determining the function coefficients.