Abstract
With the development of high-precision inertial navigation systems, the deflection of vertical (DOV), gravity disturbance, is still one of the main error sources that restrict navigation accuracy. For the DOV compensation of the Strapdown Inertial Navigation System (SINS) problem, the influences of the calculation degree of the spherical harmonic coefficient and the calculation error of the DOV on the compensation effect were studied. Based on the SINS error model, the error propagation characteristics of the DOV in SINS were analyzed. In addition, the high-precision global gravity field spherical harmonic model EIGEN-6C4 was established and the influence comparative analysis of the calculation degree of the spherical harmonic coefficient on the DOV compensation of SINS in different regions was carried out. Besides, the influence of the calculation error of the DOV on the compensation was emphatically analyzed. Finally, the vehicle experiment verified the feasibility of compensation in SINS based on the gravity field spherical harmonic model. The simulation and experiment results show that it is necessary to consider the influence of the calculation degree and the calculation error of the DOV on the compensation for long-time high-precision SINS with the position accuracy of 0.3 nm/h, while the SINS with general requirements for position accuracy can ignore the impact.
1. Introduction
Inertial navigation system (INS) is the position reckoning system that uses its inertial measurement units (IMU) to calculate the attitude, speed, and position information of the carrier [1]. It is widely used in various military and civil fields because of the independence without the external information and high reliability. In order to facilitate the navigation calculation, the normal gravity field defined in the reference ellipsoid is typically used instead of the global true gravity field for navigation calculation. That definitely ignores the difference between the normal gravity vector and the true gravity vector, which is the gravity disturbance vector [2, 3]. With the continuous improvement of the accuracy of inertial devices, the gravity disturbance compensation on the INS is particularly significant [4].
The deflection of vertical refers to the difference between the vertical line of the geoid and the normal line of the reference ellipsoid, that is, the deviation of the true gravity vector from the normal gravity vector [5]. The horizontal component of the deflection of vertical (DOV) directly determines the magnitude of the horizontal gravity disturbance. According to the INS error analysis, the DOV directly affects the initial alignment error and speed error and thus affects the position error after integration [6]. Jekeli [7] and Jekeli et al. [8] firstly studied the problem of gravity disturbance compensation for INS. Hanson [9] analyzed the influence of the initial alignment of the DOV compensation. Wang et al. [10] started with the error model of the INS and focused on analyzing the error propagation characteristics of the horizontal component of the gravity disturbance vector in the navigation system. Jin and Bian [11] simulated the INS position errors caused by three kinds of different characteristics of gravity disturbances. In recent years, some scholars had used the high-resolution spherical harmonic model (SHM) of the gravity field to compensate and analyze the DOV. Wang et al. [12] studied the accuracy of the truncated spherical harmonic model for real-time compensation of INS and proposed a simplified two-dimensional second-order polynomial model for the spherical harmonic model of the gravity field. Chang et al. [13] used the EIGEN-6C4 gravity field spherical harmonic model to compensate for the initial alignment and navigation calculation of the INS. Wu et al. [14] had researched on the calculation of the minimum update rate of gravity using the SHM in the shipboard INS. The above compensation analysis basically uses the highest degree or any degree of the spherical harmonic model for compensation analysis [15–17], which greatly reduces the calculation efficiency and does not consider the influence of the error in the INS solution position on the compensation effect during the compensation. For the INS with an accuracy of less than 100 meters, the influence of the DOV on the INS cannot be ignored, especially for the long-time INS. The position error caused by DOV definitely reduces the performance of INS. Consequently, it is necessary to analyze and compensate the influence of DOV on INS.
Starting from the SINS error model and combining the EIGEN-6C4 ultra-high-precision gravity field spherical harmonic model, this paper analyzes the compensation effect of the DOV calculated by the spherical harmonic model of different degrees. The compensation accuracies of the DOV at the true position and the INS solution position are compared and analyzed. The proposed methodology can be applied in long-time high-precision INS, for instance, the land vehicle INS and the submarine vehicle.
The contents are organized as follows: Section 1 presents the introduction; the definition of the DOV and the expression of the gravity disturbance vector (GDV) are briefly reviewed in Section 2; Section 3 describes the construction of the spherical harmonic model in the gravity field and the calculation of DOV; the analysis of DOV compensation position error of INS and compensation process are described in Section 4; in Section 5, compensation effects of SHM in different degrees in the general area and abnormal area are reported; the influence of calculation error of DOV on SINS position error and vehicle experiment are presented in Section 6; finally, conclusion is presented in Section 7.
2. The Deflection of Vertical and Gravity Disturbance Vector
The geoid is a model that describes the true shape of the Earth, which reflects the information of the global true gravity field [18–20]. Due to uneven mass distribution inside the Earth and uneven surface, the geoid is extremely difficult to describe by a mathematical model. To facilitate mathematical calculations, a reference ellipsoid model close to the geoid has been established [21]. The gravity field on the reference ellipsoid model is called the normal gravity field, and the difference between the real gravity vector and the normal gravity vector is called the gravity disturbance vector [22]. The deflection of vertical refers to the difference between the vertical line of the geoid and the normal line on the reference ellipsoid, that is, the deviation of the true gravity vector from the normal gravity vector, as shown in Figure 1. The component of the DOV in the east-west direction is , and the north-south component is , expressed by equations (1)–(3). The schematic diagram of the gravity disturbance vector in space is shown in Figure 2.where denotes the normal gravity value and means the gravity anomaly. , , and are the eastward, northward, and upward of GDV in navigation frame n, respectively.
For coordinate system definition, i means the Geocentric Inertial Coordinate System; e means the Earth-Centered Earth-Fixed Frame; n denotes the Navigation Frame with East-North-Up; select the local geographic coordinate frame; b means the body frame with Right-Forward-Upward, three axes toward the Right, Front, and Up, respectively, and the diagram of coordinate frames is shown in Figure 3.
3. Gravity Field Spherical Harmonic Model
According to the law of universal gravitation, the gravity potential of a point on the Earth’s surface is the sum of the gravitational potential and the centripetal force potential, as shown in equation (4). In the spherical coordinate frame, the gravitational potential can be represented by spherical harmonic functions, and the gravitational potential and centripetal force can be represented by equations (5) and (6), respectively. U is the gravity potential of the point; W is the gravitational potential calculated by spherical harmonic function; Q is the centripetal force potential:where G represents the gravitational constant; M represents the mass of the Earth; denotes the extra latitude of the calculation point, also known as the polar distance; denotes the longitude; r represents the distance between the calculation point and the center of the reference ellipsoid; is the Earth’s radius; N represents the maximum degree in the spherical harmonic model; and are the harmonic geopotential coefficients with degree n order m; is the rotation angular velocity of Earth; are the normalized Legendre associated functions, which can be calculated by equations (14)–(17).
Gravity is the first-order derivative of the gravity potential, as shown in the following formula:where , , and denote the unit vectors along the radial, epicenter, and longitude directions, respectively, corresponding to the up, south, and east directions on the ellipsoid, where the upward is the center of the ellipsoid pointing to the mass point direction rather than up of geographic direction.
In the spherical coordinate frame, taking derivatives of equation (4) with respect to , , and r, we can obtain
From the relationship between the spherical coordinates and the Cartesian coordinate frame, it can be calculated that the components of true gravity in the “East-North-Up” coordinate system at the mass point are shown in equations (11)–(13). The aforementioned process is the spherical harmonic model (SHM). It is obvious that the gravity can be calculated by gravity field SHM as long as the position information was given.where , , and are the gravity components in the frame, where the upward direction is from the center of the sphere points to the unit particle.
In the INS, the gravity components in the geographic coordinate frame were usually used in navigation calculation. The difference between the geographic latitude L and the geocentric latitude is , shown in Figure 4. Then the horizontal component of the true gravity in the geographic coordinate system can be obtained by the rotary transformation, which is the horizontal gravity disturbance vector. Consequently, the DOV can be calculated by formulas (14) and (15).where , , and are the gravity components in the navigation frame n.
3.1. Legendre Recursive Calculation
(1)Initial value is(2)When , ,(3)When , ,(4)When , ,
From the northward component calculation formula (8), it can be seen that solving the northward gravity disturbance requires normalizing the derivative of the Legendre function. From [15], the recursive formula for the derivative of the Legendre function is as follows.
3.2. Legendre Derivatives Recursive Calculation
(1)When ,(2)When ,
4. Analysis of DOV Compensation Position Error of INS
As we all know, the SINS specific force equation describes the relationship between the motion acceleration of carrier and the specific force, which is shown as follows:where is the motion acceleration of carrier; is the specific force of carrier measured by accelerometer; represents the component of the Earth rotation angular velocity in the navigation coordinate frame; is the angular velocity of the navigation frame relative to the Earth coordinate frame; is the gravity acceleration in navigation coordinate frame.
Differentiate the specific force equation above to obtain the velocity error equation as follows:where is the attitude error angle; and are the Earth rotation angular velocity error and navigation coordinate system update error affected by the position error and velocity error, which can be calculated by equations (24) and (25), respectively; is the zero bias error of the accelerometer; is the velocity error in navigation frame.
The velocity error equation can be sorted into the form of a matrix:
The matrices involved in (26) are given by
The attitude error equation of SINS is as follows:
The attitude error equation can be sorted into the form of a matrix:
The matrices involved in (29) are given bywhere denotes the antisymmetric matrix of .
The position error equation can be expressed as follows:
The matrices involved in (31) are given by
Equations (26), (28a), (28b), and (31) are the SINS error equations when the gravity disturbance error is considered. From the equations, it can be seen that the DOV first enters the velocity error channel and affects the velocity error of the horizontal channel, thereby affecting the SINS position error and body attitude error after integration. The schematic figure of the influence of the DOV on the SINS is shown in Figure 5. In the pure inertial navigation system, the height channel is divergent, so the coupling errors of the height channel can be ignored. This study is mainly about the compensation analysis of the DOV to the horizontal position error of SINS and the position error equation (31) is rewritten into the form of the horizontal channel as (33). The process of DOV compensation is shown in Figure 6.
5. Compensation Effects of SHM in Different Degrees
Equations (8)–(10) are the recursive spherical harmonic series of the SHM to calculate the DOV. This study utilizes the latest EIGEN-6C4 high-precision SHM of gravity field published by the German Geoscience Center in 2012. The degree of the model can reach up to 2190. It can be known from (5) that the degree N of spherical harmonic coefficients and will directly affect the accuracy of DOV calculation. This section focuses on analyzing the influence of the horizontal position error of the DOV compensation calculated by the SHM of different degrees.
5.1. DOV of the Two Areas
In order to clearly indicate the influences of the different model degrees, the coefficients of the 360-degree, 900-degree, and 2190-degree SHMs are selected to calculate the DOV for the general area and the abnormal area, where the DOV changes sharply. The general area is selected in North China province, N 34°18′45′′∼38°18′45″, E 109°7′30″–113°7′30″, and the abnormal area is selected in the Himalayas, N 26°18′45″–30°18′45″, E 89°7′30′′∼93°7′30′′. The calculation resolution of DOV is set to 5′ × 5′.
It can be seen from Figure 7 that the sizes of the DOV component and are basically similar in the selected general area. As the degree of the SHM increases, the calculated value of the DOV component in the region shows a certain increasing trend. The calculation result of the DOV in 2190-degree model is more obvious than that in the 360-degree model. It can be seen from Figure 8 that the size of has a larger value than the component in the selected abnormal area. As the degree of the SHM increases, the calculated value of the DOV increases more obviously. It can be known from the analysis that the calculation accuracy of the 2190-degree model is higher than that of the 360-degree model, but as the calculation accuracy increases, its calculation speed will decrease sharply. The calculation efficiency will be greatly reduced when there is plenty of points in a certain area. Therefore, it is necessary to analyze the influence of the spherical harmonic coefficients of different degrees to calculate the DOV to compensate for the SINS position error.
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f)
5.2. Simulation Analysis of DOV Compensation of SINS Position Error
It can be known from the analysis in Section 4 that the DOV enters the SINS speed channel first which affects the horizontal speed error and thus the position error after integration. This section mainly analyzes the influence of the DOV compensation SINS position error calculated by the SHM of three different degrees, and the compensation process is shown in Figure 6. Firstly, utilize the EIGEN-6C4 SHM of 360, 900, and 2190 degrees to calculate the DOV on the carrier’s trajectory, respectively. Then, the interpolated 2190-degree calculated data was accumulated in the normal gravity field to simulate the true gravity field data and generate the inertial measurement unit (IMU) data by simulation. Finally, DOV calculation data of the 360, 900, and 2190 degrees are used to compensate for the SINS position error.
The simulation conditions were set as follows: the speed was 40 m/s; the trajectories were set to N 34°18′45″–38°12′45″, E 109°7′30″, and N 26°18′45″–30°12′45″, E 89°7′30″, in general and abnormal areas, respectively. The gyroscope constant drift was set as 0.001°/h, and the random walk coefficient was 0.0002°/h1/2; accelerometer constant zero bias was 10 µg, and the zero bias white noise was 5 µg. The data sampling interval of the IMU was 0.01 s. The trajectory was calculated by fourth-degree Runge Kuta, and the calculation period was 0.005 s. The attitude and speed update of the INS were solved by the optimized twin-sample algorithm, the solution period was 0.02 s, and the simulation time was 10800 s.
The DOV of the trajectory in the two areas are shown in Figures 9 and 10, respectively. It can be seen from Figures 9 and 10 that the calculated value of 2190 degrees is basically consistent with the calculated value of 360 degrees and 900 degrees, and the difference between the calculated values is no more than 52″. Compared with the trajectory in the abnormal area, the 2190-degree calculated data is quite different from the 900-degree and 360-degree calculated DOV data, and the maximum difference can reach nearly 30''.
(a)
(b)
(a)
(b)
The DOV data calculated by the 360-degree, 900-degree, and 2190-degree spherical harmonic coefficients are superimposed into the normal gravity field, respectively, as the real gravity field for navigation calculation. The position error curves after 3 sets of DOV compensation are shown in Figure 11. The SINS position error after DOV compensation is significantly reduced compared to the case without compensated state, especially in the abnormal area, as shown in Figure 11. In the general area, the compensation effects of the three sets of DOV data are basically similar, and the maximum longitude error difference is only 130 m. In the abnormal area, the compensation effects of the three sets of DOV data are generally similar in latitude error and have a large difference in longitude error, with a maximum difference of nearly 300 m. The total maximum statistical data after compensation are shown in Table 1.
(a)
(b)
The simulation compensation results and statistical data show that, for long-endurance high-precision INS with a position accuracy of 0.2 nm/h, the influence of the calculation degree of the SHM on the position accuracy of the DOV compensation needs to be considered. The vertical deviation data calculated by the 2190-degree spherical harmonic coefficient are required for compensation. For the INS with general position accuracy requirements, the DOV data compensation of 360-degree calculation can satisfy the accuracy requirements.
6. The Influence of Calculation Error of DOV on SINS Position Error
In the real SINS working process, the calculated position information by SINS contains errors. In the preliminary studies, all methods of calculating the DOV based on SHM are basically the real position calculations. During practical engineering applications, the DOV data can only be calculated from the position output by SINS, so it is necessary to consider the influence of the DOV calculation error caused by the position error on the SINS compensation effect. This section mainly analyzed the effect of compensating the SINS position error by the actual position and the DOV calculated by the SINS calculated position through simulation experiments.
6.1. Simulation Analysis
The simulation condition settings in this section are consistent with Section 5.2. The true position on the trajectory in the general area and the abnormal area and the DOV data calculated by the SINS calculated position are shown in Figures 12 and 13, respectively. It can be seen from Figures 12 and 13 that the trends of the DOV calculated at the two positions are basically the same. In the general area, the calculation results of the two positions are generally the same, and the maximum difference is only about 2”. While, in the abnormal area, the calculated position of SINS was different greatly from the calculated value at the true position, the value is more obvious than . The maximum difference can reach nearly 10”.
Compensate the DOV data calculated at the true position and the DOV data calculated at the SINS calculated position into the SINS navigation calculations, respectively. Figures 14 and 15 can be obtained; Figure 14 shows the speed error curve after compensation in the two areas. As can be seen from the figure, after the DOV has been compensated, the speed error during the movement is significantly reduced, especially in the abnormal area. The speed error compensation effect at the true position is mainly the same as the calculated position. Figure 15 shows the compensated position error curves in the two areas. It can be seen from Figure 15 that the position error is definitely reduced after the DOV is compensated. In the general area, the compensation effect at the real position was basically the same as the calculated position. Meanwhile, in the abnormal area, the compensation at the two positions has basically the same effect in the latitude error channel. In the longitude error channel, the compensation at the true position is smaller than the compensation at the calculated position, and the position error is smaller. The maximum difference in compensation error between the two positions is close to 130 m. The maximum statistics position errors are shown in Table 2.
(a)
(b)
(a)
(b)
Consequently, in the general area where the DOV changes gently, the influence of the calculation error of the DOV caused by the SINS calculated position error on the position error compensation effect can be ignored. Meanwhile, in the area where the DOV changes abnormally, the INS that requires a position accuracy of 0.3 nm/h needs to consider the influence of the calculation error of the DOV.
6.2. Experiment Analysis
To further verify the reliability of the calculated DOV compensation SINS method based on the SHM of the gravity field, this section mainly conducts experiments on the vehicle-mounted SINS and performs the compensation analysis. SINS was fixedly installed on the experimental vehicle, as shown in Figure 16. The vehicle was equipped with a GPS receiver and used GPS navigation position and speed information as the reference. INS/GPS integrated navigation attitude information was regarded as the attitude reference, and SINS sampling frequency was 200 Hz. GPS provided 1 Hz speed and position reference. The gyroscope constant drifts were 0.0050, 0.0048, and 0.0045°/h and the corresponding random walk coefficients were 0.00046, 0.00044, and 0.00047°/h1/2, respectively. The accelerometer constant biases were 47, 51, and 37 μg. The starting point position was N 34°18′48″, E 109°7′4″ and the ending point was N 35°28′42″, E 110°27′42″, The true GPS navigation trajectory and SINS calculated trajectory are shown in Figure 17(a), and the attitude error during the movement is shown in Figure 17(b).
(a)
(b)
The EIGEN-6C4 SHM is used to calculate the DOV on the motion trajectory and compensated for the SINS navigation calculation. The northward speed error and eastward speed error before and after compensation are shown in Figure 18. It can be seen from Figure 18 that, after the DOV is compensated, the northward speed error is definitely reduced, while the eastward speed error is not significantly improved. The position error and the maximum statistical data before and after compensation are shown in Figure 19 and Table 3, respectively. After the DOV is compensated, the latitude position error is definitely reduced, and the longitude position error is not significantly improved. This is because of the coupling effect between the gravity disturbance and the accelerometer drift error that is not eliminated after calibration. When the accelerometer drift and the horizontal gravity disturbance have opposite signs, the position accuracy may be reduced after DOV compensation. However, with the development of rotation modulation technology, the error of IMU can be modulated and compensated. At this time, DOV compensation is more obvious to improve the position accuracy of SINS.
(a)
(b)
(a)
(b)
The land vehicle experiment with different SHM degrees compensation was carried out in this part. The DOV data calculated by the 360-degree, 900-degree, and 2190-degree SHM were compensated in navigation. The results of position errors with compensation are shown in Figure 20 and the maximum position errors are listed in Table 4.
(a)
(b)
It is apparent that the compensation effect of 2190 degrees is relatively better than that of the 360 degrees and that of the 900 degrees. The 2190-degree SHM is the optimal choice to compensate, if the calculation time is not necessary to consider in the offline compensation experiment; otherwise, the 900-degree SHM is also enough for the DOV compensation in navigation.
7. Conclusions
(1)The propagation characteristics of the influence of the DOV on the SINS are embodied as follows: the DOV first enters the speed error channel, affecting the speed error of the horizontal channel, and thus affects the SINS position error and carrier attitude error after integration.(2)In the area where the DOV changes abnormally, the SHM calculation degree has a greater impact on the compensation effect, especially in the longitude error. For long-endurance high-precision SINS with a position accuracy requirement of 0.2 nm/h, it is necessary to consider the influence of the calculation degree on DOV compensation.(3)In the general area, the influence of the calculation error of the DOV on the position accuracy after compensation can be ignored. Meanwhile, in the abnormal area, for inertial navigation systems that require a position accuracy of 0.3 nm/h, the influence of the calculation error of the vertical deviation needs to be considered.
Data Availability
The data used to support the findings of this study are owned by Xi’an Institute of Hi-Tech and so cannot be made freely available. Access to these data will be considered by the corresponding author upon request, with permission of the Research Management Department of Xi’an Institute of Hi-Tech.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was partially supported by the National Natural Science Foundation of China under Grant 61503391 and the China Postdoctoral Science Foundation under Grant 2017M613372.