Position Correction and Trajectory Optimization of Underwater Long-Distance Navigation Inspired by Sea Turtle Migration

Abstract

Accumulating evidence suggests that migrating animals store navigational “maps” in their brains, decoding location information from geomagnetic information based on their perception of the magnetic field. Inspired by this phenomenon, a novel geomagnetic inversion navigation framework was proposed to address the error constraint of a long-distance inertial navigation system. In the first part of the framework, the current paper proposed a geomagnetic bi-coordinate inversion localization approach which enables an autonomous underwater vehicle (AUV) to estimate its current position from geomagnetic information like migrating animals. This paper suggests that the combination of geomagnetic total intensity (F) and geomagnetic inclination (I) can determine a unique geographical location, and that there is a non-unique mapping relationship between the geomagnetic parameters and the geographical coordination (longitude and latitude). Then the cumulative error of the inertial navigation system is corrected, according to the roughly estimated position information. In the second part of the framework, a cantilever beam model is proposed to realize the optimal correction of the INS historical trajectory. Finally, the correctness of the geomagnetic bi-coordinate inversion localization model we proposed was verified by outdoor physical experiments. In addition, we also completed a geomagnetic/inertial navigation integrated long-distance semi-physical test based on the real navigation information of the AUV. The results show that the geomagnetic inversion navigation framework proposed in this paper can constrain long-distance inertial navigation errors and improve the navigation accuracy by 73.28% compared with the pure inertial navigation mode. This implies that the geomagnetic inversion localization will play a key role in long-distance AUV navigation correction.

1. Introduction

Underwater navigation technology is fundamental to reaching deep and pelagic resource development [1,2]. Especially for long-distance navigation, the navigation and positioning functions have a higher priority [3]. Due to its autonomy, the inertial navigation system (INS) has been an object of research for AUVs since the last century. For long-distance navigation, inertial navigation is an indispensable core device. So far, however, the rapid accumulation of errors over time is still problematic in INS [4,5,6]. Periodic errors, offsets, accelerometer noise, and other interferences make it difficult to maintain high navigation accuracy for unassisted inertial navigation [7]. At present, the most advanced portable inertial navigation system merely has a positioning accuracy of 1.11 km/h, and its price is approximately 200,000 dollars [8].

Nowadays, the most effective solution for high-precision long-distance navigation is inertial navigation, and other navigation technologies assist in correcting or constraining inertial navigation errors. For example, underwater acoustic technology and INS play crucial roles in limiting the drift of INS over time [9]. The combination of Long Base Line (LBL)\Short Base Line (SBL)\INS has been extensively studied, but it is still a challenging project to preset acoustic beacons within a range of thousands of kilometers. In addition, Doppler Velocity Log (DVL) and INS integrated navigation will also be affected by DVL installation errors and speed blind spots, resulting in a decrease of accuracy [10].

The geomagnetic field is an inherent physical field on the earth. It can be used as a stable navigation information source because of its vector and stable characteristics. Compared with underwater acoustic and inertial navigation, geomagnetic navigation has strong concealment and is not influenced by cumulative errors. Several papers relevant to this research were published, including matching and positioning techniques based on a geomagnetic map [11,12,13,14], and a target search approach for geomagnetic parameters (TSAGP) [15,16,17,18,19]. On the basis of matching algorithms, a geomagnetic matching algorithm was presented [11,20]. However, the matching technology had some difficulties. Firstly, the premise of the algorithms was that the vehicle needs to accumulate a certain length of geomagnetic information, which affected the real-time performance of geomagnetic matching navigation. Secondly, the fault-tolerant adaptation of this algorithm’s capacity is expected to be improved. Thirdly, there are few magnetic anomaly areas in the marine environment, resulting in insufficient magnetic features.

Liu et al. [16,21] carried out many studies on the TSAGP. Inspired by biological navigation behavior, the magnetic searching approach was proposed to solve geomagnetic navigation without an a priori map, and a related improved method was presented to deal with magnetic interference. Although the TSAGP does not require a priori geomagnetic map, it is challenging for multiple parameters to converge at the target point simultaneously [22].

In addition to humans, animals can also intelligently use magnetic fields to complete field navigation and positioning [23,24,25]. For example, sea turtles utilize geomagnetic navigation without the actual geomagnetic map during migration. It is not known how turtles navigate and find their homing on a voyage of thousands of kilometers. For decades, biologists and marine researchers have been studying this exciting phenomenon. Animals’ geomagnetic navigation includes a “geomagnetic map” and a “geomagnetic compass” [26]. The geomagnetic compass can provide a heading direction that animals can then use to orient in any desired direction, and the geomagnetic map tells animals their local coordinates. The most compelling explanation for the magnetic map hypothesis is based on magnetic nanoparticles of magnetite, which are sensitive to changes in the strength and inclination of the earth’s magnetic field. When the strength or direction of the magnetic field changes, the ion channels in the cell will be open or closed under the interaction between the magnetite nanoparticles and the geomagnetic field, forming chemical signals that the higher central nervous system can decode. Then, animals encode navigation information in their brains. The encoding process of geomagnetic information can be understood as inputting the geomagnetic field data to the intelligent processing center. At the same time, decoding is the use of a well-established relationship model to convert the input data into output location information. This decoding model is called geomagnetic imprinting, which informs the relationship between geomagnetic environments and geographical location.

Inspiringly, we find that the I and F contours can extend in regular directions in most regions of the Earth’s surface. It is worth noting that the intersection of two contours is almost unique, which means that the geomagnetic parameters at the intersection can uniquely represent a position on the Earth.

This paper proposes a geomagnetic bi-coordinate inversion localization (GBCIL) approach inspired by sea turtles’ homing phenomenon. We assume a nonlinear relationship between the geomagnetic bi-coordinates and the geographic coordinates, as Figure 1. Through this relationship, we will obtain artificial geomagnetic imprinting. Its function is to decode the measured geomagnetic field information into geographic location. Furthermore, we can use the calculated position information to make preliminary corrections and historical trajectory optimization for arbitrarily diverging inertial navigation.

The contribution of this paper is proposing a novel geomagnetic inversion navigation framework. This framework includes two parts, one is a geomagnetic bi-coordinate inversion localization (GBCIL) approach to realize the rough correction of the long-distance INS errors. The other is the introduction of a cantilever beam model to optimize the historical trajectory of the inertial navigation. The rest of this paper is arranged as follows: Section 2 demonstrates models, including the geomagnetic model, the GBCIL model, and the cantilever beam model. In Section 3, the physical experiment and a semi-physical simulation are carried out to verify our proposed approach. Finally, we draw some conclusions in Section 4.

2. Model

2.1. Earth Magnetic Field Geomagnetic Model

In this section, the characteristics of the geomagnetic field are described. The Earth is like a giant magnet, with directions that can be used to locate and discern a direction [27]. There are many scientific explanations for the source of the geomagnetic field. At present, the most widely accepted theory is that the magnetic field is produced by the outer core of the earth [28]. Thus, any position on Earth corresponds to a different geomagnetic vector. Seven parameters can describe the geomagnetic field vectors (see Figure 2). They are total intensity F, northerly intensity X, easterly intensity Y, vertical intensity Z, horizontal intensity H, declination angle D, and inclination I, respectively. The declination is the angle between the vectors H and X. The angle between vector F and the horizontal plane is inclination, and the I downward is positive [29].

Many geomagnetic models are used to approximately describe the Earth’s magnetic field, such as International Geomagnetic Reference Field (IGRF), High Definition Geomagnetic Model, and the World Magnetic Model (WMM). The World Magnetic Model uses a spherical harmonic representation to degree and order 12 to represent the main geomagnetic field. However, the main geomagnetic field only accounts for slightly more than 95% of the field strength at the Earth’s surface. The crustal and combined disturbance fields also play significant roles in the geomagnetic fields. In contrast, The IGRF is regarded as a standard mathematical description of the Earth’s main magnetic field by The International Association of Geomagnetism and Aeronomy (IAGA) because IGRF data combinates the magnetic field model and data from satellites, observatories and surveys around the world. Since there is less magnetic interference and fewer anomalies in the marine environment, the IGRF database can more accurately describe the true field strength of the geomagnetic field [30].

2.2. Geomagnetic Bi-Coordinate Inversion Localization Model

The main magnetic field is a potential field and can be written as the negative spatial gradient of a scalar potential. This potential can be expanded in terms of spherical harmonics using Equation (1).

$$V\left(\lambda ,{\phi }^{’},r,t\right)=a{\displaystyle \sum _{n=1}^{12}{\left(\frac{a}{r}\right)}^{n+1}{\displaystyle \sum _{m=0}^n\left(g_n^m\left(t\right)\cos \left(m\lambda \right)+h_n^m\left(t\right)\sin \left(m\lambda \right)\right){\breve{P}}_n^m\left(\sin {\phi }^{’}\right)}}$$

(1)

where, V is scalar potential of the main magnetic field, a donates the geomagnetic reference radius, $g_n^m\left(t\right)$ and $h_n^m\left(t\right)$ are the gauss coefficients of degree n and order m describing the Earth’s main magnetic field, and $\mu and {\breve{P}}_n^m$ are the Schmidt semi-normalized associated Legendre functions. The other geomagnetic vector components can be calculated based on Equation (1), and more detailed descriptions can refer to [27].

We can obtain geomagnetic data under the condition that longitude $\lambda$, latitude ${\phi }^{’}$, radius r, and time t are known. However, it is difficult for AUV to obtain accurate geographical information (longitude and latitude) in the underwater environment. In addition, the actual geomagnetic data can be measured by a magnetometer. For this reason, we attempted to consider this issue in reverse as Figure 3. The aim of Equation (1) is to obtain the geomagnetic parameters under known geographic location. The proposed GBCIL is used to calculate the geographic location under known geomagnetic parameters. They are opposite in terms of the known variables and the variables to be solved. We call this inversion mapping.

The proposed approach in this paper is a reversal solution compared with WMM. Firstly, the F and I of the planned route area are extracted from the IGRF model or the actual measurement database. A Taylor expansion was applied to transform the model into polynomial form for a linear solution for F and I. Complex numbers are introduced to express the longitude and latitude coordinates. We suppose that $\lambda$ and $\varphi$ denotes the longitude and latitude, respectively, and i is the imaginary unit. Thus, the combination of longitude with latitude is written $\lambda +\varphi \mathrm{i}$ to ensure the combination of latitude and longitude uniqueness. The model is written as Equation (2).

$$\lambda +\varphi \mathrm{i}={\displaystyle \sum _{i=0}^n{\displaystyle \sum _{j=0}^i{A}_{i-j,j}{\left(F-F_0\right)}^{i-j}{\left(I-I_0\right)}^j}}$$

(2)

where F₀ and I₀ donate expand point of variable F and I, respectively, and $A_{i-j,j}$ is polynomial coefficient. The right of the Equation (2) is written in detail,

$${\displaystyle \sum _{i=0}^n{\displaystyle \sum _{j=0}^i{A}_{i-j,j}{\left(F-F_0\right)}^{i-j}{\left(I-I_0\right)}^j}=}{A}_{0,0}{\left(F-F_0\right)}^0{\left(I-I_0\right)}^0+A_{1,0}{\left(F-F_0\right)}^{1-0}{\left(I-I_0\right)}^0+\dots +A_{n,n}{\left(F-F_0\right)}^n{\left(I-I_0\right)}^n$$

(3)

Then Equation (3) is rewritten

$${\displaystyle \sum _{i=0}^n{\displaystyle \sum _{j=0}^i{A}_{i-j,j}{\hat{F}}^{i-j}{\hat{I}}^j}=}\left[\begin{array}{cccc}{A}_{0,0},& A_{1,0},& \dots & ,A_{n,n}\end{array}\right]\left[\begin{array}{c}{\hat{F}}^0{\hat{I}}^0\\ {\hat{F}}^1{\hat{I}}^0\\ \dots \\ {\hat{F}}^n{\hat{I}}^n\end{array}\right]=AU$$

(4)

where $\hat{F}=F-F_0$, $\hat{I}=I-I_0$, $A=\left[\begin{array}{cccc}{A}_{0,0},& A_{1,0},& \dots & ,A_{n,n}\end{array}\right]$, $U={[{\hat{F}}^0{\hat{I}}^0,{\hat{F}}^1{\hat{I}}^0,\dots ,{\hat{F}}^n{\hat{I}}^n]}^{\mathrm{T}}$, the left of Equation (2) is written.

$$Y=\lambda +\varphi \mathrm{i}$$

Therefore, Equation (2) is simplified as

$$Y=AU$$

(5)

Then, the original problem is transformed into solving coefficient matrix $A$ under the condition of known Y and $U$. Equation (5) is solved by the Least Square Method (LSM), which is expressed as

$$\underset{{\rm X}\in R^{\mathrm{n}}}{\min }f\left(A\right)={\Vert AU-Y\Vert }_2$$

(6)

Making $\frac{\partial f\left(A\right)}{\partial U}=0$. Then, Equation (6) becomes

$$\frac{\partial f\left(A\right)}{\partial U}=2U^{T}UA-2U^{T}Y=0$$

(7)

The coefficient matrix $A$ link geomagnetic and geographic is

$$A={\left(U^{T}U\right)}^{-1}{U}^{T}Y$$

(8)

Therefore, the relationship between the geographic coordinates λ, $\varphi$ and the geomagnetic bi-coordinates F, I is determined in the selected area.

2.3. Trajectory Correction Model Based on Cantilever Beam

The AUV real path points are $X=\left\{x_0,x_1,\dots ,x_n\right\}$, and the INS path points are $\hat{X}=\left\{{\hat{x}}_0,{\hat{x}}_1,\dots ,{\hat{x}}_n\right\}$, where ${\hat{x}}_0=x_0$. The INS model can be described as

$$\begin{array}{l}{\hat{X}}_n={\hat{X}}_{n-1}+V_{n-1}\mathrm{\Delta }t+\frac{{\hat{a}}_{n-1}\mathrm{\Delta }{t}^2}{2}\\ ={\hat{X}}_{n-2}+V_{n-2}\mathrm{\Delta }t+\frac{{\hat{a}}_{n-2}\mathrm{\Delta }{t}^2}{2}+V_{n-1}\mathrm{\Delta }t+\frac{{\hat{a}}_{n-1}\mathrm{\Delta }{t}^2}{2}\\ ={\hat{X}}_{n-2}+\left(V_{n-2}+V_{n-1}\right)\mathrm{\Delta }t+\left({\hat{a}}_{n-2}-{\hat{a}}_{n-1}\right)\frac{\mathrm{\Delta }{t}^2}{2}\\ ={\hat{X}}_0+\left(V_0+\dots +V_{n-2}+V_{n-1}\right)\mathrm{\Delta }t+\left({\hat{a}}_0+\dots +{\hat{a}}_{n-2}-{\hat{a}}_{n-1}\right)\frac{\mathrm{\Delta }{t}^2}{2}\end{array}$$

(9)

Owing to $V_n=V_{n-1}+{\hat{a}}_{n-1}\mathrm{\Delta }t=V_0+\left({\hat{a}}_0+\dots +{\hat{a}}_{n-2}+{\hat{a}}_{n-1}\right)\mathrm{\Delta }t$, Equation (9) can be further simplified to

$${\hat{X}}_n={\hat{X}}_0+\left[n{V}_0+\left(n-1\right){\hat{a}}_0+\left(n-2\right){\hat{a}}_1+\dots +{\hat{a}}_{n-1}\right]\mathrm{\Delta }t+\left({\hat{a}}_0+\dots +{\hat{a}}_{n-2}+{\hat{a}}_{n-1}\right)\frac{\mathrm{\Delta }{t}^2}{2}$$

(10)

where ${\hat{a}}_0$ donates the acceleration measured by accelerometer and $V_0$ represents the initial speed, set to zero.

$${\hat{a}}_i=a_i+\underset{{\epsilon }_i}{\underbrace{{\mu }_i+b_i}}$$

(11)

where $a_i$ is the real acceleration of AUV, and $\mu$ and $b$ represent the random noise and bias of the accelerometer, respectively. Considering Equations (9)–(11), we obtain the following relationship.

$$\begin{array}{l}{\hat{X}}_n={\hat{X}}_0+\left[\left(n-1\right){\hat{a}}_0+\left(n-2\right){\hat{a}}_1+\dots +{\hat{a}}_{n-1}\right]\mathrm{\Delta }t+\left({\hat{a}}_0+\dots +{\hat{a}}_{n-2}+{\hat{a}}_{n-1}\right)\frac{\mathrm{\Delta }{t}^2}{2}\\ ={\hat{X}}_0+\left[\left(n-1\right)a_0+\left(n-2\right)a_1+\dots +a_{n-1}\right]\mathrm{\Delta }t+\left[\left(n-1\right){\epsilon }_0+\left(n-2\right){\epsilon }_0+\dots +{\epsilon }_{n-1}\right]\mathrm{\Delta }t+\left({\displaystyle \sum _{i=0}^{n-1}{a}_i}+{\displaystyle \sum _{i=0}^{n-1}{\epsilon }_i}\right)\frac{\mathrm{\Delta }{t}^2}{2}\\ ={\hat{X}}_0+\left[\left(n-1\right)a_0+\left(n-2\right)a_1+\dots +a_{n-1}\right]\mathrm{\Delta }t+\frac{\mathrm{\Delta }{t}^2}{2}{\displaystyle \sum _{i=0}^{n-1}{a}_i}+\left[\left(n-1\right){\epsilon }_0+\left(n-2\right){\epsilon }_0+\dots +{\epsilon }_{n-1}\right]\mathrm{\Delta }t+\frac{\mathrm{\Delta }{t}^2}{2}{\displaystyle \sum _{i=0}^{n-1}{\epsilon }_i}\\ =X_0+K+\ell \end{array}$$

(12)

where $K=\left[\left(n-1\right)a_0+\left(n-2\right)a_1+\dots +a_{n-1}\right]\mathrm{\Delta }t+\frac{\mathrm{\Delta }{t}^2}{2}{\displaystyle \sum _{i=0}^{n-1}{a}_i}$ represents the input of the AUV motion model, and $\ell =\left[\left(n-1\right){\epsilon }_0+\left(n-2\right){\epsilon }_0+\dots +{\epsilon }_{n-1}\right]\mathrm{\Delta }t+\frac{\mathrm{\Delta }{t}^2}{2}{\displaystyle \sum _{i=0}^{n-1}{\epsilon }_i}$ represents the cumulative error of the inertial navigation. Therefore,

$$\{\begin{array}{l}{X}_n=X_0+K\\ \ell ={\hat{X}}_n-X_n\end{array}$$

(13)

This paper proposes an error correction method for the cantilever beam model. The cantilever beam is a structure in structural mechanics, as shown in Figure 4. Assuming that the beam is lightweight and does not account for gravity, one end of the beam is fixed to the wall, and then an external force is applied to the beam, causing the beam to bend and deform, which is called deflection [31]. Setting the end fixed on the wall as the origin of the coordinates, the deflection of the beam at x from the origin can be described as

$$v=\frac{q{x}^2}{24EI}\left(x^2-4lx+6l^2\right)$$

(14)

where $v$ donates deflection, $l$ donates the length of the beam, $q$ represents the load acting on the beam, E is the tensile modulus of elasticity and I is the moment of inertia. According to the bending of the beam under the load, the mechanical model is converted into a motion model. The beam will maintain its original “real” shape in the mechanical model without external force. Similarly, the inertial navigation system will record the “real” trajectory of the AUV without interference from the error term. In addition, future errors of INS will all be affected by the previous error. Therefore, the error at each discrete time in the motion model can be described as

$$v_i=\frac{{\epsilon }_i{i}^2}{24}\left(i^2-4li+6l^2\right)$$

(15)

where i represents the i-th path point in the track that needs to be corrected, $v_i$ represents the corrected error values at the i-th path point, and l represents the total number of track points that need to be corrected. The modulus of elasticity and the moment of inertia are the coefficients related to the rod material, which are treated as 1 in the motion model, and the role it plays in the model is included in $v_i$. ${\epsilon }_i$ is the error of INS at each discrete time. Supposing that, at the kth time, the current position of the AUV is obtained through the GBCIL model as $X_k^{\mathrm{M}}={\left[x_k^{\mathrm{M}},y_k^{\mathrm{M}}\right]}^T$, the INS position is $\hat{X}={\left[\hat{x},\hat{y}\right]}^T$. So, from departure to time k, the cumulative error of the INS can be approximated as

$$v_k=X_k^M-{\hat{X}}_k$$

(16)

$v_k$ is a two-dimensional vector. Substituting Equation (16) into Equation (15), we can obtain

$${\epsilon }_k=\frac{8v_k}{k^4}$$

(17)

In this model, considering that the error of the INS at each moment is the same, which are equal to the error at time k, so we have

$${\epsilon }_0={\epsilon }_1=\dots ={\epsilon }_{k-1}={\epsilon }_k=\frac{8v_k}{k^4}$$

(18)

Finally, the error correction values at each discrete time can be obtained by followings

$$v_i=\frac{2v_k{i}^2}{k^2}-\frac{4v_k{i}^3}{3k^3}+\frac{v_k{i}^4}{3k^4}$$

(19)

3. Experiment and Results

Long-distance geomagnetic surveying, especially over thousands of kilometers, is a large project. In this article, in order to verify the feasibility of our proposed geomagnetic bi-coordinates to geographic bi-coordinate mapping, the GBCIL model was tested on a small scale under limited conditions. In addition, we used AUV long-distance INS data and IGRF to complete large-scale semi-physical simulations.

3.1. Physical Experiment

Aiming to verify the GBCIL model we proposed, we built a surface geomagnetic navigation experimental platform, as shown in Figure 5d. The platform includes a GPS/inertial integrated navigation system, a 3D fluxgate transmitter, a microcomputer processing module, and a mobile power supply. The integrated navigation system and the fluxgate sensor are rigidly connected to obtain the attitude of the fluxgate, as shown in Figure 5. The test was carried out in the Liquan Lake Reservoir in Xianyang, which has an average water surface width of 200 m. The one-way sailing distance of the experimental platform was about 4 km, and the round trip totaled 8 km. The sampling frequency of the fluxgate sensor was 10 HZ, and the measurement trajectory is shown in Figure 6.

In our GBCIL approach, the relationship between geographical coordination and geomagnetic data is a built-in navigation frame. However, the geomagnetic information was measured by a 3D fluxgate transmitter in the body frame. Therefore, it is necessary to employ a three-axis gyroscope to measure the angle of the body frame relative to navigation frame of the magnetic field digital transmitter. $M_b$ denotes the measured data in body frame, and $M_n$ denotes the geomagnetic data in navigation frame. A direction cosine matrix $C_b^n$ was used to transform $M_b$ to $M_n$. The notations $C_b^n$ can be written as follows:

$$C_b^n=\left(\begin{array}{ccc}\cos \psi \cos \theta & \sin \psi \cos \theta & -\sin \theta \\ \cos \psi \sin \theta \sin \gamma -\sin \psi \cos \gamma & \sin \psi \sin \theta \sin \gamma +\cos \psi \cos \gamma & \cos \theta \sin \gamma \\ \cos \psi \sin \theta \cos \gamma +\sin \psi \sin \gamma & \sin \psi \sin \theta \cos \gamma -\cos \psi \sin \gamma & \cos \theta \cos \gamma \end{array}\right)$$

(20)

where $\psi$, $\theta$, $\gamma$ denote the yaw, pitch, and roll of the 3D fluxgate transmitter, respectively. A calibration method motioned in [32] of the tri-axial magnetometer was carried out to compensate sensor output.

In the calibration procedure, we installed the fluxgate on the calibration device and then rotated around the X, Y, and Z axes of the fluxgate to measure the magnetic field data at more angles. The results are shown in Figure 7, in which the red stars represent the original measurement data, and the green circles represents the corrected data. Then, we used the obtained compensation parameters to compensate the data measured by the navigation of the experimental platform. The real measurement results are shown in Figure 8 and Figure 9.

IGRF provides a global magnetic field reference value with a resolution of $0.01°$. The magnetic inclination angle of our test result was $53°{30}^{\prime }$, and the total intensity was 52.75–53.01 μT (Figure 8 and Figure 9). Comparing our measurement results, we find that our measurement values were close to the reference values provided by IGRF (F = 52.98 μT, $I=53°{30}^{\prime }$ at 108.43 E, 34.53 N). Then, we input the measured geomagnetic field parameters and real position into the GBCIL model to solve the inversion coefficients further. We selected 24 samples from the measured data as a validation set to test the inversion accuracy of our proposed model in the test area. Except for the extracted data, the rest were used as the training set.

Figure 10 and Figure 11 show the results of our proposed GBCIL on the water surface. The black line is the actual trajectory of the float measured by GPS, and the black rectangle represents the measurement position of the geomagnetic data. The green circle represents the inversion position of the proposed approach. It can be seen from the figure that the localization result obtained by GBCIL was exceptionally close to GPS. In Figure 11, the blue bar represents the absolute error between GPS position and the inversion position. It is evident from this figure that most of the absolute errors were less than 200 m, and the most significant error was 389.8 m. The standard deviation of the proposed approach was 104.1715 m. The simulation experimental and sea trial results illustrate that the GBCIL approach could effectively obtain the localization information by measuring geomagnetic data.

3.2. Semi-Physical Simulation

In order to test the ability of our proposed geomagnetic bi-coordinate inversion localization approach to correct long-distance INS errors, this section will conduct a long-distance INS/geomagnetic fusion navigation experiment. Measuring a virtual geomagnetic library spanning thousands of kilometers is a huge challenge, and our team currently does not have this detection capability. Therefore, we extracted the geomagnetic data at the location of the AUV from the IGRF database as the actual measured values. Inertial navigation data and real navigation data are measured in actual sea trials. The AUV used in the experiment was an Orange Shark-III, 6000 mm in length, 680 mm in diameter, and 135 kg in weight, as Figure 12 shows. The Orange Shark-III was equipped with optical fiber inertial navigation (gyroscope and accelerometer bias stability of 0.01°/s and 50 μg, respectively). Under pure inertial navigation conditions, the error was 2.16 km/h.

The proposed GBCIL is more advantageous for long-distance navigation, owing that it does not produce error accumulation. First, we calculated the positioning accuracy of the GBCIL model in the AUV navigation area through simulation experiments. The relationship between the geomagnetic bi-coordinate grid and the geographic grid with size is 0.01° was obtained. Furthermore, we designed different navigation paths in two other areas to test the positioning accuracy of GBCIL. The first area was at 109.5–115 E, 16.5–22 N; we developed a lawnmower path in this area. The second area was 123–130 E, 18–23 N; we created a circular path in this area. In the simulation, according to the magnetic sensor accuracy, gaussian white noise with variances of 0.5 and 0.005 was added to the total intensity F and magnetic inclination I data as the measurement error of the sensor. The results are shown in Figure 13, Figure 14, Figure 15 and Figure 16.

It can be seen from Figure 14 and Figure 16 that the positioning accuracy of our proposed GBCIL model does not have cumulative error propagation, with an average absolute error of 714.68 m and a standard deviation of 495.93 m in the lawnmower path, and with an average absolute error 483.40 m and a standard deviation of 293.96 m in the circular path. Combining Figure 10, Figure 13 and Figure 15, it can be seen that our proposed GBCIL model can complete localization in multiple regions. It should be pointed out that the inversion accuracy is related to the grid size of the magnetic field data. When the geomagnetic grid accuracy is higher, the inversion accuracy of GBCIL will also be improved. Considering that the INS error of Orange Shark III is 2.16 km/h, we corrected the INS position every hour by utilization of the GBCIL model in the long-distance semi-physical navigation experiment, and then used the cantilever beam proposed model to optimize and compensate for the historical trajectory of INS. The results are shown in Figure 17 and Figure 18.

As shown in Figure 19, Orange Shark-III navigated about 14 h, and the INS error has accumulated to 25 km. After rough correction of GBCIL and historical trajectory optimization, the maximum error decreased to 6.7 km, and the standard deviation was 1.658 km. In Figure 15, the trough of the solid black line represents the correction point, and the INS error has apparent convergence after each correction. It can be seen from Figure 14 and Figure 15 that the corrected trajectory is very close to the GPS trajectory of the AUV. The maximum error was 1.67% of the full range, while the original inertial navigation error reached 6.25%. In contrast, our proposed method improved the navigation accuracy by 73.28%.

4. Conclusions

Inspired by animal long-distance migration navigation, this paper proposes a new geomagnetic inversion navigation framework to solve the problem of INS cumulative error during long-distance navigation. First, we demonstrated a geomagnetic bi-coordinate inversion localization approach to correct the INS error for long-distance navigation. Then, we introduced a cantilever beam to optimize the INS historical trajectory to obtain more credible navigation information. Physical experiments and long-distance semi-physical simulation verified this approach. The small-scale physical experiment proved the viability of our proposed GBCIL model. The long-distance semi-physical simulation tested the proposed method’s performance in the error correction of long-distance navigation. The results show that the proposed method can effectively converge the INS error, and that the navigation accuracy of INS improved by 73.28%.

It can be seen from the experimental results that the geomagnetic inversion navigation framework we proposed has excellent practicability. The current answer “where am I?” can be quickly obtained for measurement information. The approach proposed in this paper addresses localization information in the global navigation system. In future work, we will consider fusing other navigation information into the geomagnetic inversion navigation framework to improve the system’s robustness.