Next Article in Journal
Affimer-Based Europium Chelates Allow Sensitive Optical Biosensing in a Range of Human Disease Biomarkers
Previous Article in Journal
Distributed Sensing Network Enabled by High-Scattering MgO-Doped Optical Fibers for 3D Temperature Monitoring of Thermal Ablation in Liver Phantom
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

An Aided Navigation Method Based on Strapdown Gravity Gradiometer

Department of Navigation Engineering, Naval University of Engineering, Wuhan 430000, China
*
Author to whom correspondence should be addressed.
Sensors 2021, 21(3), 829; https://doi.org/10.3390/s21030829
Submission received: 8 December 2020 / Revised: 12 January 2021 / Accepted: 13 January 2021 / Published: 27 January 2021
(This article belongs to the Section Remote Sensors)

Abstract

:
The gravity gradient is the second derivative of gravity potential. A gravity gradiometer can measure the small change of gravity at two points, which contains more abundant navigation and positioning information than gravity. In order to solve the problem of passive autonomous, long-voyage, and high-precision navigation and positioning of submarines, an aided navigation method based on strapdown gravity gradiometer is proposed. The unscented Kalman filter framework is used to realize the fusion of inertial navigation and gravity gradient information. The performance of aided navigation is analyzed and evaluated from six aspects: long voyage, measurement update period, measurement noise, database noise, initial error, and inertial navigation system device level. When the parameters are set according to the benchmark parameters and after about 10 h of simulation, the results show that the attitude error, velocity error, and position error of the gravity gradiometer aided navigation system are less than 1 arcmin, 0.1 m/s, and 33 m, respectively.

1. Introduction

In view of the cost of gravity gradiometers and the existence of outdoor global satellite navigation systems, more attention is paid to the special application scenarios of gravity matching under the condition of satellite rejection, especially in the military application field. Submarines adopt an integrated navigation scheme based on inertial navigation and assisted by other navigation means. When global navigation satellite systems, Loran-C, and celestial navigation are used to correct an inertial navigation system (INS), the submarine needs to be in periscope navigation state, and the concealment of the submarine cannot be guaranteed. With the aid of sonar and Doppler velocity log, the velocity and altitude information of the INS can be corrected, but the acoustic signal is broadcasted, so it is easy to expose the submarine’s position to the enemy. Through building a base station, submarine positioning in a certain sea area can be realized, but the process of base station construction and maintenance is complicated, and it does not have the ability of positioning in the whole sea area. An underwater emergency positioning buoy is also a kind of submarine aided navigation means, but this method has high economic cost and is usually only used in emergency situations. So, it has been a difficult problem in the field of navigation on how to realize submarine passive autonomous, long-voyage and high-precision navigation and positioning [1,2,3,4].
The aided navigation method based on the geophysical field provides a new idea and means to solve this problem [5,6,7]. The marine gravity field and geomagnetic field are two kinds of geophysical fields that we are mainly concerned about. As the geomagnetic field belongs to weak field signals, it is easy to be interfered with by other factors, especially by the submarine shell. The geomagnetic measurement is usually carried out by a dragging mode, which provides a poor positioning accuracy of several kilometers [8,9,10] of the submarine navigation system based on geomagnetic field signals. In addition, the geomagnetic real-time measurement is more troublesome. Compared with the geomagnetic field, the gravity field signal has stronger anti-interference and stability, and is more suitable for underwater, providing aided navigation information for inertial navigation systems [11,12,13,14]. The gravity gradient is the second derivative of the gravity potential. Relative to gravity, a gravity gradiometer can measure subtle changes of gravity signals in different geographical locations [15,16]. Therefore, gravity gradient signals contain more abundant navigation information.
The gravity gradient database, gravity gradiometer, and gravity gradient matching algorithm are the three key parts affecting the gravity information aided navigation system, and most of the research work is based on the above three factors. Over the past few decades, gravity models of terrestrial planets, especially the Earth, have been improved dramatically. Derived from the combination of satellite geodetic data with high-resolution gravitational information collected from surface gravimetry, the Earth Gravitational Model 2008 (EGM2008) is complete to degree 2190 and order 2159, which is a representative of the high-resolution models. The high order gravity field model based on EGM2008 is an effective method to establish the global gravity information database. In the future, the gravity gradiometer with atomic interference technology can theoretically realize higher precision gravity gradient measurement [17]. In previous studies, the extended Kalman filter framework was used to realize gravity gradient aided navigation (GGAN) [18,19,20]. In this framework, the mathematical model needs to be linearized. Some particle filter algorithms have also been considered in the field of gravity gradient matching algorithm [21,22]. In addition, most of the previous literatures studied the mathematical model and performance analysis based on the platform gravity gradiometer, and some factors affecting performance were considered [23,24]. However, it is necessary to investigate the performance of applying the strapdown gravity gradiometer and systematically evaluate the performance of GGAN under the influence of various factors, and there is still a need to verify the effects of some other factors, such as long voyage and inertial navigation system device level.
In this study, we focus on an unscented quaternion estimator, which has a better estimation accuracy, confirmed by previous studies [25], and establish a framework of unscented Kalman filter (UKF) based on strapdown gravity gradiometer to realize GGAN. Simulation experiments were conducted for systematically evaluating the performance of GGAN under the influence of various factors, such as long voyage, measurement update period, measurement noise, database noise, initial errors (attitude, velocity, and position), and inertial navigation system device level. As a result, we drew some numerical conclusions, which can provide some suggestions for obtaining navigation results with some target accuracies in GGAN on the basis of the simulation results.
The organization of this paper is as follows. In the review of existing theories, the basic strapdown inertial navigation system (SINS) equations and the principle of gravity gradiometer are introduced. In the gravity gradient matching algorithm section, UKF is studied in detail. Finally, simulation testing is conducted and is summarized in the experimental section. The results are provided and the performance of the six factors are evaluated. A few meaningful conclusions are outlined in the summary.

2. Basic Equations of Strapdown Inertial Navigation System

The basic equations of strapdown inertial navigation system (SINS) include attitude differential equation, velocity differential equation, position differential equation, gyroscope error differential equation, and accelerometer error differential equation. Compared with the platform inertial navigation system, SINS uses the “mathematical platform” simulation navigation coordinate system. In the basic equation of SINS, “mathematical platform” is determined by the attitude transfer matrix C b n and the attitude differential equation. Different attitude expression methods correspond to different attitude differential equations. Since quaternion attitude expression is continuous and has no singularity problem, this paper adopts quaternion representation, and its differential equation is generally expressed as follows:
q ˙ = 1 2 Ξ ( q ) ω n b b = 1 2 Ω ( ω n b b ) q ,
where q = q 0 ; ρ is the quaternion, where ρ = q 1 ; q 2 ; q 3 is the vector part, and q 0 is the scalar part, and
Ω ( ω n b b ) = [ ω n b b × ] ω n b b ( ω n b b ) T 0
Other parts of the basic SINS equations are
L ˙ = v N R M + h ,
λ ˙ = v E ( R N + h ) cos L ,
h ˙ = v U ,
v ˙ E = v E ( R N + h ) cos λ + 2 ω i e e v N sin λ v E v U R N + h 2 ω i e e v U cos λ + f E ,
v ˙ N = v E ( R N + h ) cos λ + 2 ω i e e v E sin λ v E v U R M + h + f N ,
v ˙ U = v E 2 R N + h + v N 2 R M + h + 2 ω i e e v E cos λ g f U
Equations (3)–(5) are the position differential equation. p = L ; λ ; h is the position information, where L is the latitude, λ is the longitude, and h is the height; Equations (6)–(8) are the velocity differential equations. v = v E ; v N ; v U represents velocity information, where v E represents eastward velocity, v N represents northward velocity, and v U represents vertical velocity; g is the acceleration of gravity; f = [ f E ; f N ; f U ] is the specific force in the navigation coordinate system. The output of the gyroscope is ω i b b ; R M and R N represent the meridian and prime vertical radius of the Earth, as shown in the following formula:
R M = R e ( 1 e 2 ) / ( 1 e 2 sin 2 λ ) 1.5 ,
R N = R e / ( 1 e 2 sin 2 λ ) 0.5 ,
where R e = 6378137m, e = 0.0818.
The differential equation of gyro constant drift is as follows:
ω ˜ i b b = ω i b b + ε + η g v ,
ε ˙ = η g u ,
where ω ˜ i b b is the actual gyro output with gyro drift ε ; η g v and η g u are zero mean Gaussian white noise with covariance σ g v 2 and σ g u 2 , respectively.
The output equation of the accelerometer reads as follows:
f ˜ b = f b + + η a v ,
˙ = η a u ,
where f ˜ b is the actual accelerometer output with accelerometer bias , η a v and η a u are zero mean Gaussian white noise with covariance σ a v 2 and σ a u 2 , respectively.
Equations (3)–(14) together constitute the basic equation of SINS. By analyzing the characteristics of the equation, we can find that the basic equation of SINS is a nonlinear equation, if the filtering framework is used for information fusion, the integrated navigation needs to adopt nonlinear filtering processing.

3. Measurement Principle of Gravity Gradiometer

The gravity gradiometer is basically a differential accelerometer. The measurement principle is described in detail below.
Suppose that there exists an arbitrary coordinate system a, which rotates around the inertial coordinate system with angular velocity ω i a a . If the position vector of a certain point is r a and the position vector of the point in the inertial coordinate system is r i ; r i can be calculated from
r i = C a i r a ,
where C a i is the transformation matrix between coordinate system a and inertial coordinate system i . It can be obtained by time derivation of Equation (15)
r ˙ i = C a i Ω i a a r a + C a i r ˙ a ,
if ω i a a = [ ω 1 ω 2 ω 3 ] T , the antisymmetric matrix Ω i a a can be expressed as
Ω i a a = 0 ω 3 ω 2 ω 3 0 ω 1 ω 2 ω 1 0
The time derivative of Equation (16) leads to
r ¨ i = C a i r ¨ a + 2 C a i Ω i a a r ˙ a + C a i ( Ω i a a Ω i a a + Ω ˙ i a a ) r a
In inertial space the following relation holds
r ¨ i = a i + g i ,
where a i is the specific force output of the accelerometer and g i is the gravity of the Earth. Combining Equations (18) and (19) yields
a a = r ¨ a + 2 Ω i a a r ˙ a + ( Ω i a a Ω i a a + Ω ˙ i a a ) r a g a
The two accelerometers are fixed at point A and point B in the coordinate system a , respectively, and the baseline of the two accelerometers is set as ρ a = r 2 a r 1 a . According to Equation (20), we can get
a 2 a a 1 a = ρ ¨ a + 2 Ω i a a ρ ˙ a + ( Ω i a a Ω i a a + Ω ˙ i a a ) ρ a ( g 2 a g 1 a )
Since the accelerometer is attached to coordinate system a , then ρ ˙ a = ρ ¨ a = 0 . Both quantities can be eliminated from Equation (21)
a 2 a a 1 a = ( Ω i a a Ω i a a + Ω ˙ i a a ) ρ a ( g 2 a g 1 a )
The longer the baseline is, the more obvious the change of gravity signal is. Limited to the size of the gravity gradiometer, the length of baseline usually does not exceed 1 m. Therefore, it can be considered that the gravity between the two accelerometers varies linearly. The following relation can be obtained
g 2 a g 1 a = V a ( r 2 a r 1 a ) = V a ρ a ,
where V a is the total tensor matrix of the gravity gradient in coordinate system a . By substituting Equation (23) into Equation (22), we obtain
a 2 a a 1 a = ( V a + Ω i a a Ω i a a + Ω ˙ i a a ) ρ a = L a ρ a
and
L a = ( a 2 a a 1 a ) / ρ a = ( a 2 a a 1 a ) / ( r 2 a r 1 a )
Equation (25) is the basic equation of gravity gradient based on the measurement principle of differential accelerometer. It can be seen from Equations (24) and (25) that the gravity gradient tensor cannot be directly measured, and the gradients are coupled with angular velocity Ω i a a and angular acceleration Ω ˙ i a a . In order to get the current gravity gradient, the angular velocity and angular acceleration should be eliminated. The angular acceleration component can be eliminated by summation of Equation (25) and its transposition
1 2 ( L a + ( L a ) T ) = V a + Ω i a a Ω i a a
Let
L a = V a + Ω i a a Ω i a a
By substituting each component into Equation (24), we end up with
L a = ( V x x + ω y 2 + ω z 2 ) ( V x y ω x ω y ) ( V x z ω x ω z ) ( V x y ω x ω y ) ( V y y + ω x 2 + ω z 2 ) ( V y z ω y ω z ) ( V x z ω x ω z ) ( V y z ω y ω z ) ( V z z + ω x 2 + ω y 2 )
The six independent components in expression (28) are as follows:
L a = L 11 a L 22 a L 33 a L 12 a L 13 a L 23 a = ( V x x + ω y 2 + ω z 2 ) ( V y y + ω x 2 + ω z 2 ) ( V z z + ω x 2 + ω y 2 ) ( V x y ω x ω y ) ( V x z ω x ω z ) ( V y z ω y ω z )
Compared with Equation (25), Equation (27) does not need to estimate Ω ˙ i a a , which is conducive to the extraction of gravity gradient tensor. So, Equation (27) is taken as the main measurement result of gravity gradiometer.
When the gradiometer is fixed to the carrier (strapdown mounting), the body coordinate system and acceleration coordinate system are the same coordinate system, so Equation (27) can be rewritten into
L b = C n b V n ( C n b ) T Ω i b b Ω i b b ,
where C n b is the transformation matrix between the navigation coordinate system and the body coordinate system, it can be obtained by attitude calculation, and Ω i b b can be retrieved from the gyro measurements. In this paper, Equation (30) is used as the measurement equation for simulation, and the calculation formula of V n is shown in the appendix.

4. Gravity Gradiometer Aided Navigation Method (GGAN Method)

The UKF algorithm abandons the linearization process of EKF algorithm, and adopts unscented transform (UT) to avoid the error caused by linearization, reduce the complexity of the algorithm, and overcome the defects of low accuracy and poor stability of the EKF algorithm. So, it is widely used in integrated navigation [20]. In this paper, the strapdown gravity gradiometer is used to measure the total tensor gravity gradient data. Therefore, the measurement vector is
Z k = L 11 b L 22 b L 33 b L 12 b L 13 b L 23 b T ,
where the superscript b represents the projection of the gravity gradient in the body coordinate system.
The state variables are the attitude φ , velocity v , and position p of the target and the drift vector b of the 6D gyroscope and accelerometer. Therefore, the state vector is
X k = φ v p b T
At the k th moment, the system noise is Gaussian white noise w k ~ N ( 0 Q k ) , and the measurement noise is Gaussian white noise v k ~ N ( 0 R k ) . The nonlinear system of gravity gradient aided positioning can be expressed by the following formula:
X k = F ( X k 1 ) + w k 1 Z k = H ( X k ) + v k
Firstly, the root mean square S k 1 + of the error covariance matrix needs to be solved in the UKF propagation process. This step can be conducted via the Cholesky decomposition
P k 1 + = S k 1 + S k 1 + T
Calculate the sigma point from the following formula:
X k 1 + ( i ) = X ^ k 1 + + n S k 1 : , i + , i n X ^ k 1 + n S k 1 : , ( i n ) + , i > n ,
where the subscript “ : , i ” denotes the column i of the matrix. Each sigma point can be propagated through the system model:
X k ( i ) = F ( k , X k 1 + ( i ) )
After propagation, the state estimation and its error covariance are as follows:
X ^ k = 1 2 n i = 1 2 n X k ( i ) ,
P k = 1 2 n i = 1 2 n ( X k ( i ) X ^ k ) ( X k ( i ) X ^ k ) T + Q k 1
The observation update process of UKF generates new sigma points by the following formula:
P k = S k S k T X k ( i ) = X ^ k + n S k , : , i , i n X ^ k n S k , : , ( i n ) , i > n
The sigma point and mean observed innovation can be solved by the following formula:
δ Z k ( i ) = Z k H ( k , X ^ k ( i ) ) δ Z k = 1 2 n i = 1 2 n δ Z k ( i ) ,
and the covariance of the observed innovation is
C δ Z , k = 1 2 n i = 1 2 n ( δ Z k ( i ) δ Z k ) ( δ Z k ( i ) δ Z k ) T + R k
Finally, the UKF Kalman gain, state vector update, and error covariance update can be calculated as follows:
K k = [ 1 2 n i = 1 2 n ( X k ( i ) X ^ k ) ( δ Z k ( i ) δ Z k ) T ] ( C δ Z , k ) 1 ,
X ^ k + = X ^ k + K k δ Z k ,
P k + = P k K k C δ Z , k K k T
When the system noise and measurement noise are Gaussian white noise, GGAN can be realized by the standard UKF localization algorithm mentioned above.

5. Performance Analysis

Since the 21st century, with the continuous maturity and development of satellite altimetry, airborne gravimetry, and other gravity measurement technologies, the resolution and accuracy of the spherical harmonic function model of the Earth’s gravity field have been continuously improved. The most representative one is the EGM2008 ultra-high order spherical gravity field spherical harmonic model, issued by the National Geospatial-Intelligence Agency [26]. EGM2008 is a fusion of the ITG-GRACE03S model and global 5′ × 5′ gravity anomaly grid data (including land gravity survey, ocean satellite altimetry, and airborne gravity survey). In this paper, the EGM2008 model is used to simulate the gravity gradient, and the 360 order is intercepted to calculate the gravity gradient. The specific calculation formula is shown in Appendix A.
In order to evaluate the performance of gravity gradiometer aided navigation, this paper comprehensively considers the influence of long voyage, measurement update period, measurement noise, database noise, initial error, inertial device level, and other factors, in which the inertial device level mainly refers to gyro bias. In this paper, the database was calculated by the spherical harmonic model, and the calculated results were taken as the true value. To mimic a real situation, we introduced noise on the basis of the calculated values of the spherical harmonic model to simulate the actual database. The benchmark parameter settings are shown in Table 1 below.

5.1. Performance Analysis under Long Voyage Condition

In order to intuitively reflect the long voyage characteristics of inertial/gravity gradient integrated navigation, the simulation experiment designs acceleration, uniform speed, deceleration, steering, and other motion forms, and the navigation area is arbitrarily selected. The simulation step size is 0.01 s, and the total simulation time is about 6.5 h. Other parameters are set according to the benchmark parameters in Table 1. It should be noted that in the pure inertial calculation process, the altitude information is always set to zero, and the experimental results are shown in Figure 1, Figure 2 and Figure 3 and Table 2 below.
It can be seen from Figure 1, Figure 2 and Figure 3 that the error of inertial navigation system is effectively suppressed by GGAN. Both attitude and speed and position accuracy have been greatly improved, and the positioning error presents zero mean distribution. According to Table 2, after about 10 h of pure inertial calculation, the position errors in latitude and longitude directions are up to 2640 m and 2258 m, respectively. However, the position performance of GGAN method is two orders of magnitude better than that of pure inertial calculation method, and the attitude error, velocity error, and position error in all directions of GGAN system are less than 1 arcmin, 0.1 m/s, and 33 m, respectively.

5.2. Performance Analysis of Measurement Update Period

In order to evaluate the influence of different measurement update periods of gravity gradiometer on integrated navigation performance, four different measurement update periods (30, 60, 90, and 180 s) were selected for simulation experiments. Other parameters were set according to the benchmark parameters in Table 1. The simulation results are shown in the following Figure 4 and Table 3. According to Table 3, with the increase of measurement update period, the positioning error of the system gradually increases. The reason is that the observation information collected over the same time span gets less when using a larger measurement update period. When the measurement update period is less than 180 s, the position error in other directions is within 100 m, except for the longitude direction of 127 m. The 3D velocity errors for these four cases were 0.03 m/s, 0.05 m/s, 0.07 m/s, and 0.12 m/s, respectively, and the 3D position errors were 43 m, 52 m, 62 m, and 147 m.

5.3. Performance Analysis of Measurement Noise

To investigate the effects of gradiometer noise on GGAN accuracy, different levels of gradiometer noise 1 mE, 0.01 E, 0.1 E, and 1 E were simulated. The 1 E and 0.1 E noise levels represent the precision of most current generation gradiometers, such as the ARKeX’s Exploration Gravity Gradiometer and the Gedex’s High-Definition Airborne Gravity Gradiometer. The 0.01 E noise level represents the precision of the latest gradiometers, such as the GOCE’s EGG. The 0.001 E noise level represents the precision of future grade gradiometers. Other parameters were set according to the benchmark parameters in Table 1. The simulation results are shown in the following Figure 5 and Table 4. According to Table 4, with the increase of measurement noise, the positioning performance of the system decreases rapidly. When the measurement noise is 1 E, the position error in longitude direction has exceeded 1000 m. The 3D velocity errors for these four cases were 0.05 m/s, 0.05 m/s, 0.21 m/s, and 1.90 m/s, respectively, and the 3D position errors were 40 m, 45 m, 227 m, and 1533 m.

5.4. Performance Analysis of Database Noise

The accuracy of the established gravity gradient database is usually higher than that of the strapdown gradiometer. In order to reflect the impact of different database noise on the integrated navigation performance, two different database noises (0.001 and 0.01 E) were selected for our simulation experiment. Other parameters were set according to the benchmark parameters in Table 1. The simulation results are shown in the following Figure 6 and Table 5. According to Table 5, it can be seen that the gravity gradient database noise has a great impact on the performance of the GGAN system. When assuming a gradient noise of 0.001 E, the positioning performance is mostly better by one order of magnitude, compared to a noise level of 0.01 E. The 3D velocity errors for these two cases were 0.05 m/s and 1.86 m/s, respectively, and the 3D position errors were 60 m and 1509 m.

5.5. Performance Analysis of Initial Errors

In order to reflect the influence of different initial SINS errors on the GGAN performance, the error parameters were divided into attitude error, velocity error, and position error. The initial attitude error was selected as 0.1, 0.2, 0.5, and 1 degree. Other parameters were set according to the benchmark parameters in Table 1. The simulation results are shown in the following Figure 7 and Table 6. It can be seen that the simulation converges in about an hour. After about 3 h of simulation, the 3D velocity errors for these four cases were 0.06 m/s, 0.06 m/s, 0.06 m/s, and 0.6 m/s, respectively, and the 3D position errors were 51 m, 51 m, 61 m, and 577 m. Therefore, when the initial attitude error is less than 0.5 deg, the influence of attitude error is not obvious.
The initial velocity error was selected as 1, 2, 3, and 5 m/s. Other parameters were set according to the benchmark parameters in Table 1. The simulation results are shown in the following Figure 8 and Table 7. The 3D velocity errors for these four cases were 0.05 m/s, 0.05 m/s, 0.06 m/s, and 0.05 m/s, respectively, and the 3D position errors were 54 m, 45 m, 52 m, and 49 m. Therefore, under the given four cases of speed errors, the influence of speed error on the system is not obvious.
The initial position error was selected as 5, 10, 100, and 300 m. Other parameters were set according to the benchmark parameters in Table 1. The simulation results are shown in the following Figure 9 and Table 8. The 3D velocity errors for these four cases were 0.05 m/s, 0.05 m/s, 0.10 m/s, and 0.22 m/s, respectively, and the 3D position errors were 53 m, 44 m, 103 m, and 238 m.
According to Table 6, Table 7 and Table 8, we realize that the GGAN system is sensitive to the initial attitude error and initial position error of the inertial navigation system. On the other hand, the initial velocity error has little effect on the positioning performance of the system.

5.6. Performance Analysis of Inertial Navigation Level

In order to reflect the influence of different inertial navigation levels on the integrated navigation performance, the gyro bias parameters were set at 1, 0.1, 0.01, and 0.001 degree per hour, respectively. Other parameters were set according to the benchmark parameters in Table 1. The simulation results are shown in the following Figure 10 and Table 9. It can be seen from Table 9, even if an inertial navigation system with a gyro bias of 0.1 degree per hour is adopted, the position error of the GGAN system is still controlled within 100 m. The 3D velocity errors for these three cases were 0.04 m/s, 0.05 m/s, and 0.17 m/s, respectively, and the 3D position errors were 28 m, 34 m, and 89 m. Therefore, the GGAN system can effectively reduce the requirements of inertial navigation device level under the condition of satisfying certain accuracy.

6. Conclusions

With the development of antisubmarine and submarine detection technology, in view of the role of submarines in the future information warfare, working out how to realize the submarine’s passive, autonomous, long-voyage, and high-precision navigation and positioning has become an urgent problem. The GGAN system makes it possible to realize high-precision positioning of a submarine over a long voyage time. The integrated navigation method of inertial/strapdown gravity gradiometer based on quaternion unscented Kalman filter is proposed in this paper, which can effectively restrain the divergence of inertial navigation error. The experimental results show that if the GGAN system is set according to the benchmark parameters, and after about 10 h of simulation the attitude error is less than 1 arcmin and the velocity error is less than 0.1 m/s, the position error is controlled within 33 m. With the increase of the measurement update period, the positioning error of the system gradually increases. When the measurement update period is less than 180 s, the position error in other directions is within 100 m, except for the longitude direction of 127.42 m. When the measurement noise is 1 E, the position error of longitude direction is more than 1000 m. The database noise has a great impact on the performance of the GGAN system. The noise of 0.001 E database is better than that of 0.01 E database, and the positioning performance of the system is mostly better than one order of magnitude. The GGAN system is sensitive to the initial attitude error and initial position error of the inertial navigation system, but the initial velocity error has little effect on the positioning performance of the system. Using GGAN, the system can reduce the requirements for the quality of the used inertial sensor. Therefore, the integrated navigation system based on gravity gradient can provide reliable and high-precision initial navigation parameters for other weapon systems, and provide a strong guarantee for the effectiveness of the entire combat platform.

Author Contributions

Data curation: L.C.; resources: D.G.; software: D.G., L.C., B.H. and X.L.; writing—original draft: D.G.; writing—review and editing, D.G., F.Q. and L.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China (No. 61703419, No. 61873275, No. 61703419).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data sharing not applicable. No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

The normalized gravitational potential function V ( r , θ , λ ) can be expressed in the form of higher-order spherical harmonic function as follows:
V ( r , θ , λ ) = G M a n = 0 N m = 0 n ( a r ) n + 1 P ¯ n , m ( cos θ ) [ C ¯ n , m cos ( m λ ) + S ¯ n , m sin ( m λ ) ]
Here, GM is the gravitational constant multiplied by the Earth’s mass, where θ = π 2 ψ , and r , ψ , λ represent the geocentric distance, the latitude, and longitude, and a is the radius of the Earth. N is the maximum value of the degree n and order m . P ¯ n m ( cos θ ) is the normalized Legendre associated function of the degree n and order m . C ¯ n , m and S ¯ n , m are the harmonic geopotential coefficients with degree n and order m .
According to [27], the calculation formula of gravity gradients reads as follows:
V x i x j = G M 2 1 r x i x j + G M a 2 n = 2 N + 1 m = 0 n ( a r ) n + 1 P ¯ n , m ( cos θ ) [ C ¯ n , m x i x j cos ( m λ ) + S ¯ n , m x i x j sin ( m λ ) ]
where C ¯ n , m x i x i and S ¯ n , m x i x i can be obtained from [28]. The specific calculation formula of 2 1 r / x i x j is as follows:
2 1 r x i x j = 3 x i 2 r 2 r 5 if   i = j 3 x i x j r 5 if   i j
The transformation of the gravity gradient tensor matrix between the body coordinate system and the navigation coordinate system is as follows:
V b = C n b V n ( C n b ) T = C n b V n C b n

References

  1. Moryl, J.; Rice, H.; Shinners, S. The universal gravity module for enhanced submarine navigation. In Proceedings of the IEEE 1998 Position Location and Navigation Symposium (Cat. No.98CH36153), Palm Springs, CA, USA, 20–23 April 1996; pp. 324–331. [Google Scholar]
  2. Vajda, S.; Zorn, A. Survey of existing and emerging technologies for strategic submarine navigation. In Proceedings of the IEEE 1998 Position Location and Navigation Symposium (Cat. No.98CH36153), Palm Springs, CA, USA, 20–23 April 1996. [Google Scholar]
  3. Rice, H.; Mendelsohn, L.; Aarons, R.; Mazzola, D. Next generation marine precision navigation system. In Proceedings of the IEEE 2000. Position Location and Navigation Symposium (Cat. No.00CH37062), San Diego, CA, USA, 13–16 March 2000. [Google Scholar]
  4. Hays, K.M.; Schmidt, R.G.; Wilson, W.A.; Campbell, J.D.; Gokhale, M.P. A submarine navigator for the 21st century. In Proceedings of the 2002 IEEE Position Location and Navigation Symposium (IEEE Cat. No.02CH37284), Palms Springs, CA, USA, 15–18 April 2002. [Google Scholar]
  5. Karshakov, E.; Pavlov, B.; Papusha, I.; Tkhorenko, M. Promising Aircraft Navigation Systems with Use of Physical Fields: Stationary Magnetic Field Gradient, Gravity Gradient, Alternating Magnetic Field. In Proceedings of the 2020 27th Saint Petersburg International Conference on Integrated Navigation Systems (ICINS), Saint Petersburg, Russia, 25–27 May 2020; pp. 1–9. [Google Scholar]
  6. Karshakov, E. Iterated extended Kalman filter for airborne electromagnetic data inversion. Explor. Geophys. 2020, 51, 66–73. [Google Scholar] [CrossRef]
  7. Peshekhonov, V.G.; Sokolov, A.V.; Zheleznyak, L.K.; Bereza, A.D.; Krasnov, A.A. Role of Navigation Technologies in Mobile Gravimeters Development. Gyroscopy Navig. 2020, 11, 2–12. [Google Scholar] [CrossRef]
  8. Volkovitskii, A.K.; Karshakov, E.V.; Tkhorenko, M.Y.; Pavlov, B.V. Application of Magnetic Gradiometers to Control Magnetic Field of a Moving Object. Autom. Remote Control 2020, 81, 333–339. [Google Scholar] [CrossRef]
  9. Soken, H.E.; Hajiyev, C. UKF-Based Reconfigurable Attitude Parameters Estimation and Magnetometer Calibration. IEEE Trans. Aerosp. Electron. Syst. 2012, 48, 2614–2627. [Google Scholar] [CrossRef]
  10. Junji, G.; Xingle, Z.; Wubing, Z. Accurate calculating method of marine three-component geomagnetic field in shipboard measurement. In Proceedings of the 2015 Chinese Automation Congress (CAC), Wuhan, China, 27–29 November 2015; pp. 1504–1507. [Google Scholar]
  11. Tippenhauer, N.O.; Pöpper, C.; Rasmussen, K.B.; Capkun, S. On the requirements for successful GPS spoofing attacks. In Proceedings of the 18th ACM Conference on Computer Supported Cooperative Work & Social Computing; Association for Computing Machinery (ACM), Chicago, IL, USA, 17–21 October 2011; pp. 75–86. [Google Scholar] [CrossRef]
  12. Affleck, C.; Jircitano, A. Passive gravity gradiometer navigation system. In Proceedings of the IEEE Symposium on Position Location and Navigation. A Decade of Excellence in the Navigation Sciences, Las Vegas, NV, USA, 20 March 1990; pp. 60–66. [Google Scholar]
  13. Yan, Z.; Ma, J.; Tian, J. Accurate Aerial Object Localization Using Gravity and Gravity Gradient Anomaly. IEEE Geosci. Remote Sens. Lett. 2015, 12, 1214–1217. [Google Scholar] [CrossRef]
  14. Evstifeev, M.I. The state of the art in the development of onboard gravity gradiometers. Gyroscopy Navig. 2017, 8, 68–79. [Google Scholar] [CrossRef]
  15. Evstifeev, M.; Elektropribor, J.C.C. Dynamics of Onboard Gravity Gradiometers. Giroskopiya Navig. 2019, 27, 69–87. [Google Scholar] [CrossRef]
  16. Tang, J.; Hu, S.; Ren, Z.-Y.; Chen, C. Localization of Multiple Underwater Objects with Gravity Field and Gravity Gradient Tensor. IEEE Geosci. Remote Sens. Lett. 2018, 15, 247–251. [Google Scholar] [CrossRef]
  17. Boddice, D.; Metje, N.; Tuckwell, G. Quantifying the effects of near surface density variation on quantum technology gravity and gravity gradient instruments. J. Appl. Geophys. 2019, 164, 160–178. [Google Scholar] [CrossRef] [Green Version]
  18. Richeson, J.A. Gravity gradiometer aided inertial navigation within non-GNSS environments. In Proceedings of the 20th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS 2007), Fort Worth, TX, USA, 25–28 September 2007. [Google Scholar]
  19. Sun, X.; Chen, P.; Macabiau, C.; Han, C. Autonomous orbit determination via kalman filtering of gravity gradients. IEEE Trans. Aerosp. Electron. Syst. 2016, 52, 2436–2451. [Google Scholar] [CrossRef] [Green Version]
  20. Lee, J.; Kwon, J.H.; Yu, M. Performance Evaluation and Requirements Assessment for Gravity Gradient Ref-erenced Navigation. Sensors 2015, 15, 16833–16847. [Google Scholar] [CrossRef] [Green Version]
  21. Liu, F.M.; Li, F.M. Application of self-adaptive artificial physics optimized particle filter in INS/gravity gradient aided nav-igation. In Proceedings of the 2018 IEEE International Conference on Mechatronics and Automation (ICMA), Changchun, China, 5–8 August 2018; pp. 1070–1075. [Google Scholar]
  22. Liu, F.; Jing, X. INS/gravity gradient aided navigation based on gravitation field particle filter. Open Phys. 2019, 17, 709–718. [Google Scholar] [CrossRef]
  23. Gleason, D.M. Passive airborne navigation and terrain avoidance using gravity gradiometry. J. Guid. Control Dyn. 1995, 18, 1450–1458. [Google Scholar] [CrossRef]
  24. Degregoria, A. Gravity Gradiometry and Map Matching: An Aid to Aircraft Inertial Navigation Systems. Master’s Thesis, Air Force Institute of Technology, Wright-Patterson AFB, OH, USA, 2010. [Google Scholar]
  25. Zhang, M.; Li, K.; Hu, B.; Meng, C. Comparison of Kalman Filters for Inertial Integrated Navigation. Sensors 2019, 19, 1426. [Google Scholar] [CrossRef] [Green Version]
  26. Pavlis, N. An Earth Gravitational model to degree 2160: EGM08. In Proceedings of the 2008 General Assembly of the European Geosciences Union, Vienna, Austria, 13–18 April 2008. [Google Scholar]
  27. Chen, C.; Ouyang, Y.; Bian, S. Spherical Harmonic Expansions for the Gravitational Field of a Polyhedral Body with Poly-nomial Density Contrast. Surv. Geophys. 2019, 40, 197–246. [Google Scholar] [CrossRef]
  28. Petrovskaya, M.S.; Vershkov, A.N. Construction of spherical harmonic series for the potential derivatives of arbitrary or-ders in the geocentric Earth-fixed reference frame. J. Geod. 2010, 84, 165–178. [Google Scholar] [CrossRef]
Figure 1. (a) Gravity gradient aided navigation (GGAN) attitude error; (b) pure inertial calculation attitude error.
Figure 1. (a) Gravity gradient aided navigation (GGAN) attitude error; (b) pure inertial calculation attitude error.
Sensors 21 00829 g001
Figure 2. (a) GGAN speed error; (b) pure inertial calculation speed error.
Figure 2. (a) GGAN speed error; (b) pure inertial calculation speed error.
Sensors 21 00829 g002
Figure 3. (a) GGAN position error; (b) pure inertial calculation position error.
Figure 3. (a) GGAN position error; (b) pure inertial calculation position error.
Sensors 21 00829 g003
Figure 4. (a) Attitude error curve; (b) velocity error curve; (c) position error curve.
Figure 4. (a) Attitude error curve; (b) velocity error curve; (c) position error curve.
Sensors 21 00829 g004aSensors 21 00829 g004b
Figure 5. (a) Attitude error curve; (b) velocity error curve; (c) position error curve.
Figure 5. (a) Attitude error curve; (b) velocity error curve; (c) position error curve.
Sensors 21 00829 g005
Figure 6. (a) Attitude error curve; (b) velocity error curve; (c) position error curve.
Figure 6. (a) Attitude error curve; (b) velocity error curve; (c) position error curve.
Sensors 21 00829 g006
Figure 7. (a) Attitude error curve; (b) velocity error curve; (c) position error curve.
Figure 7. (a) Attitude error curve; (b) velocity error curve; (c) position error curve.
Sensors 21 00829 g007
Figure 8. (a) Attitude error curve; (b) velocity error curve; (c) position error curve.
Figure 8. (a) Attitude error curve; (b) velocity error curve; (c) position error curve.
Sensors 21 00829 g008
Figure 9. (a) Attitude error curve; (b) velocity error curve; (c) position error curve.
Figure 9. (a) Attitude error curve; (b) velocity error curve; (c) position error curve.
Sensors 21 00829 g009
Figure 10. (a) Attitude error curve; (b) velocity error curve; (c) position error curve.
Figure 10. (a) Attitude error curve; (b) velocity error curve; (c) position error curve.
Sensors 21 00829 g010
Table 1. Benchmark parameter settings.
Table 1. Benchmark parameter settings.
Benchmark ParameterMeasurement Update PeriodMeasurement NoiseDatabase NoiseInitial ErrorGyro Bias
AttitudeVelocityPosition
Value60 s0.01 E0.001 E 0.5 1 m/s10 m 0 . 01 / h
Table 2. Comparison of long voyage performance.
Table 2. Comparison of long voyage performance.
ModeAttitude Errors (arcmin)Velocity Error (m/s)Position Error (m)
PitchRollYawEastNorthUpLatitudeLongitudeHeight
Pure inertial calculation0.810.713.031.691.76/26402258/
GGAN0.280.180.720.060.090.03332917
Table 3. Simulation results under different measurement update periods.
Table 3. Simulation results under different measurement update periods.
Measurement Update Periods (s)Attitude Errors (arcmin)Velocity Error (m/s)Position Error (m)
PitchRollYawEastNorthUp3DLatitudeLongitudeHeight3D
300.240.120.890.020.030.010.0334181843
600.240.120.890.030.040.020.0540291552
900.350.181.610.040.050.030.0734491562
1800.550.272.260.090.070.050.126512736147
Table 4. Simulation results under different measurement noise.
Table 4. Simulation results under different measurement noise.
Measurement Noise (E)Attitude Errors (arcmin)Velocity Error (m/s)Position Error (m)
PitchRollYawEastNorthUp3DLatitudeLongitudeHeight3D
0.0010.080.050.350.030.030.020.0522291340
0.010.090.050.370.030.040.020.0525341445
0.10.310.141.070.120.100.140.2112617766227
10.550.301.340.890.481.611.9071810878071533
Table 5. Simulation results under different database noise.
Table 5. Simulation results under different database noise.
Database Noise (E)Attitude Errors (arcmin)Velocity Error (m/s)Position Error (m)
PitchRollYawEastNorthUp3DLatitudeLongitudeHeight3D
0.0010.240.120.890.030.040.030.0547262760
0.011.630.836.250.891.580.431.8613795592491509
Table 6. Simulation results of different initial attitude errors.
Table 6. Simulation results of different initial attitude errors.
Initial Attitude Errors
(deg)
Attitude Errors (arcmin)Velocity Error (m/s)Position Error (m)
PitchRollYawEastNorthUp3DLatitudeLongitudeHeight3D
0.10.240.120.930.030.050.020.0639242051
0.20.300.141.120.040.040.020.0644191851
0.50.120.060.420.020.050.020.0632361251
12.381.179.290.310.390.180.6314450178577
Table 7. Simulation results under different initial velocity errors.
Table 7. Simulation results under different initial velocity errors.
Initial Velocity Errors
(m/s)
Attitude Errors (arcmin)Velocity Error (m/s)Position Error (m)
PitchRollYawEastNorthUp3DLatitudeLongitudeHeight3D
10.220.110.800.030.040.020.0545222054
20.190.090.710.030.040.020.0535231645
30.230.110.830.030.050.020.0639271952
50.240.120.890.030.030.020.0538241849
Table 8. Simulation results under different initial position errors.
Table 8. Simulation results under different initial position errors.
Initial Position Errors
(m)
Attitude Errors (arcmin)Velocity Error (m/s)Position Error (m)
PitchRollYawEastNorthUp3DLatitudeLongitudeHeight3D
50.200.100.730.030.040.020.0543231953
100.170.080.620.030.040.020.0533231544
1000.540.272.080.060.070.040.10706044103
3001.250.634.890.140.110.130.22151150106238
Table 9. Simulation results of different inertial navigation levels.
Table 9. Simulation results of different inertial navigation levels.
Gyro Bias
( / h )
Attitude Errors (arcmin)Velocity Error (m/s)Position Error (m)
PitchRollYawEastNorthUp3DLatitudeLongitudeHeight3D
0.0010.080.050.270.030.020.020.0419151328
0.010.090.050.330.030.040.020.0523201434
0.10.250.140.970.100.140.030.1769541689
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Gao, D.; Hu, B.; Chang, L.; Qin, F.; Lyu, X. An Aided Navigation Method Based on Strapdown Gravity Gradiometer. Sensors 2021, 21, 829. https://doi.org/10.3390/s21030829

AMA Style

Gao D, Hu B, Chang L, Qin F, Lyu X. An Aided Navigation Method Based on Strapdown Gravity Gradiometer. Sensors. 2021; 21(3):829. https://doi.org/10.3390/s21030829

Chicago/Turabian Style

Gao, Duanyang, Baiqing Hu, Lubin Chang, Fangjun Qin, and Xu Lyu. 2021. "An Aided Navigation Method Based on Strapdown Gravity Gradiometer" Sensors 21, no. 3: 829. https://doi.org/10.3390/s21030829

APA Style

Gao, D., Hu, B., Chang, L., Qin, F., & Lyu, X. (2021). An Aided Navigation Method Based on Strapdown Gravity Gradiometer. Sensors, 21(3), 829. https://doi.org/10.3390/s21030829

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop