# Experimental Study on the Aperture of Geomagnetic Location Arrays

Authors: Xiaojun Zhang; Xiyuan Kang; Xin Chen; Zhuoshan Geng; Liming Fan; Quan Zheng; Hua Lv; Chong Kang
Journal: Journal of Sensors
 (2019)

---

## Abstract

A method of locating a magnetic target based on geomagnetic total field is proposed. In the method, a conjugate gradient algorithm is introduced to eliminate the time-varying and uneven spatial distribution of geomagnetic total field. Then a structure of the measuring array of geomagnetic total field is designed. In the measuring array, the array aperture is a primary factor for the conjugate gradient algorithm. To determine an optimal aperture, we analyze the relationship between the array aperture and the localization accuracy. According to the localization theory based on geomagnetic total field, we simulate the process of determining an optimum array aperture. Based on the simulation, we propose the basis and principle of determining the optimum array aperture. To prove it, we use optically pumped magnetometers with different array apertures to carry out the experiments of locating a car in a suburb. Through the experiment, we get the experimental relationship between apertures and location accuracy. And the relationship agrees with the theory. The result shows that the method is feasible to determine the optimum aperture.


---

## Body

## 1. Introduction

It is important to locate targets by magnetic field in geological monitoring, energy and mineral exploration, rescue of plane crash, antisubmarine detection, and medical diagnosis [1–3]. In measuring geomagnetic total field, optically pumped magnetometers have the advantages of high resolution, long detection range, and no temperature drift [4–6]. Therefore, it is feasible to locate the magnetic target by using the array of optically pumped magnetometer. In order to improve the location accuracy, the time-varying and uneven spatial distribution of geomagnetic total field must be eliminated [7, 8]. In this paper, a double gradient algorithm is introduced to eliminate the distribution. In the method, some magnetic field changes caused by targets have also been filtered out. In particular, when two sensors are close to each other, the magnetic field of targets will be filtered out more. On the one hand, all the information is missing when the array aperture is zero. On the other hand, the disturbance field cannot be filtered when the array aperture is infinite. Therefore, it is necessary to determine an optimal array aperture which filters the interference field to the maximum extent and filters the target field to the minimum extent.The double gradient algorithm comprehensively describes the actual situation of target’s magnetic field gradient, so the algorithm is not affected by the filtering of some signal gradient. The difference is only that the measured value of the decimal point has been moved forward. For the sensors with high precision and high resolution, the accuracy of locating a target will not be affected as long as the mantissa can be measured. However, if the sensor accuracy and resolution are low, the mantissa will not be detected. It means that the dual gradient algorithm cannot detect the remote or small target. An optimum aperture is an equilibrium point of the double gradient algorithm. Base on the optimum aperture, we can filter the interference field to the maximum extent and filter the target field to the minimum extent. It is impossible to eliminate the magnetic field generated by the target completely. It is only required that the desired accuracy of position is within an expected range.The optimum array aperture is related to sensor precision, resolution, noise of instrument and environment, target magnetic moment, detection position, and detection range. These quantities should be estimated and expected before designing the positioning system, so the calculation method of the optimum aperture can be given according to these estimates.


## 2. Method of Determining an Array Aperture

### 2.1. The Location Algorithm

When the distance between a target and a sensor is more than 2-3 times the scale of the target, the target can be regarded as a magnetic dipole [9–11]. To measure the magnetic field gradient of a target in a 3D space and locate it, the array can be composed of several general field sensors as shown in Figure 1. The distance D between a sensor and origin O is defined as the array aperture.Figure 1
Schematic diagram of the location array.T0 is the magnetic field vector with no target, T0 is its magnitude, and the local magnetic dip and magnetic declination are, respectively θ and φ, so the unit direction vector of T0 is e→=cosθ⋅cosφcosθ⋅sinφsinθ. For i=1,2,3…, Ti is the measurement of the magnetic sensor i at xi,yi,zi with a magnetic target. The target is at x,y,z. Then, the displacement vector from the target to the sensor i is r⇀i=xi−xi⇀+yi−yj⇀+zi−zk⇀, and ri is its magnitude. Pm is the magnetic moment vector of the target. Pm is its magnitude. α and β denote the offset angle and the tilted angle of Pm, respectively. The magnetic field vector Bi denotes the field generated by the target at the ith magnetic sensor. In the far-field condition, Bi is much less than T0, so the total field sensor measurement is Ti≈T0+e→⋅Bi [12].According to the magnetic dipole model in far field, the magnetic fieldBi generated by the target can be expressed as
(1)Bi=μ04πri33Pm⋅r→ir→iri2−Pm.Substituting(2)Pm=Pmcosαcosβcosαsinβsinα,(3)ri=x−xiy−yiz−zi,into equation (1), Bi can be expressed as
(4)Bi=μ0Pm4πri33x−xicosαcosβ+y−yicosαsinβ+z−zisinαri2x−xiy−yiz−zi−cosα⋅cosβcosα⋅sinβsinα.Then, the scalar measurement of theith magnetic sensor can be expressed as
(5)Ti=T0+μ0Pm4πri3⋅3x−xicosαcosβ+y−yicosαsinβ+z−zisinαri×x−xicosθcosφ+y−yicosθsinφ+z−zisinθri−cosθ⋅cosφcosα⋅cosβ+cosθ⋅sinφcosα⋅sinβ+sinθsinα.It can be seen that formula (5) is a scalar equation about the variables x, y, z, Pm, α, and β. x,y,z denotes the spatial position of the target. Pm can be used to estimate the size of the target. Due to high nonlinear function, it is difficult to get an analytical solution of the variables. Generally, we can obtain a numerical solution by solving the optimization problem. Therefore, we solve equation (5) fast and accurately using the software LINGO [13].


### 2.2. Double Gradient Algorithm

The time-varying and uneven spatial distribution may have an influence on locating a target. The influence must be eliminated in magnetic field measurement. The measured value ofith sensor located at xi,yi,zi is Tit0,xi,yi,zi at time t0 and Tit,xi,yi,zi at time t. The measurements of jth sensor are expressed in the same way. ΔTij is a double gradient of the time and space of geomagnetic measurements. The double gradient ΔTij is expressed as
(6)ΔTij=Tit,xi,yi,zi−Tjt,xj,yj,zj−Tit0,xi,yi,zi−Tjt0,xj,yj,zj.Through (5) and (6), we can locate the target location by ΔTij. In this way, the time-varying and uneven spatial distribution can be filtered out [12]. Therefore, we use double gradient ΔTij in the whole experiment. The magnetic field generated by the target is also filtered to some extent when using the double gradient. Therefore, the optimum array aperture is necessary to filter the interference field to the maximum extent and the target field to the minimum extent.


### 2.3. Principle of Determining the Optimum Aperture

In the linear motion of a target, its magnetic moment vector is fixed, which can be regarded as unknown constants. According to formula (6), ΔTij is a function of variables x,y,z,D. It can be expressed as
(7)ΔTij=fx,y,z,D.The uncertainty atx,y,z is expressed as
(8)ΔΔTij=∂f∂x⋅Δx+∂f∂y⋅Δy+∂f∂z⋅Δzx,y,z.In the formula,ΔΔTij cannot be less than the instrument and ambient noise ΔTmin. It means that there is a lower limit ΔTmin. Δx,Δy,Δz denote the positioning accuracy of the target in three directions. It is assumed that all the positioning accuracy in three directions is Δr, then
(9)ΔΔTij=∂f∂x+∂f∂y+∂f∂zx,y,z⋅Δr=gx,y,z,D⋅Δr≥ΔTmin.In order to improve the positioning accuracy, the functiongx,y,z,D should be extremely high. When the target is at a certain point x,y,z, gx,y,z,D is only a function of D. Therefore, there is an optimal value of the aperture D, which makes gx,y,z,D reach a maximum. It means that, when ∂gx,y,z,D/∂D=0, the optimal value of the theory is obtained. It can be expressed as
(10)∂2ΔTij∂x∂D+∂2ΔTij∂y∂D+∂2ΔTij∂z∂D=0.According to the above principle, the optimal aperture can be determined by the following procedures:(1)
Determinate the instrument and ambient noiseΔTmin
(2)
Set the location accuracyΔr of the experiment
(3)
For specificx,y,z, solve for a set of D in the function gx,y,z,D≥ΔTmin/Δr
(4)
Find the intersection ofD sets corresponding to all x,y,z points in the location space
(5)
Choose the minimum value ofD in the intersection to minimize the scale of the arrayStep 4 can be omitted when a coordinate point or a target in a local area is concerned.In theory, it is not possible in theory to find the correspondingD value for all x,y,z in the location space. Therefore, it can be calculated by dividing grid points in the measurement space, and the spacing of grid points can be set according to the locating requirement, so that the computation workload can be reduced. In practice, we can calculate the optimum aperture of several representative points in the region, or the three-dimensional space can be reduced to two dimensions or one dimension. Then the computation can be further reduced.


## 2.1. The Location Algorithm

When the distance between a target and a sensor is more than 2-3 times the scale of the target, the target can be regarded as a magnetic dipole [9–11]. To measure the magnetic field gradient of a target in a 3D space and locate it, the array can be composed of several general field sensors as shown in Figure 1. The distance D between a sensor and origin O is defined as the array aperture.Figure 1
Schematic diagram of the location array.T0 is the magnetic field vector with no target, T0 is its magnitude, and the local magnetic dip and magnetic declination are, respectively θ and φ, so the unit direction vector of T0 is e→=cosθ⋅cosφcosθ⋅sinφsinθ. For i=1,2,3…, Ti is the measurement of the magnetic sensor i at xi,yi,zi with a magnetic target. The target is at x,y,z. Then, the displacement vector from the target to the sensor i is r⇀i=xi−xi⇀+yi−yj⇀+zi−zk⇀, and ri is its magnitude. Pm is the magnetic moment vector of the target. Pm is its magnitude. α and β denote the offset angle and the tilted angle of Pm, respectively. The magnetic field vector Bi denotes the field generated by the target at the ith magnetic sensor. In the far-field condition, Bi is much less than T0, so the total field sensor measurement is Ti≈T0+e→⋅Bi [12].According to the magnetic dipole model in far field, the magnetic fieldBi generated by the target can be expressed as
(1)Bi=μ04πri33Pm⋅r→ir→iri2−Pm.Substituting(2)Pm=Pmcosαcosβcosαsinβsinα,(3)ri=x−xiy−yiz−zi,into equation (1), Bi can be expressed as
(4)Bi=μ0Pm4πri33x−xicosαcosβ+y−yicosαsinβ+z−zisinαri2x−xiy−yiz−zi−cosα⋅cosβcosα⋅sinβsinα.Then, the scalar measurement of theith magnetic sensor can be expressed as
(5)Ti=T0+μ0Pm4πri3⋅3x−xicosαcosβ+y−yicosαsinβ+z−zisinαri×x−xicosθcosφ+y−yicosθsinφ+z−zisinθri−cosθ⋅cosφcosα⋅cosβ+cosθ⋅sinφcosα⋅sinβ+sinθsinα.It can be seen that formula (5) is a scalar equation about the variables x, y, z, Pm, α, and β. x,y,z denotes the spatial position of the target. Pm can be used to estimate the size of the target. Due to high nonlinear function, it is difficult to get an analytical solution of the variables. Generally, we can obtain a numerical solution by solving the optimization problem. Therefore, we solve equation (5) fast and accurately using the software LINGO [13].


## 2.2. Double Gradient Algorithm

The time-varying and uneven spatial distribution may have an influence on locating a target. The influence must be eliminated in magnetic field measurement. The measured value ofith sensor located at xi,yi,zi is Tit0,xi,yi,zi at time t0 and Tit,xi,yi,zi at time t. The measurements of jth sensor are expressed in the same way. ΔTij is a double gradient of the time and space of geomagnetic measurements. The double gradient ΔTij is expressed as
(6)ΔTij=Tit,xi,yi,zi−Tjt,xj,yj,zj−Tit0,xi,yi,zi−Tjt0,xj,yj,zj.Through (5) and (6), we can locate the target location by ΔTij. In this way, the time-varying and uneven spatial distribution can be filtered out [12]. Therefore, we use double gradient ΔTij in the whole experiment. The magnetic field generated by the target is also filtered to some extent when using the double gradient. Therefore, the optimum array aperture is necessary to filter the interference field to the maximum extent and the target field to the minimum extent.


## 2.3. Principle of Determining the Optimum Aperture

In the linear motion of a target, its magnetic moment vector is fixed, which can be regarded as unknown constants. According to formula (6), ΔTij is a function of variables x,y,z,D. It can be expressed as
(7)ΔTij=fx,y,z,D.The uncertainty atx,y,z is expressed as
(8)ΔΔTij=∂f∂x⋅Δx+∂f∂y⋅Δy+∂f∂z⋅Δzx,y,z.In the formula,ΔΔTij cannot be less than the instrument and ambient noise ΔTmin. It means that there is a lower limit ΔTmin. Δx,Δy,Δz denote the positioning accuracy of the target in three directions. It is assumed that all the positioning accuracy in three directions is Δr, then
(9)ΔΔTij=∂f∂x+∂f∂y+∂f∂zx,y,z⋅Δr=gx,y,z,D⋅Δr≥ΔTmin.In order to improve the positioning accuracy, the functiongx,y,z,D should be extremely high. When the target is at a certain point x,y,z, gx,y,z,D is only a function of D. Therefore, there is an optimal value of the aperture D, which makes gx,y,z,D reach a maximum. It means that, when ∂gx,y,z,D/∂D=0, the optimal value of the theory is obtained. It can be expressed as
(10)∂2ΔTij∂x∂D+∂2ΔTij∂y∂D+∂2ΔTij∂z∂D=0.According to the above principle, the optimal aperture can be determined by the following procedures:(1)
Determinate the instrument and ambient noiseΔTmin
(2)
Set the location accuracyΔr of the experiment
(3)
For specificx,y,z, solve for a set of D in the function gx,y,z,D≥ΔTmin/Δr
(4)
Find the intersection ofD sets corresponding to all x,y,z points in the location space
(5)
Choose the minimum value ofD in the intersection to minimize the scale of the arrayStep 4 can be omitted when a coordinate point or a target in a local area is concerned.In theory, it is not possible in theory to find the correspondingD value for all x,y,z in the location space. Therefore, it can be calculated by dividing grid points in the measurement space, and the spacing of grid points can be set according to the locating requirement, so that the computation workload can be reduced. In practice, we can calculate the optimum aperture of several representative points in the region, or the three-dimensional space can be reduced to two dimensions or one dimension. Then the computation can be further reduced.


## 3. Simulation

### 3.1. Simulation of the Experimental Environment

As shown in Figure2, four optically pumped magnetometers T1,T2,T3,T4 form a square array.Figure 2
Schematic diagram and photo of the experimental area.
(a)
(b)The aperture of the array isD. The target car moves at a constant velocity from y=32.8m to y=−41m along a parallel line 32.02 m from the y axis. Local magnetic dip is 63.3°, and local magnetic declination is 10.34°.The curve of the magnetic field changed over time in the experimental area, as shown in Figure3. The horizontal axis represents time t, and the unit is s. The vertical axis is the double gradient value of magnetic field, and the unit is nT. The array aperture is D=6m. The red line in Figure 3 is the measurement value T1 of sensor 1 (for the convenience of making the drawing, each value minus 55315 nT). T1 includes the geomagnetic field which varies with time and the magnetic field generated by the target car which varies by car’s movement. Because the geomagnetic field itself varies greatly over time, the red line cannot be used to judge when a car passed. The blue line is the data of ΔT14, which is processed by the double gradient of formula (6). The time-varying and uneven spatial distribution of the geomagnetic field is filtered using the double gradient of the magnetic field. The sensitivity of sensor CS-L is 0.6 pT/Hz1/2@1 Hz (rms). The peak value of the instrument noise is 2 pT, and the bandwidth is 0.1 Hz. The measurement noise generated by the double gradient of local environmental magnetic field is 10 pT. The measurement value of the car magnetic moment is Pm=405Am2, α=0.592rad, and β=3.74rad.Figure 3
Curve of magnetic field variation over time in experimental area.


### 3.2. Results of the Simulation

#### 3.2.1. Results with the Optically Pumped Magnetometers

Figure4 shows the curve of the partial derivative ∂ΔT34/∂y and D in simulation that y=−40m. The horizontal axis is D value, unit m, value range (0 m, 100 m). The vertical axis is the value of ∂ΔT34/∂y, unit nT/m, variation range (0 nT/m, 0.1 nT/m). When D=32m, the curve has a maximum value. Thus, the theoretical optimum aperture D is 32 m, and the location accuracy of the target is the highest. The noise of environment and instrument is 10 pT. If the location accuracy is 1 m, it is required that ∂ΔT34/∂y≥10pT/m in the curve of Figure 4. In fact, D∈9,43.5∪47,91m satisfies the location accuracy, so Dmin=9m is chosen to minimize the array scale. Figure 5 is the curve of the second mixed partial derivative ∂2ΔT34/∂y∂D and D in simulation that y=−40m. The horizontal axis is D value, unit m, value range (0 m, 100 m). The vertical axis is the value of ∂2ΔT34/∂y∂D, unit nT/m2, variation range (0 nT/m2, 0.1 nT/m2). There are two D values corresponding to the second derivative ∂2ΔT34/∂y∂D=0, 32 m, and 59 m. There is a maximum value for the first derivative ∂ΔT34/∂y at D=32m.Figure 4
Curves of∂ΔT34/∂y and D.Figure 5
Curves of∂2ΔT34/∂y∂D and D.Figure6 is a surface graph of ∂ΔT34/∂y, y, and D. The blue horizontal axis in the horizontal plane is y value, unit m, simulation range from 0 m to -40 m. The green horizontal axis in the horizontal plane is D value, unit m, simulation range from 0 m to 100 m. The vertical axis is the value of ∂ΔT34/∂y, unit nT/m, variation range from 0 nT/m to 0.1 nT/m. It can be seen that ∂ΔT34/∂y varies with y and D. It is a three-dimensional surface graph.Figure 6
Three-dimensional curves of∂ΔT34/∂y, y, and D.Figure7 is a graph of ∂ΔT34/∂y and D with different y values. The horizontal axis is D value, unit m, value range from 0 m to 100 m. The vertical axis is the value of ∂ΔT34/∂y, unit nT/m, variation range from 0 nT/m to 0.1 nT/m. In Figure 7, the curve varies with D when the target y value is 40.00 m, 33.13 m, 26.27 m, 12.55 m, 5.69 m, -1.18 m, -14.9 m, -21.77 m, and -28.63 m.Figure 7
Curves ofD with different y values.
(a)
Curves wheny=40
(b)
Curves wheny=33.13
(c)
Curves wheny=26.27
(d)
Curves wheny=12.55
(e)
Curves wheny=5.69
(f)
Curves wheny=−1.18
(g)
Curves wheny=−14.90
(h)
Curves wheny=−21.77
(i)
Curves wheny=−28.63Therefore, the maximumD of theory changes with the target location y. When the environmental noise is 10 pT and the location accuracy is 1 m, the range of D value which meets the precision requirement is calculated according to the data in Figure 7, as shown in Table 1.Table 1
Value ofD which meets the precision requirement.
y (m)
D (m)
40.00
D∈4.5,32.5∪36,75.5
33.13
D∈3.5,25∪28.5,68.5
26.27
D∈3.5,16.5∪21,62.5
12.55
D∈2,54
5.69
D∈1,53
-1.18
D∈1,56
-14.90
D∈2.5,16.5∪21,67
-21.77
D∈2.5,25∪28.5,73.5
-28.63
D∈3.5,32.5∪36,80According to the requirements of accuracy, the intersection ofD value set in Table 1 is D∈4.5,16.5∪21,25∪28.5,32.5∪36,53m. In addition to satisfying the requirement of location accuracy, the array scale should be minimized to make it more flexible, following the principle of Section 2.3. The minimum value Dmin=4.5m is selected as the optimal aperture of CS-L array.


#### 3.2.2. Results with Fluxgate Magnetometers

We also performed the simulation with the array of HS-MS-FG3S fluxgate magnetometers made in China. The frequency-domain noise of this magnetometer is 10 pT/Hz1/2@1 Hz per channel. When we measure the total geomagnetic field by it, the total noise of three channels is 17.32 pT/Hz1/2@1 Hz, and the total noise is 20 pT/Hz1/2@1 Hz after superimposing with the white noise of geomagnetic field. When the positioning accuracy is 1 m, the range of D value satisfying the accuracy requirement is calculated according to the data in Figure 7, as shown in Table 2.Table 2
Value ofD which meets the precision requirement.
y (m)
D (m)
40.00
D∈9.5,32∪38.5,67
33.13
D∈7,24.5∪31.5,60.5
26.27
D∈6.5,16∪24,54
12.55
D∈4,43
5.69
D∈2,40.5
-1.18
D∈2,41.5
-14.90
D∈5,12.5∪21,52
-21.77
D∈4,22∪29,58.5
-28.63
D∈6.5,29.5∪36,65The intersection ofD value set in Table 2 which meets the requirements of accuracy is D∈9.5,12.5∪38.5,40.5m. In addition to satisfying the requirement of location accuracy, the array scale should be minimized to make it more flexible, following the principle of Section 2.3. The minimum value Dmin=9.5m is selected as the optimal aperture of HS-MS-FG3S array. It can be seen that the increase of instrument noise has an effect on the array aperture, which makes the optimal aperture of the array increase accordingly.


## 3.1. Simulation of the Experimental Environment

As shown in Figure2, four optically pumped magnetometers T1,T2,T3,T4 form a square array.Figure 2
Schematic diagram and photo of the experimental area.
(a)
(b)The aperture of the array isD. The target car moves at a constant velocity from y=32.8m to y=−41m along a parallel line 32.02 m from the y axis. Local magnetic dip is 63.3°, and local magnetic declination is 10.34°.The curve of the magnetic field changed over time in the experimental area, as shown in Figure3. The horizontal axis represents time t, and the unit is s. The vertical axis is the double gradient value of magnetic field, and the unit is nT. The array aperture is D=6m. The red line in Figure 3 is the measurement value T1 of sensor 1 (for the convenience of making the drawing, each value minus 55315 nT). T1 includes the geomagnetic field which varies with time and the magnetic field generated by the target car which varies by car’s movement. Because the geomagnetic field itself varies greatly over time, the red line cannot be used to judge when a car passed. The blue line is the data of ΔT14, which is processed by the double gradient of formula (6). The time-varying and uneven spatial distribution of the geomagnetic field is filtered using the double gradient of the magnetic field. The sensitivity of sensor CS-L is 0.6 pT/Hz1/2@1 Hz (rms). The peak value of the instrument noise is 2 pT, and the bandwidth is 0.1 Hz. The measurement noise generated by the double gradient of local environmental magnetic field is 10 pT. The measurement value of the car magnetic moment is Pm=405Am2, α=0.592rad, and β=3.74rad.Figure 3
Curve of magnetic field variation over time in experimental area.


## 3.2. Results of the Simulation

### 3.2.1. Results with the Optically Pumped Magnetometers

Figure4 shows the curve of the partial derivative ∂ΔT34/∂y and D in simulation that y=−40m. The horizontal axis is D value, unit m, value range (0 m, 100 m). The vertical axis is the value of ∂ΔT34/∂y, unit nT/m, variation range (0 nT/m, 0.1 nT/m). When D=32m, the curve has a maximum value. Thus, the theoretical optimum aperture D is 32 m, and the location accuracy of the target is the highest. The noise of environment and instrument is 10 pT. If the location accuracy is 1 m, it is required that ∂ΔT34/∂y≥10pT/m in the curve of Figure 4. In fact, D∈9,43.5∪47,91m satisfies the location accuracy, so Dmin=9m is chosen to minimize the array scale. Figure 5 is the curve of the second mixed partial derivative ∂2ΔT34/∂y∂D and D in simulation that y=−40m. The horizontal axis is D value, unit m, value range (0 m, 100 m). The vertical axis is the value of ∂2ΔT34/∂y∂D, unit nT/m2, variation range (0 nT/m2, 0.1 nT/m2). There are two D values corresponding to the second derivative ∂2ΔT34/∂y∂D=0, 32 m, and 59 m. There is a maximum value for the first derivative ∂ΔT34/∂y at D=32m.Figure 4
Curves of∂ΔT34/∂y and D.Figure 5
Curves of∂2ΔT34/∂y∂D and D.Figure6 is a surface graph of ∂ΔT34/∂y, y, and D. The blue horizontal axis in the horizontal plane is y value, unit m, simulation range from 0 m to -40 m. The green horizontal axis in the horizontal plane is D value, unit m, simulation range from 0 m to 100 m. The vertical axis is the value of ∂ΔT34/∂y, unit nT/m, variation range from 0 nT/m to 0.1 nT/m. It can be seen that ∂ΔT34/∂y varies with y and D. It is a three-dimensional surface graph.Figure 6
Three-dimensional curves of∂ΔT34/∂y, y, and D.Figure7 is a graph of ∂ΔT34/∂y and D with different y values. The horizontal axis is D value, unit m, value range from 0 m to 100 m. The vertical axis is the value of ∂ΔT34/∂y, unit nT/m, variation range from 0 nT/m to 0.1 nT/m. In Figure 7, the curve varies with D when the target y value is 40.00 m, 33.13 m, 26.27 m, 12.55 m, 5.69 m, -1.18 m, -14.9 m, -21.77 m, and -28.63 m.Figure 7
Curves ofD with different y values.
(a)
Curves wheny=40
(b)
Curves wheny=33.13
(c)
Curves wheny=26.27
(d)
Curves wheny=12.55
(e)
Curves wheny=5.69
(f)
Curves wheny=−1.18
(g)
Curves wheny=−14.90
(h)
Curves wheny=−21.77
(i)
Curves wheny=−28.63Therefore, the maximumD of theory changes with the target location y. When the environmental noise is 10 pT and the location accuracy is 1 m, the range of D value which meets the precision requirement is calculated according to the data in Figure 7, as shown in Table 1.Table 1
Value ofD which meets the precision requirement.
y (m)
D (m)
40.00
D∈4.5,32.5∪36,75.5
33.13
D∈3.5,25∪28.5,68.5
26.27
D∈3.5,16.5∪21,62.5
12.55
D∈2,54
5.69
D∈1,53
-1.18
D∈1,56
-14.90
D∈2.5,16.5∪21,67
-21.77
D∈2.5,25∪28.5,73.5
-28.63
D∈3.5,32.5∪36,80According to the requirements of accuracy, the intersection ofD value set in Table 1 is D∈4.5,16.5∪21,25∪28.5,32.5∪36,53m. In addition to satisfying the requirement of location accuracy, the array scale should be minimized to make it more flexible, following the principle of Section 2.3. The minimum value Dmin=4.5m is selected as the optimal aperture of CS-L array.


### 3.2.2. Results with Fluxgate Magnetometers

We also performed the simulation with the array of HS-MS-FG3S fluxgate magnetometers made in China. The frequency-domain noise of this magnetometer is 10 pT/Hz1/2@1 Hz per channel. When we measure the total geomagnetic field by it, the total noise of three channels is 17.32 pT/Hz1/2@1 Hz, and the total noise is 20 pT/Hz1/2@1 Hz after superimposing with the white noise of geomagnetic field. When the positioning accuracy is 1 m, the range of D value satisfying the accuracy requirement is calculated according to the data in Figure 7, as shown in Table 2.Table 2
Value ofD which meets the precision requirement.
y (m)
D (m)
40.00
D∈9.5,32∪38.5,67
33.13
D∈7,24.5∪31.5,60.5
26.27
D∈6.5,16∪24,54
12.55
D∈4,43
5.69
D∈2,40.5
-1.18
D∈2,41.5
-14.90
D∈5,12.5∪21,52
-21.77
D∈4,22∪29,58.5
-28.63
D∈6.5,29.5∪36,65The intersection ofD value set in Table 2 which meets the requirements of accuracy is D∈9.5,12.5∪38.5,40.5m. In addition to satisfying the requirement of location accuracy, the array scale should be minimized to make it more flexible, following the principle of Section 2.3. The minimum value Dmin=9.5m is selected as the optimal aperture of HS-MS-FG3S array. It can be seen that the increase of instrument noise has an effect on the array aperture, which makes the optimal aperture of the array increase accordingly.


## 3.2.1. Results with the Optically Pumped Magnetometers

Figure4 shows the curve of the partial derivative ∂ΔT34/∂y and D in simulation that y=−40m. The horizontal axis is D value, unit m, value range (0 m, 100 m). The vertical axis is the value of ∂ΔT34/∂y, unit nT/m, variation range (0 nT/m, 0.1 nT/m). When D=32m, the curve has a maximum value. Thus, the theoretical optimum aperture D is 32 m, and the location accuracy of the target is the highest. The noise of environment and instrument is 10 pT. If the location accuracy is 1 m, it is required that ∂ΔT34/∂y≥10pT/m in the curve of Figure 4. In fact, D∈9,43.5∪47,91m satisfies the location accuracy, so Dmin=9m is chosen to minimize the array scale. Figure 5 is the curve of the second mixed partial derivative ∂2ΔT34/∂y∂D and D in simulation that y=−40m. The horizontal axis is D value, unit m, value range (0 m, 100 m). The vertical axis is the value of ∂2ΔT34/∂y∂D, unit nT/m2, variation range (0 nT/m2, 0.1 nT/m2). There are two D values corresponding to the second derivative ∂2ΔT34/∂y∂D=0, 32 m, and 59 m. There is a maximum value for the first derivative ∂ΔT34/∂y at D=32m.Figure 4
Curves of∂ΔT34/∂y and D.Figure 5
Curves of∂2ΔT34/∂y∂D and D.Figure6 is a surface graph of ∂ΔT34/∂y, y, and D. The blue horizontal axis in the horizontal plane is y value, unit m, simulation range from 0 m to -40 m. The green horizontal axis in the horizontal plane is D value, unit m, simulation range from 0 m to 100 m. The vertical axis is the value of ∂ΔT34/∂y, unit nT/m, variation range from 0 nT/m to 0.1 nT/m. It can be seen that ∂ΔT34/∂y varies with y and D. It is a three-dimensional surface graph.Figure 6
Three-dimensional curves of∂ΔT34/∂y, y, and D.Figure7 is a graph of ∂ΔT34/∂y and D with different y values. The horizontal axis is D value, unit m, value range from 0 m to 100 m. The vertical axis is the value of ∂ΔT34/∂y, unit nT/m, variation range from 0 nT/m to 0.1 nT/m. In Figure 7, the curve varies with D when the target y value is 40.00 m, 33.13 m, 26.27 m, 12.55 m, 5.69 m, -1.18 m, -14.9 m, -21.77 m, and -28.63 m.Figure 7
Curves ofD with different y values.
(a)
Curves wheny=40
(b)
Curves wheny=33.13
(c)
Curves wheny=26.27
(d)
Curves wheny=12.55
(e)
Curves wheny=5.69
(f)
Curves wheny=−1.18
(g)
Curves wheny=−14.90
(h)
Curves wheny=−21.77
(i)
Curves wheny=−28.63Therefore, the maximumD of theory changes with the target location y. When the environmental noise is 10 pT and the location accuracy is 1 m, the range of D value which meets the precision requirement is calculated according to the data in Figure 7, as shown in Table 1.Table 1
Value ofD which meets the precision requirement.
y (m)
D (m)
40.00
D∈4.5,32.5∪36,75.5
33.13
D∈3.5,25∪28.5,68.5
26.27
D∈3.5,16.5∪21,62.5
12.55
D∈2,54
5.69
D∈1,53
-1.18
D∈1,56
-14.90
D∈2.5,16.5∪21,67
-21.77
D∈2.5,25∪28.5,73.5
-28.63
D∈3.5,32.5∪36,80According to the requirements of accuracy, the intersection ofD value set in Table 1 is D∈4.5,16.5∪21,25∪28.5,32.5∪36,53m. In addition to satisfying the requirement of location accuracy, the array scale should be minimized to make it more flexible, following the principle of Section 2.3. The minimum value Dmin=4.5m is selected as the optimal aperture of CS-L array.


## 3.2.2. Results with Fluxgate Magnetometers

We also performed the simulation with the array of HS-MS-FG3S fluxgate magnetometers made in China. The frequency-domain noise of this magnetometer is 10 pT/Hz1/2@1 Hz per channel. When we measure the total geomagnetic field by it, the total noise of three channels is 17.32 pT/Hz1/2@1 Hz, and the total noise is 20 pT/Hz1/2@1 Hz after superimposing with the white noise of geomagnetic field. When the positioning accuracy is 1 m, the range of D value satisfying the accuracy requirement is calculated according to the data in Figure 7, as shown in Table 2.Table 2
Value ofD which meets the precision requirement.
y (m)
D (m)
40.00
D∈9.5,32∪38.5,67
33.13
D∈7,24.5∪31.5,60.5
26.27
D∈6.5,16∪24,54
12.55
D∈4,43
5.69
D∈2,40.5
-1.18
D∈2,41.5
-14.90
D∈5,12.5∪21,52
-21.77
D∈4,22∪29,58.5
-28.63
D∈6.5,29.5∪36,65The intersection ofD value set in Table 2 which meets the requirements of accuracy is D∈9.5,12.5∪38.5,40.5m. In addition to satisfying the requirement of location accuracy, the array scale should be minimized to make it more flexible, following the principle of Section 2.3. The minimum value Dmin=9.5m is selected as the optimal aperture of HS-MS-FG3S array. It can be seen that the increase of instrument noise has an effect on the array aperture, which makes the optimal aperture of the array increase accordingly.


## 4. The Location Experiment and Data Analysis

We set four CS-L cesium optically pumped magnetometers in the suburb as shown in Figure2. The target moved along the planned trajectory parallel to the y axis direction in the horizontal plane, x=32.02m. The horizontal axis in Figure 8 is the values of y, y∈−41,32.8m. The vertical axis is the double gradient values of total geomagnetic field, which is with the unit nT, and t0 is the moment y=32.8m. The experiment is divided into four situations: D=1m, 2 m, 4 m, and 6 m. Red symbol ○ in Figure 8 represents the y value of double gradient ΔT14,ΔT24,ΔT34 measured by experiment. The blue line in Figure 8 is the theoretical curve of simulation with MathCAD and the target location y. Within the scope of the car moving, the curve of the simulation is consistent with the experimental data. It proves that the far-field model of magnetic dipole is correct in describing the magnetic field model of the car. In the experimental process, the time-varying and uneven spatial distribution of geomagnetic total field is eliminated, which proves that formula (5) and formula (6) are correct.Figure 8
Experimental data and theoretical curves ofΔTij.
(a)
D=1m, curves of ΔT14 and y
(b)
D=1m, curves of ΔT24 and y
(c)
D=1m, curves of ΔT34 and y
(d)
D=2m, curves of ΔT14 and y
(e)
D=2m, curves of ΔT24 and y
(f)
D=2m, curves of ΔT34 and y
(g)
D=4m, curves of ΔT14 and y
(h)
D=4m, curves of ΔT24 and y
(i)
D=4m, curves of ΔT34 and y
(j)
D=6m, curves of ΔT14 and y
(k)
D=6m, curves of ΔT24 and y
(l)
D=6m, curves of ΔT34 and yAccording to formulae (5) and (6), the position (x, y) of the target at each time is calculated using the experiment data, as shown in Figure 8. Formula (6) is a high-order nonlinear equation, which is difficult to solve. The software LINGO is used for the numerical solution, and the location results based on the magnetic field are shown in Table 3.Table 3
Results of location experiment.
No.
D=1m, x0=9.10m, Pm=406α=0.87, β=4.25
D=2m, x0=32.02m, Pm=512α=0.92, β=4.335
D=4m, x0=32.02m, Pm=476α=0.85, β=4.21
D=6m, x0=32.02m, Pm=405α=0.592, β=3.74
y0 (m)
x (m)
y (m)
∆r (m)
y0 (m)
x (m)
y (m)
∆r (m)
y0 (m)
x (m)
y (m)
∆r (m)
y0 (m)
x (m)
y (m)
∆r (m)
1
19.52
2.12
-2.22
22.82
19.42
28.08
14.21
6.53
19.42
33.49
19.92
1.55
18.76
30.61
17.22
2.09
2
12.81
2.04
-2.12
16.51
12.55
34.56
14.57
3.25
12.55
33.88
13.75
2.21
11.56
33.71
13.02
2.23
3
6.10
7.11
4.91
2.32
5.69
34.81
7.42
3.29
5.69
32.09
5.76
0.10
4.36
33.19
5.51
1.64
4
-0.61
8.48
-0.09
0.81
-1.18
19.70
-3.85
12.61
-1.18
29.61
-1.97
2.53
-2.84
23.90
-4.36
8.26
5
-7.32
2.00
-2.69
8.48
-8.04
23.52
-9.91
8.70
-8.04
31.56
-8.10
0.47
-10.04
35.51
-8.93
3.66
6
-14.03
16.37
-2.58
13.56
-14.90
24.68
-18.03
7.98
-14.90
27.86
-16.03
4.31
-17.24
35.39
-17.27
3.37
7
-20.74
2.84
-7.21
14.91
-21.77
40.00
-20.80
8.04
-21.77
39.06
-20.49
7.16
-24.44
21.01
-26.17
11.15
8
-27.45
2.00
-2.85
25.60
-28.63
20.00
7.22
37.82
-28.63
33.16
-29.01
1.20
-31.64
28.80
-36.01
5.43
9
-34.16
2.00
-2.70
32.25
-35.50
30.82
-36.77
1.75
-35.50
33.80
-32.79
3.25
-38.84
31.20
-41.00
2.31
Δr¯ (m)
15.25
9.99
2.53
4.46In Table3, (x0, y0) is the actual location of the target car, and (x, y) is the target location calculated using the geomagnetic measurements, with the unit m. Δr is the distance between the measuring point and the actual location of the target; Δr¯ is the average of the distance, with the unit m. The smaller the Δr¯, the higher the positioning accuracy. Based on the proposed method of determining the optimum aperture, the minimum array aperture Dmin=4.5m satisfied the requirement in the simulation when the location accuracy is 1 m. In the table, it can be seen that Δr¯ decreases from 15.25 m to 2.53 m when D increases from 1 m to 4 m. Then the location accuracy is getting higher, which meet with the simulation results. However, when D increases to 6 m, the location accuracy is not higher, and Δr¯=4.46m. It is because of the multivalueness of equation (5), and LINGO does not converge to the real optimal solution. In order to solve this problem, the additional conditions or criteria should be added to the algorithm. One way is to increase the number of linearly independent magnetic field sensor.


## 5. Conclusion

The method of constructing an array of total field magnetometers to locate a target is proposed. This method is based on a magnetic dipole model of far field. We filtered out the time-varying and uneven spatial distribution by the conjugate gradient algorithm. Then the method of determining an array aperture is proposed. We simulate the location algorithm and the process of determining an array aperture. The location experiment is carried out in the suburb. We choose 1 m, 2 m, 4 m, and 6 m, respectively, as the aperture. The experimental data is consistent with the theoretical curve. In the location algorithm, interactive linear and general optimization solver LINGO is used to solve the nonlinear equations. The location accuracy of four apertures is calculated, respectively, which basically agrees with the simulation. In the experiment, the CS-L cesium optically pumped magnetometers are used. The CS-L cesium optically pumped magnetometers have a high resolution. Also, the measurement is not affected by the temperature. The attitude of the magnetometers is not rigidly calibrated. In the location experiment, the diagonal of the car is 4.64 m and the magnetic moment is about405~512Am2. The vertical distance between the array and the car is 32.02 m, and the parallel distance is -41 m to 32.8 m. The experimental location results are shown in Figure 9. When the array aperture is 4 m, the average deviation of the magnetic moment in locating the car is 2.53 m. In a word, the method of locating the target and determining the aperture is feasible.Figure 9
Results of location experiment.
(a)
WhenD=1, results of location experiment
(b)
WhenD=2, results of location experiment
(c)
WhenD=4, results of location experiment
(d)
WhenD=6, results of location experiment



---

*Source: 7491871-2019-03-07.xml*