2. GEOMAGNETIC GRADIENT TENSOR AND ITS MEASUREMENT

The geomagnetic gradient tensor G is the spatial variation rate of the geomagnetic vector along three directions in a Cartesian coordinate system, namely

(1)$$G = \left[ {\matrix{ {\displaystyle{{\partial H_x} \over {\partial x}}} & {\displaystyle{{\partial H_x} \over {\partial y}}} & {\displaystyle{{\partial H_x} \over {\partial z}}} \cr {\displaystyle{{\partial H_y} \over {\partial x}}} & {\displaystyle{{\partial H_y} \over {\partial y}}} & {\displaystyle{{\partial H_y} \over {\partial z}}} \cr {\displaystyle{{\partial H_z} \over {\partial x}}} & {\displaystyle{{\partial H_z} \over {\partial y}}} & {\displaystyle{{\partial H_z} \over {\partial z}}} \cr}} \right] = \left[ {\matrix{ {g_{xx}} & {g_{xy}} & {g_{xz}} \cr {g_{yx}} & {g_{yy}} & {g_{yz}} \cr {g_{zx}} & {g_{zy}} & {g_{zz}} \cr}} \right]$$

where, H _x, H_y and H _z are the three components in x, y and z directions, respectively. We can obtain g _xx + g _yy + g _zz = 0, g _xy = g _yx, g_xz = g _zx, g_yz = g _zy because the divergence and rotation of the geomagnetic field are all zero. Therefore, the magnetic gradient tensor is a 3 × 3 symmetrical matrix whose trace is zero and there exist five independent components selected as g _xx, g_yy, g_xy, g_yz and g _zx respectively.

To measure the geomagnetic gradient tensor accurately, a type of sensor configuration scheme called a ten single-axis configuration is adopted and shown in Figure 1, where the arrows represent the sensitive axes of sensors. Two three-axis magnetometers are located at position A to measure the magnetic field denoted as (h ₁, h ₂, h ₃) and B to measure the magnetic field denoted as (h ₄, h ₅, h ₆), respectively. Similarly, the magnetic fields (h ₇, h ₈) and (h ₉, h ₁₀) are measured by two two-axis magnetometers located respectively at positions C and D. The baseline lengths of gradient measurement in the directions x and y are respectively l _x and l _y. Then, the geomagnetic gradient can be given by

(2)$$\left\{ \matrix{g_{xx} = \displaystyle{{h_4 - h_1} \over {l_x}}, g_{yy} = \displaystyle{{h_9 - h_7} \over {l_y}}, g_{yx} = \displaystyle{{h_5 - h_2} \over {l_x}} \hfill \cr g_{zy} = \displaystyle{{h_{10} - h_8} \over {l_y}}, g_{zx} = \displaystyle{{h_6 - h_3} \over {l_x}} \hfill} \right.$$

Figure 1. Diagram of geomagnetic gradient tensor measurement configuration.

The local perturbations or anomalies from rocks, terrain, man-made objects, etc can be viewed as the components of Earth's magnetic gradient tensor, and these tensor measurements can be accurately surveyed by a precise gradient tensor system through calibration, in which resolution can reach 0·52nT/m (Sui et al., Reference Sui, Li, Wang and Lin2014).

3. ATTITUDE MEASUREMENT PRINCIPLE BASED ON MAGNETIC GRADIENT TENSOR

As shown in Figure 2, the navigation coordinate system is represented by ox _ny_nz_n, whose three axes are respectively east, north and up, and ox _by_bz_b is the vehicle coordinate system. The three angles ψ, θ and γ are the heading, pitch and roll angles. Consequently, the relationship between the navigation coordinate and vehicle coordinate can be described as follows: firstly, ox ₁y ₁z ₁ is obtained by rotating ox _ny_nz_n with ψ around −oz _n; secondly, Ox ₂y ₂z ₂ is obtained by rotating Ox ₁y ₁z ₁ with θ around ox ₁; finally, Ox _by_bz_b can be given by rotating Ox ₂y ₂z ₂ around oy ₂.

Figure 2. Schematic diagram for attitude measurement principle.

The transformation matrix C _n^b from the navigation coordinate system to the vehicle coordinate system can be written as

(3)$$C_n^b = \left[ {\matrix{ {\cos \gamma \cos \psi + \sin \gamma \sin \psi \sin \theta} & { - \cos \gamma \sin \psi + \sin \gamma \cos \psi \sin \theta} & { - \sin \gamma \cos \theta} \cr {\sin \psi \cos \theta} & {\cos \psi \cos \theta} & {\sin \theta} \cr {\sin \gamma \cos \psi - \cos \gamma \sin \psi \sin \theta} & { - \sin \gamma \sin \psi - \cos \gamma \cos \psi \sin \theta} & {\cos \gamma \cos \theta} \cr}} \right]$$

The quantities Gⁿ and G^b show the magnetic gradient tensors in the navigation coordinate system and vehicle coordinate system, respectively. Thus,

(4)$$G^b = C_n^b G^n (C_n^b )^T $$

where, T expresses the matrix transpose. The quantity G^b can be obtained by a magnetic gradient measurement device because it is generally strapped on the vehicle. On the other hand, Gⁿ can be measured in advance and stored in the navigation computer as the navigation geophysical database. From Equation (4), we can obtain five independent equations for the attitude angles (ψ, θ, γ). To avoid trigonometric function operation of the optimal estimation, the quaternion Q = q ₀ + q ₁i ₀ + q ₂j ₀ + q ₃k ₀ is introduced to determine the transform matrix, namely

(5)$$C_b^n = (C_b^n )^T = \left[ {\matrix{ {q_0^2 + q_1^2 - q_2^2 - q_3^2} & {2(q_1 q_2 - q_0 q_3 )} & {2(q_1 q_3 + q_0 q_2 )} \cr {2(q_1 q_2 + q_0 q_3 )} & {q_0^2 - q_1^2 + q_2^2 - q_3^2} & {2(q_2 q_3 - q_0 q_1 )} \cr {2(q_1 q_3 - q_0 q_2 )} & {2(q_2 q_3 + q_0 q_1 )} & {q_0^2 - q_1^2 - q_2^2 + q_3^2} \cr}} \right]$$

From Equations (4) and (5), we can deduce the following five independent equations with regard to q ₀, q ₁, q ₂ and q ₃

(6)$$\eqalign{ g_{xx}^b & = g_{xx}^n {\left( {q_0^2 + q_1^2 - q_2^2 - q_3^2 } \right)^2} + 4g_{yy}^n {\left( {{q_1}{q_2} + {q_0}{q_3}} \right)^2} \cr &\quad + 4g_{xy}^n \left( {q_0^2 + q_1^2 - q_2^2 - q_3^2 } \right)\left( {{q_1}{q_2} + {q_0}{q_3}} \right) \cr & \quad+ 8g_{zy}^n \left( {{q_1}{q_2} + {q_0}{q_3}} \right)\left( {{q_1}{q_3} - {q_0}{q_2}} \right) \cr & \quad+ 4g_{xz}^n \left( {q_0^2 + q_1^2 - q_2^2 - q_3^2 } \right)\left( {{q_1}{q_3} - {q_0}{q_2}} \right) \cr &\quad - 4\left( {g_{xx}^n + g_{yy}^n } \right){\left( {{q_1}{q_3} - {q_0}{q_2}} \right)^2}} $$

(7)$$\eqalign{ g_{yy}^b & = 4g_{xx}^n {\left( {{q_1}{q_2} - {q_0}{q_3}} \right)^2} + g_{yy}^n {\left( {q_0^2 - q_1^2 + q_2^2 - q_3^2 } \right)^2} \cr &\quad+ 4g_{yx}^n \left( {q_0^2 - q_1^2 + q_2^2 - q_3^2 } \right)\left( {{q_1}{q_2} - {q_0}{q_3}} \right) \cr & \quad+ 4g_{zy}^n \left( {q_0^2 - q_1^2 + q_2^2 - q_3^2 } \right)\left( {{q_2}{q_3} + {q_0}{q_1}} \right) \cr &\quad + 8g_{zx}^n \left( {{q_1}{q_2} - {q_0}{q_3}} \right)\left( {{q_2}{q_3} + {q_0}{q_1}} \right) \cr & \quad- 4\left( {g_{xx}^n + g_{yy}^n } \right){\left( {{q_2}{q_3} - {q_0}{q_1}} \right)^2}} $$

(8)$$\eqalign{& g_{yx}^b = 2g_{xx}^n \left( {q_0^2 + q_1^2 - q_2^2 - q_3^2} \right)\left( {q_1 q_2 - q_0 q_3} \right){\kern 1pt} {\kern 1pt} {\kern 1pt} + 2g_{yy}^n \left( {q_0^2 - q_1^2 + q_2^2 - q_3^2} \right)\left( {q_1 q_2 + q_0 q_3} \right) \cr & \quad{\kern 1pt} {\kern 1pt} {\kern 1pt} \matrix{ {} & {} \cr} + g_{yx}^n \left[ {\left( {q_0^2 + q_1^2 - q_2^2 - q_3^2} \right)\left( {q_0^2 - q_1^2 + q_2^2 - q_3^2} \right) + 4\left( {q_1 q_2 - q_0 q_3} \right)\left( {q_1 q_2 + q_0 q_3} \right)} \right] \cr & {\kern 1pt} {\kern 1pt} \matrix{ {} & {} \cr} \quad{\kern 1pt} + 2g_{zy}^n \left[ {2\left( {q_1 q_2 + q_0 q_3} \right)\left( {q_2 q_3 + q_0 q_1} \right) + \left( {q_0^2 - q_1^2 + q_2^2 - q_3^2} \right)\left( {q_1 q_3 - q_0 q_2} \right)} \right] \cr & \quad{\kern 1pt} {\kern 1pt} {\kern 1pt} \matrix{ {} & {} \cr} + 2g_{zx}^n \left[ {2\left( {q_0^2 + q_1^2 - q_2^2 - q_3^2} \right)\left( {q_2 q_3 + q_0 q_1} \right) + 2\left( {q_1 q_3 - q_0 q_2} \right)\left( {q_1 q_2 - q_0 q_3} \right)} \right]{\kern 1pt} {\kern 1pt} {\kern 1pt} \cr & \quad\matrix{ {} & {} \cr} - 4\left( {g_{xx}^n + g_{yy}^n} \right)\left( {q_1 q_3 - q_0 q_2} \right)\left( {q_2 q_3 + q_0 q_1} \right)} $$

(9)$$\eqalign{& g_{zy}^b = 4g_{xx}^n \left( {q_1 q_2 - q_0 q_3} \right)\left( {q_1 q_3 + q_0 q_2} \right){\kern 1pt} {\kern 1pt} {\kern 1pt} + 2g_{yy}^n \left( {q_0^2 - q_1^2 + q_2^2 - q_3^2} \right)\left( {q_2 q_3 + q_0 q_1} \right) \cr & \quad{\kern 1pt} {\kern 1pt} {\kern 1pt} \matrix{ {} & {} \cr} + 2g_{yx}^n \left[ {2\left( {q_1 q_2 - q_0 q_3} \right)\left( {q_2 q_3 - q_0 q_1} \right) + \left( {q_0^2 - q_1^2 + q_2^2 - q_3^2} \right)\left( {q_1 q_3 + q_0 q_2} \right)} \right] \cr & \quad{\kern 1pt} {\kern 1pt} {\kern 1pt} \matrix{ {} & {} \cr} + g_{zy}^n \left[ {\left( {q_0^2 - q_1^2 + q_2^2 - q_3^2} \right)\left( {q_0^2 - q_1^2 - q_2^2 + q_3^2} \right) + 4\left( {q_0 q_1 + q_2 q_3} \right)\left( {q_2 q_3 - q_0 q_1} \right)} \right] \cr & \quad{\kern 1pt} {\kern 1pt} {\kern 1pt} \matrix{ {} & {} \cr} + 2g_{zx}^n \left[ {2\left( {q_0^2 - q_1^2 - q_2^2 + q_3^2} \right)\left( {q_1 q_2 - q_0 q_3} \right) + 2\left( {q_2 q_3 + q_0 q_1} \right)\left( {q_1 q_3 - q_0 q_2} \right)} \right] \cr & \quad{\kern 1pt} {\kern 1pt} {\kern 1pt} \matrix{ {} & {} \cr} - 2\left( {g_{xx}^n + g_{yy}^n} \right)\left( {q_0^2 - q_1^2 - q_2^2 + q_3^2} \right)\left( {q_2 q_3 + q_0 q_1} \right)} $$

(10)$$\eqalign{& g_{zx}^b = 2g_{xx}^n \left( {q_0^2 + q_1^2 - q_2^2 - q_3^2} \right)\left( {q_1 q_3 + q_0 q_2} \right){\kern 1pt} {\kern 1pt} + 4g_{yy}^n \left( {q_1 q_2 + q_0 q_3} \right)\left( {q_2 q_3 - q_0 q_1} \right) \cr & \quad{\kern 1pt} {\kern 1pt} {\kern 1pt} \matrix{ {} & {} \cr} + 2g_{yx}^n \left[ {2\left( {q_1 q_2 + q_0 q_3} \right)\left( {q_1 q_3 + q_0 q_2} \right) + \left( {q_0^2 + q_1^2 - q_2^2 - q_3^2} \right)\left( {q_2 q_3 + q_0 q_1} \right)} \right] \cr & \quad{\kern 1pt} {\kern 1pt} {\kern 1pt} \matrix{ {} & {} \cr} + 2g_{zy}^n \left[ {\left( {q_0^2 - q_1^2 - q_2^2 + q_3^2} \right)\left( {q_1 q_2 + q_0 q_3} \right) + 2\left( {q_1 q_3 - q_0 q_2} \right)\left( {q_2 q_3 - q_0 q_1} \right)} \right] \cr & \quad{\kern 1pt} {\kern 1pt} {\kern 1pt} \matrix{ {} & {} \cr} + g_{zx}^n \left[ {2\left( {q_0^2 + q_1^2 - q_2^2 - q_3^2} \right)\left( {q_0^2 - q_1^2 - q_2^2 + q_3^2} \right) + 4\left( {q_1 q_3 + q_0 q_2} \right)\left( {q_1 q_3 - q_0 q_2} \right)} \right] \cr & {\kern 1pt} {\kern 1pt} {\kern 1pt} \matrix{ {} & {} \cr}\quad - 2\left( {g_{xx}^n + g_{yy}^n} \right)\left( {q_0^2 - q_1^2 - q_2^2 + q_3^2} \right)\left( {q_1 q_3 - q_0 q_2} \right)} $$

In addition there still exists the following relation for the quaternion

(11)$$q_0^2 + q_1^2 + q_2^2 + q_3^2 = 1$$

The above six equations can be replaced by

(12)$$f_i (q_0, q_1, q_2, q_3 ) = 0\,(i = 1,2,\ldots,6)$$

Equation (12) is a group of nonlinear equations, which is traditionally solved using an iterative method, the Newton method, the steepest descent method or conjugate gradient method, etc. According to characteristics of the nonlinear equations shown as Equation (12), the Newton Downhill method is chosen, which is a fast and effective method for solving nonlinear equations. The iterative formula of Newton Downhill method is

(13)$$x^{n + 1} = x^n - \omega \left( {F'\left( {x^n} \right)} \right)^{ - 1} F\left( {x^n} \right)$$

where 0 < ω ≤ 1. To guarantee convergence, ω generally meets

(14)$$\left\Vert {F\left( {x^{n + 1}} \right)} \right\Vert \lt \left\Vert {F\left( {x^n} \right)} \right\Vert$$

We adopt the successive halving method to determine ω. To lighten the burden of calculation, partial derivatives in the Jacobian matrix are replaced by difference approximation.

The solutions $\hat q_0 $, $\hat q_1 $, $\hat q_2 $, $\hat q_3 $ of Equation (12) are given using the above iterative algorithm, and then matrix elements $\hat c_{ij} $ of C _bⁿ are obtained (i, j = 1, 2, 3) through substituting $\hat q_0 $,$\hat q_1 $,$\hat q_2 $,$\hat q_3 $ into Equation (5). The following relations can be given using Equations (3) and (4),

(15)$$\left\{ \matrix{\hat \gamma _{main} = \arctan \left( { - \displaystyle{{\hat c_{31}} \over {\hat c_{33}}}} \right) \hfill \cr \theta = \arcsin \left( {\hat c_{32}} \right) \hfill \cr \hat \psi _{main} = \arctan \left( {\displaystyle{{\hat c_{12}} \over {\hat c_{22}}}} \right) \hfill} \right.$$

where, $\hat \theta $ is the estimated pitch angle; $\hat \psi _{main} $ and $\hat \gamma _{main} $ are respectively the main values of estimated heading and roll angles. The estimated heading and roll angles are decided by the rules shown in Tables 1 and 2, respectively.

Table 1. The truth table of heading angle.

Table 2. The truth table of roll angle.

The detailed steps for underwater full attitude determination based on geomagnetic gradient tensor measurement are described as follows.

Step 1: Calculate the five components of geomagnetic gradient tensor in the vehicle coordinate system according to the sensor output from the magnetic tensor measurement device and Equation (2);

Step 2: Extract the five components of magnetic gradient tensor in the navigation coordinate system from the geomagnetic navigation database according to the position indicated by the inertial navigation system;

Step 3: Solve the nonlinear Equation (12) and obtain $\hat q_0 $, $\hat q_1 $, $\hat q_2 $ and $\hat q_3 $ using the Newton Downhill method;

Step 4: Substitute $\hat q_0 $, $\hat q_1 $, $\hat q_2 $ and $\hat q_3 $ into Equation (5) to calculate $\hat c_{ij} $, and then calculate three attitude angles according to Equation (15) and Tables 1 and 2.

4. NUMERICAL EXPERIMENTS

We use a single magnetic dipole to produce a geomagnetic anomaly gradient, whose magnetic moments in three directions of the navigation coordinate system are m _x = 10⁹ A·m², m _y = 2 × 10⁸ A·m² and m _z = 10⁸ A·m², respectively. The spatial positions of the magnetic dipole relative to the magnetic gradient tensor measurement device are x = 100 m, y = 50 m and z = 20 m in the navigation coordinate system, respectively. The baseline lengths of magnetic gradient measurement are both 1 m, namely l _x = l _y = 1 m. The iterative termination condition of the Newton Downhill algorithm is selected as ε = ||x _k+1− x _k|| ≤ 10⁻⁶.

Firstly, the influence of the initial solution on the algorithm accuracy is investigated. The noise of each magnetometer is assumed to be an independent Gaussian noise process with zero mean and variance of σ = 2nT.

Figure 3 shows the convergence of the Newton Downhill method for different initial solutions, where δΦ represents the difference of initial solutions deviating from the true solution. When the difference is very obvious, for example, δΦ = 20^o and δΦ = 30^o, the error e(t) will fluctuate markedly. But when δΦ = 20^o, the Newton Downhill method still shows convergence after a few iterations, which means that the error is less than the given minimum.

Figure 3. The convergence algorithm under different initial solution conditions.

The Newton Downhill method generally shows good convergence through not more than 20 iteration steps when δΦ ≤ 20^o, which is proved by many simulations. This shows the full attitude determination algorithm costs less computing time and possesses potential application value.

Figure 4 gives the estimated attitude angles in the case of different initial angle errors. From Figure 4, we can see the full attitude determination algorithm is not sensitive to initial angle error, and can estimate effectively the attitude angles as long as initial angle error is less than 20°, where the estimation errors of θ, ψ and γ are about 0·47°, 0·15°, 0·12°, respectively. This further indicates that the full attitude determination algorithm possesses potential application value. The output attitude error of the inertial navigation system is less than 10°, which can be used as the initial attitude of the attitude determination algorithm. The error of the full attitude determination algorithm based on a geomagnetic gradient tensor is less than 0·5°, which remarkably improves the precision of the navigation system.

Figure 4. The estimated attitude angles for different initial solutions.

Next, the influence of the sensor noise on the algorithm accuracy is investigated with other invariable simulation parameters. The resolution and absolute precision of a good fluxgate magnetometer such as the RM100 are respectively 0·1 nT and 1 nT. The resolution of a superconducting vector magnetometer such as the 755R can reach 2 × 10⁻³ nT. The resolution and relative accuracy of a laser optical pump helium vector magnetic sensor such as the P-2000 are better than 0·01 nT and 0·1 %, respectively. Figure 5 shows the convergence of the algorithm for different sensor noise expressed by σ, which is lower than 10 nT.

Figure 5. Algorithm convergence under different noise levels.

From Figure 5, we can see that convergence of the algorithm is not obviously influenced by the sensor noise as long as it is less than 10 nT, when the maximum number of iterations is 18. However, the convergence speed slows significantly when the noise is higher than 10 nT.

Figure 6 gives the estimated attitude angles in conditions of different magnetometer noise. When the noise is less than 10 nT, the estimation errors of θ, ψ and γ are about 0·35°, 0·11°, 0·13°, respectively. The noise of a low-cost three-axis fluxgate sensor is less than 10 nT, and therefore it is feasible for the algorithm to use low-cost three-axis fluxgate sensors to construct a geomagnetic gradient tensor measurement device.

Figure 6. Attitude angle estimation under different noise levels.

During surveying, the device straps down the inertial platform with gyro of 0·01^o/h integrated with other sensors, and then its horizontal misalignment is better than 1′. The errors in platform attitude are viewed directly as attitude determination error, so its effect on attitude determination can be ignored. On the other hand, the errors in location during surveying will also lead to gradient tensor error, and the relationship between them in the surveying plane can be expressed as

(16)$$\Delta g_{ij} = \displaystyle{{\partial g_{ij}} \over {\partial x}}\Delta x + \displaystyle{{\partial g_{ij}} \over {\partial y}}\Delta y$$

where, i, j = x, y, z, Δg _ij is the gradient tensor error and Δx and Δy are the location errors. According the above simulations, the attitude determination precision is acceptable when the noise of fluxgate sensor is less than 10 nT. Correspondingly, the errors of gradient tensor are about 10 nT/m. So, Δg _ij should be less than 10 nT/m and the errors in location are determined by Equation (16).