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

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.


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 Figure 1. The definition of gravity disturbance vector, gravity disturbance and deflection of vertical, and this figure is adapted from [11].
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.

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 n E 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 γ:

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: 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: Usually the geocentric colatitude ϑ is used in the SHM calculation: Then the horizontal gravity disturbances can be obtained as: The disturbed gravitational potential T is the solution of Laplace equation in the ellipsoidal coordinate frame, which can be represented as follows [11] 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, *

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: 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: Usually the geocentric colatitude ϑ is used in the SHM calculation: Then the horizontal gravity disturbances can be obtained as: The disturbed gravitational potential T is the solution of Laplace equation in the ellipsoidal coordinate frame, which can be represented as follows [11]: a r n C * nm cos mλ + S nm sin mλ P nm (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 * nm and S nm are the coefficients of the SHM, P nm (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) a r n m −C * nm · sin mλ + S nm · cos mλ P nm (cos ϑ) (9) Substituting Equations (8) and (9) into Equations (5) and (6), the formulas of calculating horizontal gravity disturbance are obtained:

Reference Coordinate Frames
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 ω e ie = [0 0 Ω] T , Ω is the Earth's rotation angular velocity, as shown in Figure 3a.

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 x n −y n −z n axes respectively point towards North-Up-East, as shown in Figure 3a. It should be noted that z n is collinear with the normal gravity vector which points towards the center of the reference ellipsoid.

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.

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 { x n ,y n ,z n } and x n also points towards north, but z n is collinear with the true gravity vector and y n can be determined based on the right-hand rule, as shown in Figure 3b.

Reference Coordinate Frames
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 = 0 0 Ω , Ω is the Earth's rotation angular velocity, as shown in Figure 3a.

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.

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.

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

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. 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: where n r is the projection of arbitrary vector r in n coordinate frame, b r is the projection of vector r in b coordinate frame.

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]: where b nb ω is the body angular rate with respect to the navigation frame n and is given by [2]: b ib ω is the body angular rate with respect to inertial frame and is measured by the gyroscopes. n ie ω is the earth rotational rate and is given by [2]: 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]: where ω b nb is the body angular rate with respect to the navigation frame n and is given by [2]: ω b ib is the body angular rate with respect to inertial frame and is measured by the gyroscopes. ω n ie is the earth rotational rate and is given by [2]: ω n ie 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 ω n ie is given by [2]: R M and R N are the meridian and transverse radius of the ellipsoid curvature, respectively. v n E 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]: where v n e 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: 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]:

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]: where δϕ is the attitude error vector, δα is the roll error, δφ is the yaw error, δβ is the pitch error and δω b ib 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): 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): where v n e is the calculated velocity-containing error. ω n ie and ω n en are angular velocities containing error, γ n is the normal gravity vector, and δf n is the white noise in accelerometer output. C n b is the DCM containing error and is defined as follows: The velocity error equation without accelerometer bias is derived by subtracting Equation (22) from Equation (23): where δv n e is the velocity error vector, δv n N , δv n U and δv n E 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 n U and f n E are the components of f n in the n frame. Initial alignment is usually implemented when the INS is at static or mooring, f n can be regarded as the opposite of the true gravity vector: 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: 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) .
The new velocity error equation is obtained by subtracting Equation (22) from Equation (32): where b n a is the accelerometer bias and b n a,N , b n a,U and b n a,E 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: 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): 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.

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): 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 n e is the difference of velocity provided by GNSS, ω n ie and ω n en can be calculated by substituting the velocity and position provided by GNSS into Equation (15) and Equation (16). C n b 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: where ∆g n N,E is the true value of horizontal gravity disturbance, ∆ g n N,E 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: The measurement error of horizontal gravity disturbance can be obtained by subtracting the measured value from the true value: where δg n N and δg n E 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 b n a,N andb n a,E 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.

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): 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 Earth's gravity field, but the contribution of measurement noise to measurement error may be associated with the gravity field.

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): It should be noted that the y-axis of n′ coordinate frame is collinear with the true gravity vector, the horizontal components of n′ g 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.
Rotation Axis Figure 5. The geometric relationship between the two navigation coordinate frames.
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 satisfies four constraints: (1) In the plane n n x z − ; (2) Pass through the origin of the two coordinate frames; (3) Be orthogonal to the plane n n y y ′ − ; (4) u is unit vector; Based on constraint 1: Based on constraint 2, k is a constant to be determined: Based on constraint 3: Based on constraint 4: 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: Based on constraint 2, k is a constant to be determined: Based on constraint 3: Based on constraint 4: According to Equations (46)-(49), the rotation axis u can be determined: ϑ is the angle between γ n and g n , and can be calculated based on vector product: Based on the rotation axis and the rotation angle, the corresponding quaternion and DCM can be obtained [2]: Equations (54) and (55) can be simplified due to u y being zero: Based on the Equations (50)-(57), it can be seen that C n n is a function of the horizontal gravity disturbance: and Equation (18) can be approximately rearranged as: Now, when only considering the measurement noise, the DCM based on the measured values is supposed to has the following form: Substituting Equation (60) into Equation (59): The expression of measurement noise is obtained: Suppose δC n n has the following form: The model of measurement noise is obtained, this is the foundation of estimating accelerometer bias in this paper: Substituting Equation (64) into Equation (43), the model of measurement error is rearranged as:

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 1 and δq 3 are the errors caused only by the measurement noise. Substituting Equation (66) into Equation (57), we can get: The model of measurement error is rearranged as: The next work is to derivate the expression of δq i,i=0,1,3 through the following partial differential equations: According to the Equations (50)-(57), these partial derivatives are obtained:

Estimation Model of the Accelerate Bias
Finally, based on the model of measurement noise, the model of measurement error is rearranged as follows: δξ 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 .

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: where z 2N×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 2N rows and one column. H 2N×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 n a,N andb n a,E are the estimations of the accelerometer biases:

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.
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 2N rows and one column.
is the observation matrix which has 2N rows and four columns,

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

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 as 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 as is constant throughout an IMU operating period but varies from run to run. The dynamic component denoted by b ad is also known as the bias instability and varies over time and temperature: 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 Figures 10 and 11. Comparing Figures 10 and 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,  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 Figures 10 and 11. Comparing Figures 10 and 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) 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 Figures 10 and 11. Comparing Figures 10 and 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,  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: 10 20 log magnitude of accelerometer bias SNR dB magnitude of measurement noise 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. Table 3. Simulations on the survey lines.

Simulation No. Description
Simulation I Only bias repeatability is considered and is equal to the mean value of all the survey lines Simulation II Only bias repeatability is considered and is equal to the median value of all the survey lines Simulation III 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 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 Figures 12 and 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 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: 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. Table 3. Simulations on the survey lines.

Simulation No. Description
Simulation I Only bias repeatability is considered and is equal to the mean value of all the survey lines Simulation II Only bias repeatability is considered and is equal to the median value of all the survey lines Simulation III 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 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 Figures 12 and 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: where t is time,   In Simulation III, some modifications are needed on Equations (83)-(86) to estimate the bias repeatability and drift rate simultaneously: where t is time,  In Simulation III, some modifications are needed on Equations (83)-(86) to estimate the bias repeatability and drift rate simultaneously: x = δξ δη b n as,N b n as,E k n ad,N k n ad,E T (93) where t is time,b n as,N andb n as,E are the estimations of bias repeatability,k ad,N andk ad,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 Figures 14 and 15.
In Simulation III, the accelerometer bias is time varying, so the x-axis is measurement noise instead of SNR. And from Figures 14 and 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. 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 Figures 14 and 15.
In Simulation III, the accelerometer bias is time varying, so the x-axis is measurement noise instead of SNR. And from Figures 14 and 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. 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 Figures 14 and 15.
In Simulation III, the accelerometer bias is time varying, so the x-axis is measurement noise instead of SNR. And from Figures 14 and 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. 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.  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 Figures 17 and 18, it can be seen that the maximum estimation error of accelerometer bias is less than 1 10 − % 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 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 Figures 17 and 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 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  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.

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

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

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.