Next Article in Journal
Unannounced Meals in the Artificial Pancreas: Detection Using Continuous Glucose Monitoring
Next Article in Special Issue
Height Error Correction for Shoe-Mounted Inertial Sensors Exploiting Foot Dynamics
Previous Article in Journal
A New 3D Object Pose Detection Method Using LIDAR Shape Set
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Model of Gravity Vector Measurement Noise for Estimating Accelerometer Bias in Gravity Disturbance Compensation

1
College of Mechatronics and Automation, National University of Defense Technology, Changsha 410073, China
2
Department of Navigation Engineering, Naval University of Engineering, Wuhan 430000, China
*
Author to whom correspondence should be addressed.
Sensors 2018, 18(3), 883; https://doi.org/10.3390/s18030883
Submission received: 3 February 2018 / Revised: 13 March 2018 / Accepted: 15 March 2018 / Published: 16 March 2018
(This article belongs to the Special Issue Inertial Sensors and Systems 2018)

Abstract

:
Compensation of gravity disturbance can improve the precision of inertial navigation, but the effect of compensation will decrease due to the accelerometer bias, and estimation of the accelerometer bias is a crucial issue in gravity disturbance compensation. This paper first investigates the effect of accelerometer bias on gravity disturbance compensation, and the situation in which the accelerometer bias should be estimated is established. The accelerometer bias is estimated from the gravity vector measurement, and a model of measurement noise in gravity vector measurement is built. Based on this model, accelerometer bias is separated from the gravity vector measurement error by the method of least squares. Horizontal gravity disturbances are calculated through EGM2008 spherical harmonic model to build the simulation scene, and the simulation results indicate that precise estimations of the accelerometer bias can be obtained with the proposed method.

1. Introduction

An inertial navigation system (INS) is an instrument which can autonomously determine the ship’s attitude, velocity and position based on its self-contained gyroscopes and accelerometers [1,2]. Because of the self-localization feature, INS is not vulnerable to any type of external interference and has wide application in military fields. However, as a dead reckoning approach, the precision of INS will drift with time due to its inherent error sources. The inherent error sources of INS not only include the errors of inertial sensors, but also include the gravity disturbance [3,4]. In recent years, the significant improvement of inertial sensors has left the gravity disturbance as the most important error source in high precision inertial navigation [5,6]. In the future the inertial-sensors-induced position error would be reduced to only a few meters per hour with cold atom interferometry gyroscopes [6], and the gravity disturbance compensation should be considered [7,8,9,10].
The definition of gravity disturbance is illustrated in Figure 1 [11]. According to the potential theory, gravity vector is the perpendicular line of equipotential surface of gravity. The Earth’s equipotential surface of gravity is very complex, and the equipotential surface of reference ellipsoid model such as WGS-84, is used to approximate the Earth’s equipotential surface of gravity. As shown in Figure 1, g is the true gravity vector of point P and γ is the normal gravity vector of point P, the gravity disturbance vector is the difference between the true gravity vector and the normal gravity vector [11]. The difference in magnitude is the gravity disturbance and the difference in direction is the deflection of vertical (DOV). Due to DOV, there are some projection components of the true gravity vector in the horizontal plan, which are named horizontal gravity disturbance.
Because of the instability of the INS’s vertical channel, INS generally only provides horizontal velocity, latitude and longitude of the ship, and horizontal gravity disturbance is the main error source of gravity-induced error. The study of compensation of horizontal gravity disturbance can be traced back to the sixties of last century, and works have been done on the analysis of error induced by horizontal gravity disturbance. In 1962, Kayton demonstrated that the practical limit of measurement of inertial acceleration is the knowledge of gravitational field, and any differences between the actual field and the model will cause a navigation error and the use of accelerometers which are more sensitive than 10−5 g is probably ineffective unless acquiring detailed knowledge of Earth’s gravity field [12]. Based on the work of Kayton, Levine and Gelb quantitatively analyzed the effect of horizontal gravity disturbances on INS, where the horizontal gravity disturbance is modeled as a first-order Gauss-Markov process and the effect on INS is evaluated with the steady-state solution of the covariance matrix differential equation [13]. After that, based on the covariance analysis method proposed by Levine and Gelb, Jordan analyzed the effect of horizontal gravity disturbance with a new model, named self-consistent statistical model. Different from the first-order Gauss-Markov model built by Levine and Gelb, this model is built based on the potential theory, so this model is more consistent with the true covariance of horizontal gravity disturbance [14]. These works indicate that the horizontal gravity disturbance will affect the accuracy of the initial alignment and velocity calculation. In order to compensate the effect on INS, gravity gradiometers were developed to measure the horizontal gravity disturbance along the ship’s trajectory, such as the universal gravity module (UGM) developed by Lockheed Martin Federal Systems [15]. Nowadays, with the release of ultra-high degree global Earth gravitational models, such as the EGM2008, horizontal gravity disturbances can be precisely calculated based on these spherical harmonic models (SHMs), and the effect of horizontal gravity disturbance on INS can be compensated. The compensation method based on SHM is preferable, because only some software updates of INS are needed to obtain the horizontal gravity disturbance, instead of using a costly gradiometer [16,17]. In this paper, calculation of the horizontal gravity disturbance based on SHM is described in detail.
There is a crucial issue in horizontal gravity disturbance compensation, namely the fact that the compensation effect is associated with the accelerometer bias. Because of the tight coupling between accelerometer bias and horizontal gravity disturbance, the method of compensation is ambiguous, especially when considering the initial alignment and INS calculation synthetically. In 1959, the horizontal gravity component from the mission data was used in the MH-311 system for the Army’s AN/USD-5 surveillance. After three years, the derivative of the MH-311 compensates the horizontal gravity disturbance during self-alignment, while removes such compensation in navigation mode. In 1982, an “improved” DOV compensation procedure was used in Mini-GEANS. In particular, the system is firstly aligned to the local gravity vector and after takeoff, the alignment matrix is rotated to the reference ellipsoid. In our early work [18], the compensation is studied in both the initial alignment and navigation calculation. The inertial navigation experiment demonstrates that compensation only in navigation calculation and not in initial alignment is more preferred. Reference [19] reported a military standard ring laser inertial navigation unit named LN-93E, which is an enhanced derivative of its earlier version, LN-93. One of the reasons for the performance improvement of LN-93E is just the DOV compensation at the align location. In [20], the authors analyze why the DOV compensation in initial alignment is not universal. It is pointed out that this is mainly because the correlations between the accelerometer bias and horizontal gravity disturbance. Unfortunately, the corresponding conclusion in [20] is qualitative and the proposed improved compensation procedure is also case-dependent and not universal.
If the accelerometer bias can be accurately estimated, the compensation method can be worked out, and the effect of accelerometer bias on compensation can be eliminated. INS/GNSS (GNSS, global navigation satellite system) integrated is the usual method to estimate the accelerometer bias [21]. However, the precision of this method is limited by the poor observability of accelerometer bias [22,23], and the estimation is the combination of accelerometer bias and horizontal gravity disturbance. Although a statistical model of the horizontal gravity disturbance can be used in the filter to separate the accelerometer bias from the horizontal gravity disturbance, it’s hard to build a universal statistical model of horizontal gravity disturbance [24].
In some cases, the use of GNSS is not possible and the positioning can only depend on INS. To improve the accuracy of INS in long-endurance navigation, the compensation of gravity disturbance is very necessary. However, the effect of compensation will decrease due to the accelerometer bias. One possible solution to this problem is estimating the accelerometer bias with GNSS before the long-endurance inertial navigation. In this paper, we try to do the estimation in the method of gravity vector measurement [25,26], we attempt to separate the accelerometer bias from the measurement error of gravity vector. The model of measurement noise is built to implement the estimation. The contents are organized as follows: the horizontal gravity disturbance calculation using SHM is described in Section 2. In Section 3, the reason why accelerometer bias influences compensation is analyzed. In Section 4, our proposed method of estimating accelerometer bias is introduced in detail and simulation results are provided in Section 5. Finally, conclusions are drawn in Section 6.

2. Horizontal Gravity Disturbance and Spherical Harmonic Model

2.1. Definition of Horizontal Gravity Disturbance

DOV has two components as shown in Figure 2, where g is the true gravity vector, and γ is the normal gravity vector obtained from the reference ellipsoid model such as WGS-84. The DOV, which is a vector quantity, is usually decomposed into two mutually perpendicular components: a north-south or meridional component ξ, which is reckoned positive northward, and an east-west or prime vertical component η, which is reckoned positive eastward [27]. In other words, the deflection components are positive if the direction of the gravity vector points further south and further west than the corresponding ellipsoidal normal [28], or the level surface is rising to the south or west, respectively, with respect to the ellipsoid [29].
The north component of horizontal gravity disturbance Δ g N n is associated with ξ, and the east component of horizontal gravity disturbance Δ g E n is associated with η. As shown in Figure 2, Equation (1) describes the connections between DOV and horizontal gravity disturbance where γ is the norm of the normal gravity vector γ:
tan ξ ξ = Δ g N n γ tan η η = Δ g E n γ

2.2. Calculation of Horizontal Gravity Disturbance Based on SHM

According to the potential theory [11], DOV is the partial derivative of the disturbed gravitational potential T, and r is the distance from the center of reference ellipsoid to the calculated point:
ξ = 1 γ · r T L η = 1 γ · r · cos L T λ
Substituting Equation (2) into Equation (1), the connections between the disturbed gravitational potential and the horizontal gravity disturbance can be built, and L is the geocentric latitude of calculated point and λ is longitude of the calculated point:
Δ g N n = 1 r T L Δ g E n = 1 r · cos L T λ
Usually the geocentric colatitude ϑ is used in the SHM calculation:
ϑ = π 2 L
Then the horizontal gravity disturbances can be obtained as:
Δ g N n = 1 r T ϑ
Δ g E n = 1 r · sin ϑ T λ
The disturbed gravitational potential T is the solution of Laplace equation in the ellipsoidal coordinate frame, which can be represented as follows [11]:
T = G M r n = 2 n m a x m = 0 n ( a r ) n ( C ¯ n m * cos m λ + S ¯ n m sin m λ ) P ¯ n m ( cos ϑ )
where G is the gravitational constant, M is the mass of the Earth, a is the major semi-axis length of the reference ellipsoid, n and m are called the degree and order of the SHM, C ¯ n m * and S ¯ n m are the coefficients of the SHM, P ¯ n m ( cos ϑ ) is the fully normalized Legendre functions of degree n and order m. The partial derivatives of the disturbed gravitation potential are Equations (8) and (9) [11]:
T ϑ = G M r n = 2 n m a x m = 0 n ( a r ) n ( C ¯ n m * · cos m λ + S ¯ n m · sin m λ ) d P ¯ n m ( cos ϑ ) d ϑ
T λ = G M r n = 2 n m a x m = 0 n ( a r ) n [ m ( C ¯ n m * · sin m λ + S ¯ n m · cos m λ ) ] P ¯ n m ( cos ϑ )
Substituting Equations (8) and (9) into Equations (5) and (6), the formulas of calculating horizontal gravity disturbance are obtained:
Δ g N n = G M r 2 n = 2 n m a x m = 0 n ( a r ) n ( C ¯ n m * · cos m λ + S ¯ n m · sin m λ ) d P ¯ n m ( cos ϑ ) d ϑ
Δ g E n = G M sin ϑ · r 2 n = 2 n m a x m = 0 n ( a r ) n [ m ( C ¯ n m * · sin m λ + S ¯ n m · cos m λ ) ] P ¯ n m ( cos ϑ )

3. The Effect of Accelerometer Bias on Horizontal Gravity Disturbance Compensation

3.1. Reference Coordinate Frames

3.1.1. Earth-Centered-Earth-Fixed Frame e

The origin of this coordinate frame is at center of the Earth, whose z-axis points in the direction of the North pole, x-axis points towards the Greenwich Meridian, and y-axis completes the right-handed orthogonal frame. This frame rotates with the Earth with the rate ω i e e = [ 0   0   Ω ] T , Ω is the Earth’s rotation angular velocity, as shown in Figure 3a.

3.1.2. Navigation Coordinate Frame with North-Up-East Definition n

This frame is a local geodetic north-oriented, local-level coordinate frame, the origin of this frame is at the position of the ship, and its xn−yn−zn axes respectively point towards North-Up-East, as shown in Figure 3a. It should be noted that zn is collinear with the normal gravity vector which points towards the center of the reference ellipsoid.

3.1.3. Body Coordinate Frame with Forward-Upward-Right b

This frame is defined based on the input axes of inertial sensors. Its axes respectively point towards Right-Forward-Upward of the ship.

3.1.4. Navigation Coordinate Frame with True Vertical n’

This frame is similar to n coordinate frame. Based on the true gravity vector, whose axes are denoted by { xn′,yn′,zn′} and xn′ also points towards north, but zn′ is collinear with the true gravity vector and yn′ can be determined based on the right-hand rule, as shown in Figure 3b.
Direction cosine matrix (DCM) is used to express a rotation in three dimensions as a mathematical transformation. The DCM between navigation frame and body frame is a function of Euler angles:
r n = C b n r b
where r n is the projection of arbitrary vector r in n coordinate frame, r b is the projection of vector r in b coordinate frame.

3.2. Mathmatical Formulation of INS

Inertial navigation is an integration algorithm based on Newton’s second law. It can be decomposed into two successive steps as shown in Figure 4. Step I is the initial alignment in which the initial values of the integration algorithm are obtained, including initial attitudes, initial velocities, and initial positions. Step II is the integration calculation named navigation calculation, including attitude calculation, velocity calculation and position calculation. In fact, the implementation of step I contains step II. Navigation calculation is first executed in step I, then Kalman Filter (KF) recursion is performed following the one-step navigation calculation. After the KF measurement update, the estimated state can be fed back to fix the corresponding navigation calculation errors. Finally, the precise initial values will be obtained in the KF recursion.
Mechanization is very important for inertial navigation and the analysis in this paper. The traditional INS strapdown mechanization is chosen in this paper for two reasons. Firstly, as introduced in Section 1, the background of this paper is compensation of horizontal gravity disturbance for high precision INS. Secondly, the accelerometer bias is estimated from the strapdown gravity vector measurement, and in the field of strapdown gravity vector measurement, the traditional INS strapdown mechanization is appropriate [30,31,32].
The navigation calculation is performed in the navigation frame and all the vectors should be transformed into this frame before they can be used. Navigation coordinate frame with north-up-east definition is usually chosen as the coordinate frame in which the navigation calculation is implemented.
The attitude kinematical equation with DCM parameterization is given by [2]:
C ˙ b n = C b n [ ω n b b × ]
where ω n b b is the body angular rate with respect to the navigation frame n and is given by [2]:
ω n b b = ω i b b C n b ( ω i e n + ω e n n )
ω i b b is the body angular rate with respect to inertial frame and is measured by the gyroscopes. ω i e n is the earth rotational rate and is given by [2]:
ω i e n = [ Ω cos L Ω s i n L 0 ] T
ω i e n is the angular rate of navigation frame with respect to the Earth frame, which is caused by the linear motion of the ship on the ellipsoidal surface. The formulation of ω i e n is given by [2]:
ω e n n = [ v E n R N + h v E n R N + h tan L v N n R M + h ] T
R M and R N are the meridian and transverse radius of the ellipsoid curvature, respectively. v E n and v N n are the east and north components of the velocity, respectively. h is the height of the ship relative to the reference ellipsoid, and it should be noted that this paper focuses on the ship mounted INS and the height can be set to zero.
The velocity kinematical equation in navigation frame is given by [2]:
v ˙ e n = C b n f b ( 2 ω i e n + ω e n n ) × v e n + g n
where v e n is the velocity relative to the Earth, f b is the specific force measured by the accelerometers, and it’s emphasized here that g n is the true gravity vector at the position of ship. Δ g n is the horizontal gravity disturbance described in Section 2, and γ n is the normal gravity vector obtained from reference ellipsoid model:
g n = [ Δ g N n γ Δ g E n ] = [ Δ g N n 0 Δ g E n ] + [ 0 γ 0 ] = Δ g n + γ n
The definition of horizontal gravity disturbance compensation can be defined here, that is horizontal gravity disturbance compensation means that the gravity vector used in initial alignment and navigation calculation is g n rather than γ n .
The position kinematical equation is given here and the kinematical equation of height is not considered here for ship mount INS [2]:
L ˙ = 1 R M + h v N n λ ˙ = sec L R N + h v E n

3.3. Effect of Acceleromter Bias on Compensation

Attitude error equation and velocity error equation are the foundations of KF recursion in initial alignment and are also the key point of the analysis. The attitude error equation is given in [2]:
δ φ = [ δ α δ ϕ δ β ] T
δ φ ˙ = ( ω i e n + ω e n n ) × δ φ + δ ω i n n C b n δ ω i b b
where δ φ is the attitude error vector, δ α is the roll error, δ ϕ is the yaw error, δ β is the pitch error and δ ω i b b is the gyroscope noise which can be usually regarded as white noise.
The velocity error equation is derived as follows. The true velocity kinematic equation is given in [2], as Equation (22):
v ˙ e n = C b n f b ( 2 ω i e n + ω e n n ) × v e n + γ n + Δ g n
When the navigation calculation is implemented without the gravity disturbance compensation and no accelerometer bias exists, the practical velocity kinematic equation in navigation calculation is Equation (23):
v ˜ ˙ e n = C ˜ b n f b ( 2 ω ˜ i e n + ω ˜ e n n ) × v ˜ e n + γ n + δ f n
where v ˜ e n is the calculated velocity-containing error. ω ˜ i e n and ω ˜ e n n are angular velocities containing error, γ n is the normal gravity vector, and δ f n is the white noise in accelerometer output. C ˜ b n is the DCM containing error and is defined as follows:
C ˜ b n = ( I 3 × 3 [ δ φ × ] ) C b n
[ δ φ × ] = [ 0 δ β δ ϕ δ β 0 δ α δ ϕ δ α 0 ]
The velocity error equation without accelerometer bias is derived by subtracting Equation (22) from Equation (23):
δ v e n = v ˜ e n v e n = [ δ v N n δ v U n δ v E n ] T
[ f n × ] = [ 0 f E n f U n f E n 0 f N n f U n f N n 0 ]
δ v ˙ e n = [ f n × ] δ φ + δ f n Δ g n
where δ v e n is the velocity error vector, δ v N n , δ v U n and δ v E n are the north, upward and east component of velocity error vector. f n is the specific force measured at the initial point and projected in the n frame. f N n , f U n and f E n are the components of fn in the n frame.
Initial alignment is usually implemented when the INS is at static or mooring, fn can be regarded as the opposite of the true gravity vector:
f n = g n = [ Δ g N n γ Δ g E n ] T
The Kalman filter state is updated with a new observation, velocity error is usually used as observation in static initial alignment. Because the initial alignment is implemented when the ship is at static or moored, the true value of velocity is approximately equal to zero, then the non-zero velocity output of INS is the velocity error. KF recursion of initial alignment will converge when the velocity error reaches zero, and the attitude estimation errors can be obtained by setting Equation (28) equal to zero:
δ α Δ g E n γ = η · γ γ = η
δ β Δ g N n γ = γ · ξ γ = ξ
From Equations (30) and (31), it can be seen that horizontal gravity disturbance will give rise to the attitude errors in the initial alignment. Compensation for horizontal gravity disturbance is necessary in initial alignment.
However, if some accelerometer bias exists, the coupling between accelerometer bias and horizontal gravity disturbance may decrease the compensation effect as analyzed below. When navigation calculation is implemented without gravity disturbance compensation and accelerometer bias exists, the velocity kinematic equation is Equation (32)
v ˜ ˙ e n = C ˜ b n f b ( 2 ω ˜ i e n + ω ˜ e n n ) × v ˜ e n + γ n + δ f n + b a n
The new velocity error equation is obtained by subtracting Equation (22) from Equation (32):
b a n = [ b a , N n b a , U n b a , E n ] T
δ v ˙ e n = [ f n × ] δ φ + δ f n Δ g n + b a n
where b a n is the accelerometer bias and b a , N n , b a , U n and b a , E n are the north, upward and east component of the accelerometer bias. The estimation errors of attitude are obtained as Equations (35) and (36) which clearly describe the connection among attitude error, horizontal gravity disturbance and accelerometer bias:
δ α Δ g E n + b a , E n γ = η · γ γ + b a , E n γ = η + b a , E n γ
δ β Δ g N n + b a , N n γ = γ · ξ γ + b a , N n γ = ξ + b a , N n γ
In the case of the accelerometer bias being much larger than the horizontal gravity disturbance, whether the horizontal gravity disturbance being compensated will not obviously improve the precision of initial alignment, because the accelerometer bias is the main error source. The accelerometer bias can be estimated through INS/GNSS integrated Kalman filter [21].
In the case of the accelerometer bias being much smaller than the horizontal gravity disturbance, the compensation of horizontal gravity disturbance which is the dominant error source will reliably improve the accuracy of initial alignment, and there is no need to estimate accelerometer bias in practice.
In the case of accelerometer bias being at the same order with horizontal gravity disturbance, this is the hardest situation to handle. Taking Equation (35) for example, if the sign and magnitude of east accelerometer bias are the same as those of east gravity disturbance, the effect of horizontal gravity disturbance on initial alignment is counteracted with the accelerometer bias. If the gravity disturbance is still compensated in initial alignment, then a new attitude error will arise due to the compensation. Therefore, estimation of accelerometer bias is necessary in this situation, otherwise whether do the gravity disturbance compensation in initial alignment will be ambiguous as mentioned in the introduction.
The above analysis also shows the effect of horizontal gravity disturbance on navigation calculation. As in Equation (32), the accelerometer bias is coupled with the horizontal gravity disturbance in the velocity calculation. When the horizontal gravity disturbance is compensated, the practical kinematic equation is Equation (37):
v ˜ ˙ e n = C ˜ b n f b ( 2 ω ˜ i e n + ω ˜ e n n ) × v ˜ e n + γ n + Δ g n + δ f n + b a n
Comparing Equation (37) with Equation (22), it can be found that the effect of compensation is counteracted with the accelerometer bias. In the limit situation of the accelerometer bias being the opposite of horizontal gravity disturbance, the compensation will be inefficient.

4. Model of Gravity Vector Measurement Noise

From the analysis in Section 3, when the accelerometer bias is at the same order with horizontal gravity disturbance, accurate estimation of accelerometer bias is the crucial issue in gravity disturbance compensation. Usually the accelerometer bias can be estimated using INS/GNSS integrated Kalman filter in which the accelerometer bias is modeled as a constant. The disadvantage to this approach is that the estimation is a combination of the accelerometer bias and horizontal gravity disturbance. Maybe we can try to do the estimation from another angle. Inertial sensors not only can be used to navigate but also can be used to measure the gravity vector. If the gravity information is the input, position information is obtained. On the contrary, gravity vector can be measured as position and velocity information being the input, that’s the principle of the strapdown gravimeter [25].
Horizontal gravity disturbance can be measured by subtracting the normal gravity vector from the measured value of gravity vector as Equation (38):
Δ g n + w = v ˙ e n C b n f b + ( 2 ω i e n + ω e n n ) × v e n γ n
Equation (38) is the transformation of Equation (22). In gravity vector measurement, some terms of Equation (38) is calculated based on GNSS information, and some are provided by inertial sensors. v ˙ e n is the difference of velocity provided by GNSS, ω i e n and ω e n n can be calculated by substituting the velocity and position provided by GNSS into Equation (15) and Equation (16). C b n is calculated with the output of gyroscopes and Equation (13). f b is the output of the accelerometers.
The precise measurement of specific force is crucial for gravity vector measurement and INS is aided with GNSS to improve attitude accuracy, thus high precise measurement of specific force can be obtained. As described in [33], there are four main data fusion schemes for INS/GNSS. In the traditional strapdown gravity vector measurement, feedback correction is the main data fusion scheme of INS/GNSS which is also the data fusion scheme adopted in this paper. And a new and compact data fusion scheme proposed in [33] is a novel and worthwhile approach in the advancement of strapdown gravity vector measurement.
In addition, it should be noted that w is the measurement noise of gravity vector, which is a composite term due to some error sources except accelerometer bias, including noise of inertial sensors and error of GNSS. The measured value of horizontal gravity disturbance can be expressed as follows:
Δ g N , E n = [ Δ g N n Δ g E n ] T w = [ w N w E ] T
Δ g ˜ N , E n = Δ g N , E n + w
where Δ g N , E n is the true value of horizontal gravity disturbance, Δ g ˜ N , E n is the measured value containing error, w is the measurement noise. When the accelerometer bias exits, the accelerometer bias is added to the measurement equation and the measured value also changes:
Δ g n + w + Δ f n = v ˙ e n C b n f b + ( 2 ω i e n + ω e n n ) × v e n γ n
Δ g ˜ N , E n = Δ g N , E n + w + [ b a , N n b a , E n ] T
The measurement error of horizontal gravity disturbance can be obtained by subtracting the measured value from the true value:
δ g N , E n = [ δ g N n δ g E n ] = [ b a , N n b a , E n ] + [ w N w E ]
where δ g N n and δ g E n are the north and east component of measurement error of horizontal gravity disturbance. Because the background of this paper is using the calculated horizontal gravity disturbance based on SHM to compensate INS, the horizontal gravity disturbance from SHM is regarded as the true value, the measurement error can be calculated by subtracting the true value from the measured value.
In fact, we can have estimations of accelerometer biases by averaging the measurement error as follows:
b ^ a , N n = 1 N i = 1 N [ δ g N n ( i ) ] b ^ a , E n = 1 N i = 1 N [ δ g E n ( i ) ]
b ^ a , N n and b ^ a , E n will be precise estimations if and only if the measurement noise is zero-mean, but it is a too rigorous condition for practice. Based on the Equation (44), when the mean value of measurement noise isn’t zero, is it possible to have a precise estimation of the accelerometer bias? The conjecture is that accelerometer bias is the inherent error of accelerometer which is not related to the Earth’s gravity field, but the contribution of measurement noise to measurement error may be associated with the gravity field.

4.1. The Measurement Noise of Horizontal Gravity Disturbance

In this subsection, the model of measurement noise will be derived. The true gravity vector projected in n coordinate frame is Equation (18), while the true gravity vector projected in n coordinate frame is Equation (45):
g n = [ 0 γ 0 ] T
It should be noted that the y-axis of n coordinate frame is collinear with the true gravity vector, the horizontal components of g n are zeroes. The DCM between n coordinate frame and n coordinate frame can be determined based on the geometric relationship between the two navigation coordinate frames, as shown in Figure 5.
The coordinate frame n′ can be obtained through rotating the coordinate frame n by ϑ along the rotational axis u. It is obvious that the rotational axis and rotational angle are associated with the DOV components and can be determined based on some constraints. The rotation axis u n = [ u x u y u z ] T satisfies four constraints:
(1)
In the plane x n z n ;
(2)
Pass through the origin of the two coordinate frames;
(3)
Be orthogonal to the plane y n y n ;
(4)
u is unit vector;
Based on constraint 1:
u y = 0
Based on constraint 2, k is a constant to be determined:
u x = k · u z
Based on constraint 3:
u n · g n = [ u x u y u z ] T · [ ε · γ γ η · γ ] T = 0
Based on constraint 4:
u x 2 + u y 2 + u z 2 = 1
According to Equations (46)–(49), the rotation axis u can be determined:
u x = ( η ) / ξ 2 + η 2 u y = 0 u z = ξ / ξ 2 + η 2
u = [ u x 0 u z ] = [ η / ξ 2 + η 2 0 ξ / ξ 2 + η 2 ] T
ϑ is the angle between γ n and g n , and can be calculated based on vector product:
γ n = [ 0 γ 0 ] T g n = [ δ g N γ δ g E ] cos ϑ = ( g ¯ · g ) / ( g ¯ · g )
ϑ = arccos ( 1 / 1 + ξ 2 + η 2 )
Based on the rotation axis and the rotation angle, the corresponding quaternion and DCM can be obtained [2]:
Q = [ q 0 q 1 q 2 q 3 ] = [ cos ( ϑ / 2 ) u x · cos ( ϑ / 2 ) u y · cos ( ϑ / 2 ) u z · cos ( ϑ / 2 ) ]
C n n = [ q 0 2 + q 1 2 q 2 2 q 3 2 2 ( q 1 · q 2 q 0 · q 3 ) 2 ( q 1 · q 3 + q 0 · q 2 ) 2 ( q 1 · q 2 + q 0 · q 3 ) q 0 2 q 1 2 + q 2 2 q 3 2 2 ( q 2 · q 3 q 0 · q 1 ) 2 ( q 1 · q 3 q 0 · q 2 ) 2 ( q 2 · q 3 + q 0 · q 1 ) q 0 2 q 1 2 q 2 2 + q 3 2 ]
Equations (54) and (55) can be simplified due to uy being zero:
Q = [ cos ( ϑ / 2 ) u x · cos ( ϑ / 2 ) 0 u z · cos ( ϑ / 2 ) ]
C n n = [ q 0 2 + q 1 2 q 3 2 2 ( q 0 · q 3 ) 2 ( q 1 · q 3 ) 2 ( q 0 · q 3 ) q 0 2 q 1 2 q 3 2 2 ( q 0 · q 1 ) 2 ( q 1 · q 3 ) 2 ( q 0 · q 1 ) q 0 2 q 1 2 + q 3 2 ]
Based on the Equations (50)–(57), it can be seen that C n n is a function of the horizontal gravity disturbance:
C n n = F ( Δ g N n , Δ g E n )
and Equation (18) can be approximately rearranged as:
C n n g n = C n n [ 0 γ 0 ] = g n = [ Δ g N n γ Δ g E n ]
Now, when only considering the measurement noise, the DCM based on the measured values is supposed to has the following form:
F ( Δ g N n + w N , Δ g E n + w E ) = C ˜ n n = ( C n n + δ C n n )
Substituting Equation (60) into Equation (59):
[ Δ g N n + w N γ Δ g E n + w E ] C ˜ n n [ 0 γ 0 ] = ( C n n + δ C n n ) [ 0 γ 0 ] = C n n [ 0 γ 0 ] + δ C n n [ 0 γ 0 ]
The expression of measurement noise is obtained:
[ w N 0 w E ] = δ C n n [ 0 γ 0 ]
Suppose δ C n n has the following form:
δ C n n = [ δ c 11 δ c 12 δ c 13 δ c 21 δ c 22 δ c 23 δ c 31 δ c 32 δ c 33 ]
The model of measurement noise is obtained, this is the foundation of estimating accelerometer bias in this paper:
w N = δ c 12 · ( γ ) w E = δ c 32 · ( γ )
Substituting Equation (64) into Equation (43), the model of measurement error is rearranged as:
[ δ g N n δ g E n ] = [ b a , N n δ c 12 · ( γ ) b a , E n δ c 32 · ( γ ) ]

4.2. Parameters of the Model of Measurement Noise

δ c 12 and δ c 32 are derived below. According Equations (50)–(57), the elements of the quaternion are also the functions of measured values of horizontal gravity disturbance. These elements will also contain some error due to the measurement noise as follows:
q ˜ 0 = q 0 + δ q 0 q ˜ 1 = q 1 + δ q 1 q ˜ 3 = q 3 + δ q 3
δ q 0 , δ q 1 and δ q 3 are the errors caused only by the measurement noise. Substituting Equation (66) into Equation (57), we can get:
δ c 12 = 2 ( q 0 δ q 3 + δ q 0 q 3 ) δ c 32 = 2 ( q 0 δ q 1 + δ q 0 q 1 )
δ C n n = [ 2 ( q 0 · δ q 0 + q 1 · δ q 1 q 3 · δ q 3 ) 2 ( q 0 δ q 3 + δ q 0 q 3 ) 2 ( q 1 δ q 3 + δ q 1 q 3 ) 2 ( q 0 δ q 3 + δ q 0 q 3 ) 2 ( q 0 · δ q 0 q 1 · δ q 1 q 3 · δ q 3 ) 2 ( q 0 δ q 1 + δ q 0 q 1 ) 2 ( q 1 δ q 3 + δ q 1 q 3 ) 2 ( q 0 δ q 1 + δ q 0 q 1 ) 2 ( q 0 · δ q 0 q 1 · δ q 1 + q 3 · δ q 3 ) ]
The model of measurement error is rearranged as:
[ δ g N n δ g E n ] = [ b a , N n + [ 2 ( q 0 δ q 3 + δ q 0 q 3 ) ] · ( γ ) b a , E n [ 2 ( q 0 δ q 1 + δ q 0 q 1 ) ] · ( γ ) ]
The next work is to derivate the expression of δ q i , i = 0 , 1 , 3 through the following partial differential equations:
δ q 0 = ( q 0 ξ ) δ ξ + ( q 0 η ) δ η δ q 1 = ( q 1 ξ ) δ ξ + ( q 1 η ) δ η δ q 3 = ( q 3 ξ ) δ ξ + ( q 3 η ) δ η
According to the Equations (50)–(57), these partial derivatives are obtained:
f ( ξ , η ) = 1 / 1 + ξ 2 + η 2 g ( f ) = arccos ( f )
q 0 ξ = ( sin ( g ( f ) / 2 ) · ξ 2 ( 1 + η 2 + ξ 2 ) η 2 + ξ 2 )
q 0 η = sin ( g ( f ) / 2 ) · η 2 ( 1 + η 2 + ξ 2 ) η 2 + ξ 2
q 1 ξ = ( ξ · η ) ( η 2 + ξ 2 ) [ cos ( g ( f ) / 2 ) 2 ( 1 + η 2 + ξ 2 ) + sin ( g ( f ) / 2 ) ξ 2 + η 2 ]
q 1 η = cos ( g ( f ) / 2 ) · η 2 2 ( 1 + ξ 2 + η 2 ) ( ξ 2 + η 2 ) sin ( g ( f ) / 2 ) · ( ξ 2 ) ( ξ 2 + η 2 ) 3 / 2
q 3 ξ = cos ( g ( f ) / 2 ) · ξ 2 2 ( 1 + η 2 + ξ 2 ) · ( η 2 + ξ 2 ) + sin ( g ( f ) ) · ( η 2 ) ( ξ 2 + η 2 ) 3 / 2
q 3 η = cos ( g ( f ) / 2 ) ξ · η 2 ( 1 + ξ 2 + η 2 ) ( ξ 2 + η 2 ) sin ( g ( f ) / 2 ) · ( ξ · η ) ( ξ 2 + η 2 ) 3 / 2

4.3. Estimation Model of the Accelerate Bias

Finally, based on the model of measurement noise, the model of measurement error is rearranged as follows:
[ δ g N n δ g E n ] = [ 2 γ · Φ 1 2 γ · Φ 2 1 0 2 γ · Φ 3 2 γ · Φ 4 0 1 ] [ δ ξ δ η b a , N n b a , E n ]
Φ 1 = { cos ( g ( f ) 2 ) ( q 3 ξ ) + sin ( g ( f ) 2 ) · ( ξ ξ 2 + η 2 ) ( q 0 ξ ) }
Φ 2 = { cos ( g ( f ) 2 ) ( q 3 η ) + sin ( g ( f ) 2 ) · ( ξ ξ 2 + η 2 ) ( q 0 η ) }
Φ 3 = { cos ( g ( f ) 2 ) ( q 1 ξ ) + sin ( g ( f ) 2 ) · ( η ξ 2 + η 2 ) · ( q 0 ξ ) }
Φ 4 = { cos ( g ( f ) 2 ) ( q 1 η ) + sin ( g ( f ) 2 ) · ( η ξ 2 + η 2 ) · ( q 0 η ) }
δ ξ and δ η can be regarded as the measurement errors of DOV only due to the measurement noise w N and w E , and the contribution of δ ξ and δ η to the measurement error of horizontal gravity disturbance is associated with the variation of gravity field, these connections are represented by Φ i , i = 1 , 2 , 3 , 4 .

5. Simulation Results

5.1. Estimation Methods of Accelerometer Bias

The model of measurement error, as Equation (78), can be utilized in the method of least squares to estimate the accelerometer bias, the elements Φ i , i = 1 , 2 , 3 , 4 change with the variation of gravity field, then uncorrelated observation equations are built to ensure the numerical stability of solving the pseudo-inverse in the least squares. Assuming that N times measurements of horizontal gravity disturbance are made, then there are 2N observation equations built as below:
z 2 N × 1 = H 2 N × 4 · x
x = [ x 1 x 2 x 3 x 4 ] T = [ δ ξ δ η b a , N n b a , E n ] T
z 2 N × 1 = [ δ g N n ( 1 ) δ g E n ( 1 ) δ g N n ( N ) δ g E n ( N ) ]
H 2 N × 4 = [ 2 γ ( 1 ) · Φ 1 ( 1 ) 2 γ ( 1 ) · Φ 2 ( 1 ) 1 0 2 γ ( 1 ) · Φ 3 ( 1 ) 2 γ ( 1 ) · Φ 4 ( 1 ) 0 1 2 γ ( N ) · Φ 1 ( N ) 2 γ ( N ) · Φ 2 ( N ) 1 0 2 γ ( N ) · Φ 3 ( N ) 2 γ ( N ) · Φ 4 ( N ) 0 1 ]
where z 2 N × 1 is the measurement error of horizontal gravity disturbance which is calculated by subtracting the true values from the measurement values, and the subscript means the vector has 2 N rows and one column. H 2 N × 4 is the observation matrix which has 2N rows and four columns, Φ i , i = 1 , 2 , 3 , 4 can be determined by substituting true values of DOV into Equations (79)–(82). γ ( i ) is the norm of the normal gravity vector.
The accelerometer bias can be estimated by least squares as follows. b ^ a , N n and b ^ a , E n are the estimations of the accelerometer biases:
x = ( H 2 N × 4 T · H 2 N × 4 ) 1 · H 2 N × 4 T · z 2 N × 1
b ^ a , N n = x 3 b ^ a , E n = x 4

5.2. Simulation of Horizontal Gravity Disturbance

From the model of measurement error, it can be known that the contribution of measurement noise to measurement error changes with the variation of gravity field. In this simulation, an area with moderate variation of gravity field is chosen whose latitude range is from 1° N to 5° N and longitude range is from 76° E to 80° E, as shown in Table 1. The horizontal gravity disturbance in this area is calculated based on SHM as described in Section 2, the SHM adopted here is EGM2008 whose the maximum degree is 2190. The horizontal gravity disturbance in the simulation area is calculated with maximum degree and the grid spacing is 5 n mile, the horizontal gravity disturbance in this simulation area is shown in Figure 6.
The ship is assumed to sail along some straight lines in this area, these straight lines are called survey lines. The gravity disturbance measurement is implemented along these survey lines. Five latitude lines and five longitude lines in this area are taken as survey lines as shown in Figure 7 and the horizontal gravity disturbance on the survey lines is shown in Figure 8.

5.3. Simulation Results

In this subsection, the efficiency of the proposed estimation method will be verified by the simulations. The true values of horizontal gravity disturbances are the calculation results of SHM, and the measured values are constructed as Equation (42), then the measurement error can be obtained by subtracting the true values from the measurement values.
The bias of accelerometer b a is sometimes split into static component and dynamic component [34]. The static component denoted by b a s is also known as the bias repeatability and comprises the run-to-run variation of instrument bias plus the residual fixed bias remaining after calibration. b a s is constant throughout an IMU operating period but varies from run to run. The dynamic component denoted by bad is also known as the bias instability and varies over time and temperature:
b a = b a s + b a d
The value of bias repeatability is chosen with reference to the datasheet of QA3000, which is the highest navigation-grade accelerometer from Honeywell®. The one-year composite repeatability of QA3000 is less than 40 mGal, this value can be regarded as the upper bound of accelerometer bias. According to the analysis in the end of Section 3, when the accelerometer bias is at the same order with horizontal gravity disturbance, estimation of accelerometer bias is a crucial issue in the horizontal gravity disturbance compensation, so the magnitude of accelerometer bias is carefully chosen based on the magnitude of horizontal gravity disturbance in the simulation area. The statistical values of the horizontal gravity disturbance on the survey lines are list in Table 2.
The values of bias instability are chosen with reference to the measured data of a high-quality accelerometer. A long term static experiment for investigating the bias instability characteristics of this accelerometer is carried out in December 2016. The output of the accelerometer is recorded for 450 h and the accelerometer is mounted on a stable platform. The sampling frequency is 200 Hz and the mean values of outputs are calculated every one hundred seconds as shown in Figure 9a. The drift of the accelerometer output can be obtained by subtracting the first mean value form the subsequent mean values as shown in Figure 9b.
It can be seen from Figure 9 that the characteristics of bias instability changes over time. From the beginning to the 150th hour, the bias instability of accelerometer is a linear drift, and the drift rate is approximately 1.31 mGal/day. From the 150th hour to the 250th hour, the bias instability changes to a second-order curve. From the 250th hour to the end, the bias instability changes to a linear drift again, and the drift rate is approximately 0.59 mGal/day.
Maintaining a nearly constant and stable temperature of the IMU improves its accuracy during calibration and operation, as temperature stability is directly related to the sensor accuracy. A multi-point thermal control is used in this experiment to improve the precision of the accelerometer, and the temperatures of the accelerometer and the environment are recorded for the first 140 h. The connections between the outputs of the accelerometer and the measured temperatures are shown in Figure 10 and Figure 11. Comparing Figure 10 and Figure 11, there is no cyclic change in the temperature of the accelerometer while the temperature of the environment changing day and night, which means under the control the temperature of the accelerometer is stable. And there is no obvious correlation between the output of the accelerometer and the temperature of the accelerometer. This indicates that the bias instability of accelerometer is approximately irrelevant to the temperature in this experiment, and can be considered to be a linear drift over time, then the model of bias instability can be built as Equation (90), where t is time, k a d , N n is the drift rate of north component and k a d , E n is the drift rate of the east component:
b a d , N n = k a d , N n · t b a d , E n = k a d , E n · t
The bias repeatability is supposed to be equal to the mean value or median of all the survey lines, and the bias instability is supposed to be equal to the first drift rate in the above static experiment, 1.31 mGal/day. The measurement noise ranges from 1 mGal to 20 mGal, which is a reasonable range for high precision inertial sensors and GNSS receiver. The signal-to-noise ratio (SNR) is defined as:
S N R = 20 log 10 ( m a g n i t u d e   o f   a c c e l e r o m e t e r   b i a s m a g n i t u d e   o f   m e a s u r e m e n t   n o i s e ) d B
Three simulations are designed here to verify the validity of the proposed method as described in Table 3, and the differences among these simulations are the types and scales of the accelerometer bias.
The estimation errors of accelerometer bias in Simulation I are shown in Figure 12 and the estimation errors of accelerometer bias in Simulation II are shown in Figure 13. From Figure 12 and Figure 13, the precise estimations of accelerometer bias are obtained in the proposed method. Even if the SNR is −15 dB, the estimation error is less than 10 2 %. For estimating the north component of accelerometer bias, the estimation error of survey line 10 is significantly greater than the errors of other survey lines, this may be explained from Table 2, the mean value and median of survey line 10 are both large values which indicates the variation of north component of horizontal gravity disturbance on line 10 is the most drastic. The more drastic variation of gravity field, the more contribution of measurement noise to measurement error, so the weight of accelerometer bias in measurement noise is less. The decrease of the effective information in the observed quantity leads to the decrease of the estimation accuracy. This interpretation is also suitable for survey line 1, on which the variation of east component of horizontal gravity disturbance is the most drastic, and the estimation error of east component of accelerometer bias is significantly greater.
In Simulation III, some modifications are needed on Equations (83)–(86) to estimate the bias repeatability and drift rate simultaneously:
z 2 N × 1 = H 2 N × 6 · x
x = [ δ ξ δ η b a s , N n b a s , E n k a d , N n k a d , E n ] T
H 2 N × 6 = [ 2 γ ( 1 ) · Φ 1 ( 1 ) 2 γ ( 1 ) · Φ 2 ( 1 ) 1 0 t 0 2 γ ( 1 ) · Φ 3 ( 1 ) 2 γ ( 1 ) · Φ 4 ( 1 ) 0 1 0 t 2 γ ( N ) · Φ 1 ( N ) 2 γ ( N ) · Φ 2 ( N ) 1 0 t 0 2 γ ( N ) · Φ 3 ( N ) 2 γ ( N ) · Φ 4 ( N ) 0 1 0 t ]
x = ( H 2 N × 6 T · H 2 N × 6 ) 1 · H 2 N × 6 T · z 2 N × 1
b ^ a s , N n = x 3 b ^ a s , E n = x 4 k ^ a d , N n = x 5 k ^ a d , E n = x 6
where t is time, b ^ a s , N n and b ^ a s , E n are the estimations of bias repeatability, k ^ a d , N and k ^ a d , E are the estimations of drift rates.
In order to obtain more observations to improve the accuracy of estimation, the speed of the ship is assumed to be 5 knot, and it can be seen from Table 2 that the sailing time is 48 h. The estimation error in Simulation III are shown in Figure 14 and Figure 15.
In Simulation III, the accelerometer bias is time varying, so the x-axis is measurement noise instead of SNR. And from Figure 14 and Figure 15, the estimation errors of bias repeatability are less than 10−2%, and the estimation errors of drift rate are also less than 10−2%. The results of simulation III indicate that the proposed method can accurately estimate the bias repeatability and bias instability simultaneously.
It should be noted that, in the above simulations, although a precise estimation of accelerometer bias is obtained, the length of the survey line is 240 n miles, which is too long and not suitable for practical applications. Further simulations are carried out for the typical application scenario.
The ship’s speed is generally 20 knots, and the time for estimating accelerometer bias is assumed to be one hour. The sample period of strapdown gravimeter is generally 160~300 s [25], then the distance between the measurement points is supposed to be 1 n mile. Based on the three assumed conditions, the length of survey line is supposed to be 20 n miles and the coordinate range of these short survey lines are list in Table 4 and the horizontal gravity disturbances on the short survey lines are shown in Figure 16.
Three simulations are designed here to verify the validity of the proposed method on the short survey lines as described in Table 5, and the differences among these simulations are also the types and scales of the accelerometer bias. The measurement noise is consistent with the previous setting.
The estimation errors of accelerometer bias on the short survey lines in Simulation IV are shown in Figure 17 and the estimation errors of accelerometer bias on the short survey lines in Simulation V are shown in Figure 18. From Figure 17 and Figure 18, it can be seen that the maximum estimation error of accelerometer bias is less than 10 1 % in the typical application scenario. It can be known from Table 4 that there are only twenty-one measured values of the horizontal gravity disturbance in this application scenario. Compared the estimation errors with those of the long survey lines, the reduction of estimation accuracy is due to the reduction of observation.
The estimation errors of bias repeatability in Simulation VI are shown in Figure 19 and the estimation errors of bias instability in Simulation VI are shown in Figure 20. It can be found that the estimation errors of bias repeatability are less than 10−1%, and the estimation errors of drift rate are also less than 10−1%. The results in Simulation VI indicate that the proposed method can accurately estimate the bias repeatability and bias instability simultaneously in the typical application scenario. Also due to the reduction of the observations, the estimation errors of bias repeatability and bias instability in simulation VI are both larger than those in Simulation III.
It should be noted that in the above simulations the bias of the gyroscope is assumed to be zero. In practice, gyroscope bias may affect the estimation accuracy of the proposed method, as gyroscope bias will produce some specific force measurement error. In other words, the estimation may be a combination of accelerometer bias and gyroscope bias. As described in [35], observability is crucial for estimation of gyroscope bias in INS/GNSS. From the observability analysis, biases of the horizontal gyroscopes can be accurately estimated. But it’s difficult to accurately estimate the bias of the vertical gyroscope due to the poor observability. Thus, the estimation accuracy of the proposed approach may be affected by the bias of the vertical gyroscope.

6. Conclusions

Horizontal gravity disturbance produces a restriction on the achievable accuracy of INS. Although based on SHM horizontal gravity disturbances can be accurately obtained and compensated, the effect of compensation will decrease due to the accelerometer bias. Through the analysis, it is found that estimating the accelerometer bias is a necessary condition to ensure the effect of horizontal gravity disturbance compensation, especially when the accelerometer bias is of the same order as the horizontal gravity disturbance. The method of estimating the accelerometer bias from the gravity vector measurement is proposed, and the model of measurement noise is derived to separate the accelerometer bias from the measurement error. Simulation results verify the effect of the proposed method, and a precise estimation of accelerometer bias is obtained in the typical application scenario.

Acknowledgments

This work is supported by the National Key R&D Program of China, contract No. 2016YFC0303002; and Youth Funds of Ministry of Land and Resources, the Key Laboratory of Aeronautical Geophysics and Remote Sensing Geology, contract No. 2016YPL06.

Author Contributions

J.T. and J.C. conceived the theory and designed the simulation experiments; M.W. and L.C. did the work of literature search; S.C. analyzed the data and J.L. contributed the data analysis.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Britting, K.R. Inertial Navigation Systems Analysis; Artech House: Norwood, MA, USA, 2010. [Google Scholar]
  2. Titterton, D.; Weston, J.L. Strapdown Inertial Navigation Technology; Institution of Engineering and Technology: Stevenage, UK, 2004. [Google Scholar]
  3. Kwon, J.H.; Jekeli, C. Gravity requirements for compensation of ultra-precise inertial navigation. J. Navig. 2005, 58, 479–492. [Google Scholar] [CrossRef]
  4. Harriman, D.W.; Van Dam, C.C.P.G. Gravity-induced errors in airborne inertial navigation. J. Guid. Control Dyn. 1986, 9, 419–426. [Google Scholar]
  5. Jing, W.; Gongliu, Y.; Xiangyun, L.; Xiao, Z. Application of the spherical harmonic gravity model in high precision inertial navigation systems. Meas. Sci. Technol. 2016, 27, 95103. [Google Scholar] [CrossRef]
  6. Xiao, Z.; Gongliu, Y.; Jing, W.; Jing, L. An improved gravity compensation method for high-precision free-INS based on MEC–BP–AdaBoost. Meas. Sci. Technol. 2016, 27, 125007. [Google Scholar] [CrossRef]
  7. Dai, D.; Wang, X.; Zhan, D.; Huang, Z.; Xiong, H. An improved method for gravity disturbances compensation in INS/GPS integrated navigation. In Proceedings of the 12th International Conference on Signal Processing (ICSP), Hangzhou, China, 19–23 October 2014; pp. 148–153. [Google Scholar]
  8. Cong, L.; Zhao, Z.; Yang, X. On gravity disturbance compensation technology of high-precision SINS based on B-spline method. In Proceedings of the 2014 IEEE Chinese Guidance, Navigation and Control Conference, Yantai, China, 8–10 August 2014; pp. 16–18. [Google Scholar]
  9. Hao, W.; Xuan, X.; Zhi-Hong, D.; Meng-Yin, F. The influence of gravity disturbance on high-precision long-time INS and its compensation method. In Proceedings of the Fourth International Conference on Instrumentation and Measurement, Computer, Communication and Control, Harbin, China, 18–20 September 2014; pp. 104–108. [Google Scholar]
  10. Fang, J.; Chen, L.; Yao, J. An accurate gravity compensation method for high-precision airborne POS. IEEE Trans. Geosci. Remote Sens. 2014, 52, 4564–4573. [Google Scholar] [CrossRef]
  11. Hofmann-Wellenhof, B.; Moritz, H. Physical Geodesy; Springer Vienna: Vienna, Austria, 2009. [Google Scholar]
  12. Kayton, M. Fundamental limitations on inertial measurements. In Guidance and Control; Farrior, J.S., Ed.; Academic Press: Cambridge, MA, USA, 1962; Volume 8, pp. 367–394. [Google Scholar]
  13. Gelb, A.; Levine, S.A. Effect of deflections of the vertical on the performance of a terrestrial inertial navigation system. J. Spacecr. Rockets 1969, 6, 978–984. [Google Scholar] [CrossRef]
  14. Jordan, S.K. Self-consistent statistical models for the gravity anomaly, vertical deflections, and undulation of the geoid. J. Geophys. Res. 1972, 77, 3660–3670. [Google Scholar] [CrossRef]
  15. 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, Palm Springs, CA, USA, 20–23 April 1998; pp. 324–331. [Google Scholar]
  16. Wang, J.; Yang, G.; Li, J.; Zhou, X. An online gravity modeling method applied for high precision free-INS. Sensors 2016, 16, 1541. [Google Scholar] [CrossRef] [PubMed]
  17. Wu, R.; Wu, Q.; Han, F.; Liu, T.; Hu, P.; Li, H. Gravity compensation using EGM2008 for high-precision long-term inertial navigation systems. Sensors 2016, 16, 2177. [Google Scholar] [CrossRef] [PubMed]
  18. Tie, J.; Wu, M.; Cao, J.; Lian, J.; Cai, S. The impact of initial alignment on compensation for deflection of vertical in inertial navigation. In Proceedings of the IEEE 8th International Conference on Cybernetics and Intelligent Systems, Robotics Automation and Mechatronics, Ningbo, China, 19–21 November 2017. [Google Scholar]
  19. Gilster, G. High accuracy performance capabilities of the military standard ring laser gyro inertial navigation unit. In Proceedings of the Position Location and Navigation Symposium, Las Vegas, NV, USA, 11–15 April 1994; pp. 464–473. [Google Scholar]
  20. Hanson, P.O. Correction for deflections of the vertical at the runup site. In Proceedings of the IEEE Position Location and Navigation Symposium, Orlando, FL, USA, USA, 29 November–2 December 1988; pp. 288–296. [Google Scholar]
  21. Huang, Y.; Zhang, K.; Wu, M.; Ma, H. Estimation of accelerometer bias of airborne gravity while aircraft sliding. In Proceedings of the International Conference on Measuring Technology and Mechatronics Automation, Changsha, China, 13–14 March 2010; pp. 153–157. [Google Scholar]
  22. Sinpyo, H.; Man Hyung, L.; Ho-Hwan, C.; Sun-Hong, K.; Speyer, J.L. Observability of error states in GPS/INS integration. IEEE Trans. Veh. Technol. 2005, 54, 731–743. [Google Scholar]
  23. Goshen-Meskin, D.; Bar-Itzhack, I.Y. Observability analysis of piece-wise constant systems with application to inertial navigation. In Proceedings of the 29th IEEE Conference on Decision and Control, Honolulu, HI, USA, 5–7 December 1990; Volume 822, pp. 821–826. [Google Scholar]
  24. Shaw, L.; Paul, I.; Henrikson, P. Statistical models for the vertical deflection from gravity anomaly models. J. Geophys. Res. 1969, 74, 4259–4265. [Google Scholar] [CrossRef]
  25. Cai, S.K.; Tie, J.B.; Zhang, K.D.; Cao, J.L.; Wu, M.P. Marine gravimetry using the strapdown gravimeter SGA-WZ. Mar. Geophys. Res. 2017, 38, 325–340. [Google Scholar] [CrossRef]
  26. Cai, S.K.; Zhang, K.D.; Wu, M.P. Improving airborne strapdown vector gravimetry using stabilized horizontal components. J. Appl. Geophys. 2013, 98, 79–89. [Google Scholar] [CrossRef]
  27. Featherstone, W. The Use and Abuse of Vertical Deflections. In Proceedings of the 40th Australian Surveyors Congress, Fremantle, Australia, 1–6 November 1999; pp. 85–97. [Google Scholar]
  28. Vanícek, P.K.; Krakiwsky, E.J. Geodesy: The Concepts; Elsevier Science: Amsterdam, The Netherlands, 1986. [Google Scholar]
  29. Bomford, G. Geodesy, 4th ed.; Clarendon Press: Oxford, UK, 1980. [Google Scholar]
  30. Cai, S.; Zhang, K.; Wu, M.; Huang, Y.; Yang, Y. An iterative method for the accurate determination of airborne gravity horizontal components using strapdown inertial navigation system/global navigation satellite system. Geophysics 2015, 80, G119–G129. [Google Scholar] [CrossRef]
  31. Cai, S.; Wu, M.; Zhang, K.; Cao, J.; Tuo, Z.; Huang, Y. The first airborne scalar gravimetry system based on SINS/DGPS in China. China Earth Sci. 2013, 56, 2198–2208. [Google Scholar] [CrossRef]
  32. Bruton, A.M. Improving the Accuracy and Resolution of SINS/DGPS Airborne Gravimetry; University of Calgary: Calgary, AB, Canada, 2000. [Google Scholar]
  33. Wagner, J.F.; Wieneke, T. Integrating Satellite and Inertial Navigation–Conventional and New Fusion Approaches. Control Eng. Pract. 2003, 11, 543–550. [Google Scholar] [CrossRef]
  34. Groves, P.D. Principles of GNSS, Inertial, and Multisensor Integrated Navigation Systems, 2nd ed.; Artech House: Norwood, MA, USA, 2013. [Google Scholar]
  35. Wagner, J. GNSS/INS Integration: Still an Attractive Candidate for Automatic Landing Systems. GPS Solut. 2005, 9, 179–193. [Google Scholar] [CrossRef]
Figure 1. The definition of gravity disturbance vector, gravity disturbance and deflection of vertical, and this figure is adapted from [11].
Figure 1. The definition of gravity disturbance vector, gravity disturbance and deflection of vertical, and this figure is adapted from [11].
Sensors 18 00883 g001
Figure 2. The definition of deflection of vertical and horizontal gravity disturbance.
Figure 2. The definition of deflection of vertical and horizontal gravity disturbance.
Sensors 18 00883 g002
Figure 3. The definition of navigation coordinate frame, (a) is the diagram of n and (b) is the diagram of n′.
Figure 3. The definition of navigation coordinate frame, (a) is the diagram of n and (b) is the diagram of n′.
Sensors 18 00883 g003
Figure 4. The schematic of inertial navigation.
Figure 4. The schematic of inertial navigation.
Sensors 18 00883 g004
Figure 5. The geometric relationship between the two navigation coordinate frames.
Figure 5. The geometric relationship between the two navigation coordinate frames.
Sensors 18 00883 g005
Figure 6. The horizontal gravity disturbance in the simulation area. (a) is the north component of the horizontal gravity disturbance and (b) is the east component of the horizontal gravity disturbance. (1 mGal = 10−5 m/s2).
Figure 6. The horizontal gravity disturbance in the simulation area. (a) is the north component of the horizontal gravity disturbance and (b) is the east component of the horizontal gravity disturbance. (1 mGal = 10−5 m/s2).
Sensors 18 00883 g006
Figure 7. The survey lines in the simulation area.
Figure 7. The survey lines in the simulation area.
Sensors 18 00883 g007
Figure 8. The horizontal gravity disturbance on the survey lines. (a) is the north component of horizontal gravity disturbance on survey lines and (b) is the east component of horizontal gravity disturbance on survey lines.
Figure 8. The horizontal gravity disturbance on the survey lines. (a) is the north component of horizontal gravity disturbance on survey lines and (b) is the east component of horizontal gravity disturbance on survey lines.
Sensors 18 00883 g008
Figure 9. (a) is the mean values of the accelerometer outputs every one hundred seconds, and (b) is the bias instability of the accelerometer.
Figure 9. (a) is the mean values of the accelerometer outputs every one hundred seconds, and (b) is the bias instability of the accelerometer.
Sensors 18 00883 g009
Figure 10. Comparing the temperature of the accelerometer with the output of the accelerometer.
Figure 10. Comparing the temperature of the accelerometer with the output of the accelerometer.
Sensors 18 00883 g010
Figure 11. Comparing the temperature of the environment with the output of the accelerometer.
Figure 11. Comparing the temperature of the environment with the output of the accelerometer.
Sensors 18 00883 g011
Figure 12. The estimation error of accelerometer bias in simulation I. (a) is the estimation error of the north component of the accelerometer bias, and (b) is the estimation error of the east component of the accelerometer bias.
Figure 12. The estimation error of accelerometer bias in simulation I. (a) is the estimation error of the north component of the accelerometer bias, and (b) is the estimation error of the east component of the accelerometer bias.
Sensors 18 00883 g012
Figure 13. The estimation error of accelerometer bias in simulation II. (a) is the estimation error of the north component of the accelerometer bias, and (b) is the estimation error of the east component of the accelerometer bias.
Figure 13. The estimation error of accelerometer bias in simulation II. (a) is the estimation error of the north component of the accelerometer bias, and (b) is the estimation error of the east component of the accelerometer bias.
Sensors 18 00883 g013
Figure 14. The estimation error of accelerometer bias in simulation III. (a) is the estimation error of the north component of the bias repeatability, and (b) is the estimation error of the east component of the bias repeatability.
Figure 14. The estimation error of accelerometer bias in simulation III. (a) is the estimation error of the north component of the bias repeatability, and (b) is the estimation error of the east component of the bias repeatability.
Sensors 18 00883 g014
Figure 15. The estimation error of drift rate in simulation III. (a) is the estimation error of the drift rate of the north component, and (b) is the estimation error of the drift rate of the east component.
Figure 15. The estimation error of drift rate in simulation III. (a) is the estimation error of the drift rate of the north component, and (b) is the estimation error of the drift rate of the east component.
Sensors 18 00883 g015
Figure 16. The horizontal gravity disturbance on the short survey lines. (a) is the north component of horizontal gravity disturbance on the short survey lines and (b) is the east component of horizontal gravity disturbance on the short survey lines.
Figure 16. The horizontal gravity disturbance on the short survey lines. (a) is the north component of horizontal gravity disturbance on the short survey lines and (b) is the east component of horizontal gravity disturbance on the short survey lines.
Sensors 18 00883 g016
Figure 17. The estimation error of accelerometer bias on the short survey lines in simulation IV. (a) Is the estimation error of the north component of the accelerometer bias, and (b) is the estimation error of the east component of the accelerometer bias.
Figure 17. The estimation error of accelerometer bias on the short survey lines in simulation IV. (a) Is the estimation error of the north component of the accelerometer bias, and (b) is the estimation error of the east component of the accelerometer bias.
Sensors 18 00883 g017
Figure 18. The estimation error of accelerometer bias on the short survey lines in simulation V. (a) Is the estimation error of the north component of the accelerometer bias, and (b) is the estimation error of the east component of the accelerometer bias.
Figure 18. The estimation error of accelerometer bias on the short survey lines in simulation V. (a) Is the estimation error of the north component of the accelerometer bias, and (b) is the estimation error of the east component of the accelerometer bias.
Sensors 18 00883 g018
Figure 19. The estimation error of the bias repeatability on the short survey lines in simulation VI. (a) Is the estimation error of the north component of the bias repeatability, and (b) is the estimation error of the east component of the bias repeatability.
Figure 19. The estimation error of the bias repeatability on the short survey lines in simulation VI. (a) Is the estimation error of the north component of the bias repeatability, and (b) is the estimation error of the east component of the bias repeatability.
Sensors 18 00883 g019
Figure 20. The estimation error of drift rate on the short survey lines in simulation VI. (a) is the estimation error of the drift rate of the north component, and (b) is the estimation error of the drift rate of the east component.
Figure 20. The estimation error of drift rate on the short survey lines in simulation VI. (a) is the estimation error of the drift rate of the north component, and (b) is the estimation error of the drift rate of the east component.
Sensors 18 00883 g020
Table 1. The coordinate range of the survey lines.
Table 1. The coordinate range of the survey lines.
Line No.Latitude RangeLongitude RangeGrid Spacing
Line 11° N~5° N76° E5 n mile
Line 21° N~5° N77° E5 n mile
Line 31° N~5° N78° E5 n mile
Line 41° N~5°79° E5 n mile
Line 51° N~5° N80° E5 n mile
Line 61° N76° E~80° E5 n mile
Line 72° N76° E~80° E5 n mile
Line 83° N76° E~80° E5 n mile
Line 94° N76° E~80° E5 n mile
Line 105° N76° E~80° E5 n mile
Table 2. The statistical values of horizontal gravity disturbance on the survey lines.
Table 2. The statistical values of horizontal gravity disturbance on the survey lines.
Line No.Mean Value/mGalMedian/mGal
North ComponentEast ComponentNorth ComponentEast Component
1−1.21−10.36−1.55−10.57
2−1.61−6.87−2.42−6.64
3−2.61−3.20−2.971.25
4−7.07−6.27−3.43−3.23
5−5.502.58−8.653.32
6−6.42−3.62−5.70−2.17
7−2.24−2.62−2.01−3.21
8−4.96−2.58−4.44−3.10
9−3.61−5.79−2.92−3.50
106.84−7.928.35−6.84
Mean value of all survey linesMedian of all survey lines
−2.84−4.67−2.97−3.23
Table 3. Simulations on the survey lines.
Table 3. Simulations on the survey lines.
Simulation No.Description
Simulation IOnly bias repeatability is considered and is equal to the mean value of all the survey lines
Simulation IIOnly bias repeatability is considered and is equal to the median value of all the survey lines
Simulation IIIBias repeatability and bias instability are both considered, where bias repeatability is equal to the mean value of all the survey lines and drift rate is 1.31 mGal/day
Table 4. The coordinate range of the short survey lines.
Table 4. The coordinate range of the short survey lines.
Line No.Latitude RangeLongitude RangeSpacing
Line 11° N~1°20′ N76° E1 n mile
Line 21° N~1°20′ N77° E1 n mile
Line 31° N~1°20′ N78° E1 n mile
Line 41° N~1°20′ N79° E1 n mile
Line 51° N~1°20′ N80° E1 n mile
Line 61° N76° E~76°20′ E1 n mile
Line 72° N76° E~76°20′ E1 n mile
Line 83° N76° E~76°20′ E1 n mile
Line 94° N76° E~76°20′ E1 n mile
Line 105° N76° E~76°20′ E1 n mile
Table 5. Simulations on the short survey lines.
Table 5. Simulations on the short survey lines.
Simulation No.Description
Simulation IVOn the short survey lines, only bias repeatability is considered and is equal to the mean value of all the survey lines
Simulation VOn the short survey lines, only bias repeatability is considered and is equal to the median value of all the survey lines
Simulation VIOn the short survey lines, bias repeatability and bias instability are both considered, where bias repeatability is equal to the mean value of all the survey lines and drift rate is 1.31 mGal/day

Share and Cite

MDPI and ACS Style

Tie, J.; Cao, J.; Chang, L.; Cai, S.; Wu, M.; Lian, J. A Model of Gravity Vector Measurement Noise for Estimating Accelerometer Bias in Gravity Disturbance Compensation. Sensors 2018, 18, 883. https://doi.org/10.3390/s18030883

AMA Style

Tie J, Cao J, Chang L, Cai S, Wu M, Lian J. A Model of Gravity Vector Measurement Noise for Estimating Accelerometer Bias in Gravity Disturbance Compensation. Sensors. 2018; 18(3):883. https://doi.org/10.3390/s18030883

Chicago/Turabian Style

Tie, Junbo, Juliang Cao, Lubing Chang, Shaokun Cai, Meiping Wu, and Junxiang Lian. 2018. "A Model of Gravity Vector Measurement Noise for Estimating Accelerometer Bias in Gravity Disturbance Compensation" Sensors 18, no. 3: 883. https://doi.org/10.3390/s18030883

APA Style

Tie, J., Cao, J., Chang, L., Cai, S., Wu, M., & Lian, J. (2018). A Model of Gravity Vector Measurement Noise for Estimating Accelerometer Bias in Gravity Disturbance Compensation. Sensors, 18(3), 883. https://doi.org/10.3390/s18030883

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