An Aided Navigation Method Based on Strapdown Gravity Gradiometer

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.


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.

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 n b 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: 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 Other parts of the basic SINS equations are . . . . .
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 ω b ib ; R M and R N represent the meridian and prime vertical radius of the Earth, as shown in the following formula: where R e = 6378137m, e = 0.0818. The differential equation of gyro constant drift is as follows: .
where ω b ib is the actual gyro output with gyro drift ε; η gv and η gu are zero mean Gaussian white noise with covariance σ 2 gv and σ 2 gu , respectively. The output equation of the accelerometer reads as follows: .
where f b is the actual accelerometer output with accelerometer bias ∇, η av and η au are zero mean Gaussian white noise with covariance σ 2 av and σ 2 au , 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.

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 ω a ia . 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 where C i a is the transformation matrix between coordinate system a and inertial coordinate system i. It can be obtained by time derivation of Equation (15) .
if ω a ia = [ ω 1 ω 2 ω 3 ] T , the antisymmetric matrix Ω a ia can be expressed as The time derivative of Equation (16) leads to ..
In inertial space the following relation holds ..
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 .
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 a 2 − r a 1 . According to Equation (20), we can get Since the accelerometer is attached to coordinate system a, then .
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 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 and 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 Ω a ia and angular acceleration . Ω a ia . 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 Let By substituting each component into Equation (24), we end up with The six independent components in expression (28) are as follows: Compared with Equation (25), Equation (27) does not need to estimate . Ω a ia , 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 where C b n is the transformation matrix between the navigation coordinate system and the body coordinate system, it can be obtained by attitude calculation, and Ω b ib 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.

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 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 At the kth moment, the system noise is Gaussian white noise w k ∼ N(0Q k ), and the measurement noise is Gaussian white noise v k ∼ N(0R k ). The nonlinear system of gravity gradient aided positioning can be expressed by the following formula: 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 Calculate the sigma point from the following formula: where the subscript ":, i" denotes the column i of the matrix. Each sigma point can be propagated through the system model: After propagation, the state estimation and its error covariance are as follows: The observation update process of UKF generates new sigma points by the following formula: The sigma point and mean observed innovation can be solved by the following formula: and the covariance of the observed innovation is Finally, the UKF Kalman gain, state vector update, and error covariance update can be calculated as follows: When the system noise and measurement noise are Gaussian white noise, GGAN can be realized by the standard UKF localization algorithm mentioned above.

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.

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 Figures 1-3 and Table 2 below.
In order to intuitively reflect the long voyage characteristics of inertial/gravity gradi-ent 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 Figures 1-3 and Table 2 below.   ent 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 Figures 1-3 and Table 2 below.   ent 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 Figures 1-3 and Table 2 below.    It can be seen from Figures 1-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.

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.  It can be seen from Figures 1-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.

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.

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

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.

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 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.

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.

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.

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 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 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  According to Tables 6-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.

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.  According to Tables 6-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.

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.

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 lon-

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.  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: 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 m n (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: where C x i x i n,m and S x i x i n,m can be obtained from [28]. The specific calculation formula of ∂ 2 1 r /∂x i x j is as follows: The transformation of the gravity gradient tensor matrix between the body coordinate system and the navigation coordinate system is as follows: