Latitude Determination and Error Analysis for Stationary SINS in Unknow-Position Condition

The initial geographic latitude information is the key to the self-alignment of the strapdown inertial navigation system (SINS), but how to determine the latitude when the latitude cannot be obtained directly or in a short time? The latitude determination (LD) methods are introduced, including magnitude method, geometric method, and analytical methods 1 and 2, to solve this situation only by the output of the SINS itself. Simulation and experimental test results validate the efficiency of these LD methods. In order to improve the accuracy of the LD, the error of the LD method is derived through comparative analysis. Based on the relationship between LD error and inertial measurement unit (IMU) bias. Partial bias estimation method is introduced and executed during latitude determination. After compensating the estimated IMU bias, the accuracy of the LD will be further improved. Latitude errors are also affected by the latitude where SINS is located. Comprehensive simulation and experimental tests verify the effectiveness of the method. The IMU determined latitude can not only be used to achieve the self-alignment of the SINS, but also to correct the navigation latitude of the long-term SINS, thereby improving the autonomy and positioning accuracy of the navigation system.

forests, underwater, canyons, and other places that interfere with satellite information, and other application environments where it is difficult to obtain accurate geographic latitudes.
Whereas, the commonly used navigation coordinate frame (n-frame) is established based on the SINS location (latitude, longitude, and altitude), and its axes x n , y n , and z n point to the east, north and up (ENU), known as the local geographic coordinate frame. In addition, the amplitude of gravity acceleration changes with the change of geographic latitude. The gravity acceleration of P point on surface of the Earth can be obtained by g p = 9.78049(1 + 0.0052884sin 2 L − 0.0000059sin 2 2L), where L is the geographic latitude [5]. Because the gravity acceleration is a major function of latitude that latitude information is very important for calculating the magnitude of local gravity acceleration. Moreover, the geographic latitude can be used to decompose the Earth rate vector into n-frame. In n-frame, the projection of the Earth rate vector in three axes is [ 0, Ω cos L, Ω sin L ] T (Ω is the amplitude of the Earth rate vector) respectively, while the acceleration vector of gravity is only in the z-axis direction, and the amplitude is equal to the local gravity [6]. Therefore, the initial geographic latitude is significant and indispensable for INS, so how to determine the geographic latitude in these applications mentioned above is a core task. When there is only an INS, can we only use the static measurement information of the IMU to obtain the geographic latitude?
Inspired by the analytical CA method (also known as the three-axis attitude determination (TRIAD)) [7][8][9], the gyros triads measure the Earth's rate vector while the accelerometers triads measure the local gravity acceleration vector in body frame (b-frame) when the SINS is static [10,11]. Therefore, the geographic latitude is determined by using the static constraints of SINS and the measurement of IMU, and the four geographic latitude determination (LD) methods are elaborated in the paper. In [12], the derivation process of geometric LD is given. Yan's work is mainly to estimate the latitude for initial alignment, the error analysis of LD is not his interest. Aiming at the alignment problem of SINS without latitude, Refs. [13,14] provide a method based on the quaternion calculation of vector of earth rotational angular velocity to achieve initial alignment for both stationary and swaying base. Even if the initial alignment without latitude can be achieved, without the latitude information, the system cannot continue the navigation task. Refs. [15,16] proposed a method for latitude estimation and initial alignment on swaying base. These works did not analyze the latitude estimation error and the latitude determination algorithm is not comprehensive. This paper mainly expends four LD methods for the stationary SINS in unknown-position condition. By the error analysis of LD methods, it is concluded that the estimated latitude error is mainly affected by the biases of the accelerometers and gyros, or some certain direction gyros and accelerometers bias in navigation frame. The analytical CA and gyro bias estimation method by two position alignment was proposed in [17]. Using the two different attitudes of INS, the three-axis gyros bias can be calculated. The estimation and compensation of the gyros biases can effectively reduce the aligned heading misalignment angle and improve the navigation accuracy. Silva analyzes normality and orthogonality error of the analytical CA DCM from b-frame to n-frame C n b then gives the bias of northward and upward gyros and upward accelerometer bias estimation [18]. After the LD process and initial alignment can be performed simultaneously, the estimated latitude can be used synchronously for the INS to quickly start without the initial position. Based on the LD error analysis, the CA bias estimation method (CA-CBE) is introduced. After compensating the bias, the accuracy of latitude estimation is improved obviously. Furthermore, the relationship between latitude error and latitude of four methods is analyzed, so that the best LD method can be selected in different locations. In addition, the determined latitude can correct the position error of the inertial navigation system that works for a long time.
This paper is organized as follows. Section 2 gives four LD methods, whilst the comprehensive error analysis of each LD method is proposed in Section 3. The IMU partial bias estimation algorithm is drawn into Section 4. Sections 5 and 6 present simulation and experiment results. Finally, the conclusion is devoted in Section 7.

Analytic Coarse Alignment Algorithm
The purpose of CA of the SINS is to calculate the attitude matrix, which relates the body frame and the navigation frame. The analytic CA performed based on computing a set of three non-collinear T in the body frame and a set of three non-collinear vectors [g n , ω n ie , (g × ω ie ) n ] T in the navigation frame [3,8]. The body frame is denoted by b with its axes x b , y b and z b pointing to the right, forward, and upward directions (RFU), respectively. The attitude matrix can be expressed as It is assumed that the INS is located at point P on the surface of the Earth. Where , Ω is the magnitude of the earth rate, L is the latitude of P and g p is the gravity acceleration of P somewhere on the earth, C n b can also be expressed as: The DCM of the attitude matrix C n b is given.
where cx represents cosx and sx represents sinx, where x = θ, γ, ψ. Combining Equation (2) with Equation (3), pitch θ, roll γ, and azimuth ψ are calculated by the following equation The analytic CA process is to calculate the coarse attitude for the FA process, which is an indispensable stage for the initial alignment of the static INS. It can be seen from the analytic CA, the geographic latitude is indispensable, because the representation of ω n ie in the n frame has a direct relationship with the geographic latitude; in addition, different latitudes and different heights will affect the local gravitational acceleration g p [5], and the relationship is shown in the following equation. g p = 9.78049(1 + 0.0052884sin 2 L − 0.0000059sin 2 2L+) − 0.000003086h (5) where h is the height of SINS. The main research purpose of this paper is to determine the latitude, and the height can be obtained by barometer.

Magnitude Latitude Determination Method
The measurement of IMU can be expressed into n frame by C n b when the INS is perfectly stationary, as well as the error and noise of inertial devices are not considered [5,6].
The inner product of the two vectors can be obtained as follows.
The latitude L can be can be expressed by the IMU outputs as follows.
Due to the unknown position of the IMU, it is considered that the gravity acceleration at point P cannot be calculated according to the gravity model as Equation (5). Therefore, the global average gravity accelerationḡ = 9.80665 is used here instead of the gravity acceleration g p at point P. Then latitude L can be can be expressed as follows.

Geometric Latitude Determination Method
Under the assumption of the spherical model of the Earth, as shown in Figure 1, the geographic latitude is defined as the complementary angle between the Earth's rate vector and the gravity vector. The geometric method determines the latitude based on this definition. According to [12], the latitude L at P point on the Earth's surface can be expressed by the following equation.
Simultaneously performing sines operators on both sides of the above equation Suppose there is a stationary SINS without considering its error and noise, that is f b = −g b and where |x| denotes the length of the vector x, that is equal to the 2-norm of the vector.

Analytic Latitude Determination Method
The gyros triad measures the vector of the Earth's rate, while the accelerometers triad measures the vector of local gravity acceleration when the SINS is static, that is where C b n is the DCM from n frame to b frame, and C b n = C n b T . Expanding Equation (14) can obtain the following equation.
where g p ≈ f 2 x + f 2 y + f 2 z , and expanding Equation (7) can obtain the following equation.
The initial pitch θ and roll γ can be obtained from Equation (16).
Bringing θ and γ into Equation (16c), then latitude can be calculated as follows, that is called the analytical 1 method.

Analytic 2 Method
The azimuth ψ can be computed combine the Equation (16a) with Equation (17) as follows.
The latitude can be calculated by dividing both sides of Equation (16c) by Equation (16b).
where θ and γ can be obtained from Equation (17).

Latitude Determination Error Analysis
The previous section derived the LD formula based on the IMU error-free assumption, but for the actual INS, the inertial sensor error usually includes bias, random walk, scale factor error and bias drift [1,2]. However, it takes only a few minutes to complete the LD, and calculating the average value can effectively reduce the influence of noise, therefore the in-run biases of the IMU are considered as the main error source for LD. In addition, when the INS is static, the IMU's scale factor error is negligible.
The in-run biases of gyros and accelerometers are δω b and δ f b , respectively, that is To unify the error analysis of the latitude estimation, the equivalent biases in the n frame are analyzed as the error source, since the b frame can be arbitrarily oriented with respect to the n frame. For example, The expressions become where the subscripts N, U represent the amount of error in the ENU direction.

Magnitude Latitude Determination Method
The Equation (10) of magnitude LD method can be rewritten as follows Substituting Equation (21) into Equation (23) yields Assuming that the latitude error is a small angle error, where δL is the latitude error, expand the above equation and ignore the second-order error term. The above equation combined with Equation (22) can be simplified to obtain the following equation.
where k g = g/ḡ is the gravity acceleration model error scaling factor. The main error source of magnitude LD method is gravity acceleration model error, which is caused by local gravity replaced by the average gravity acceleration. The northward and upward bias of the accelerometer and the upward bias of the gyros also affect the LD error, as well as the location itself.

Geometric Latitude Determination Method
Substituting Equation (21) into Equation (13) yields The following is the derivation and approximation of the first square root reciprocal of Equation (26). In the process of error derivation, only the first-order error terms are retained and the second-order error terms are ignored.
Similarly, the derivation and approximation of the second square root reciprocal of Equation (26) can be obtained as follows.
Then the error of geometric LD method becomes By substituting Equation (22) into Equation (29), the following geometric LD error equation can be obtained.
By error analyzing the geometric LD method, according to Equation (30), it can be found that the latitude error is mainly affected by the upward and northward bias of gyros and northward bias of accelerometer, as well as being modulated by the latitude itself.

Analytic 1 Latitude Determination Method
From Equation (15), it is easy to get sine and cosine expressions of pitch θ and roll γ related to IMU output, and substitute these expressions into Equation (18). The LD equation of analytic 1 method can be rewritten as follows.
Then the error of latitude can be written as follows By rearranging Equation (32) in combination with Equations (22) and (27), the latitude error equation of analytic 1 LD method can be expressed as follows According to Equation (33), both the northward bias of accelerometer and the upward bias of gyros contribute to the latitude error of the analytic 1 LD method.

Analytic 2 Latitude Determination Method
Similar to the derivation process of latitude error in analytic 1 method, it is easy to obtain the sine and cosine expressions between azimuth ψ and IMU output from Equation (19). Associating Equations (15), (19) and (20), the LD equation of analytical 2 method can be rewritten as follows: where The derivation of the latitude error is rather complicated. The LD error equation of Analytic 2 are not carried out. As a suggestion for future work, we intend to derive error equations similar to other LD methods. So far, the error analysis is carried out by using the subsequent simulation and experimental results.

Imu Partial Bias Estimation Algorithm
Inspired by the IMU bias estimation method in [18], the CA-CBE method performed after the latitude determination. The specific calculation process is to calculate the analytical CA matrixĈ n b according to Equation (2) whereÊ s is a symmetric matrix representing the attitude matrix normality error vectorη and orthogonality error vectorsô, respectively, that iŝ According to the bias estimation Equations (45)-(47) in [18], the upward bias of accelerometer, northward bias and upward bias of the gyros can be computed as follows.
From the previous error analysis of the LD method, it can be found that the LD error mainly originates from the IMU bias, and the northward bias of accelerometer, upward and northward biases of gyro are the main error sources. By estimating and compensating for the partial bias of the IMU, the accuracy of the LD can be improved.

Simulation Results
In order to verify the proposed LD methods and error analysis, corresponding simulation experiments were carried out. The simulation experiments included four aspects as follows.

Multiple Latitude Determination Tests at the Same Location
For the purpose of the test, a navigation grade INS with three gyros and accelerometers was taken as a candidate simulation SINS. The accelerometer data sets were corrupted by 100 µg of constant biases and 10 µg/ √ h of random walk, and the gyros data sets were added 0.01 • /h of constant biases and 0.001 • / √ h of random walk, respectively [19,20]. The initial position of the simulation SINS was located in Beijing of China with a latitude of 39.97 • (N), longitude of 116.34 • (E) and altitude of 50 m. Then, the 5 min data sets sampled with 100 Hz were generated.
The 500 Monte Carlo LD simulation tests were performed at the same location. The statistical results are listed in Table 1. The latitude estimated convergence curves by four different LD methods are shown in Figure 2. According to Table 1, it can be seen that the latitude error calculated by the magnitude method at the simulated latitude was 9.91 arcmin, while the latitude error calculated by the other three methods was relatively small (latitude error was less than 3.36 arcmin). The simulated latitude error was very close to the error calculated by the LD error equations. The latitude error obtained by the magnitude method was the largest, and the geometrical method can obtain the smallest LD error in this simulation condition. As shown in Figure 2, the estimated latitude converged at 50 s, and the latitudes estimated by the three LD methods eventually converged to the reference latitude within 5 min. It shows that the LD methods can effectively determine latitude using static IMU output.

The IMU Biases Estimation and Biases Compensation Latitude Determination Test
The CA-CBE algorithm was used to simultaneously estimate the northward, upward gyros bias, and upward accelerometer bias according to Equation (37), and the simulation conditions were the same as Section 5.1. The bias estimated curve is shown in Figures 3-5. The final value of the estimated bias is listed in Table 2. As shown in Table 2, the estimated upward biases of the accelerometer by three LD methods are 96.97 µg, 92.98 µg and 94.84 µg. Similarly, this also can estimate the north and upward gyro bias. However, only the latitude estimated by the geometric LD method can be used to estimate the upward and northward bias of the gyros, as shown in bold values in Table 2. While the latitude calculated by analytic methods cannot be estimated. This also verified that the latitude error estimated by the geometric method was the smallest at the simulated position.
In order to further verify the bias estimation algorithm, a comparison test of LD method with uncompensated and compensated IMU bias was performed. Table 3 shows the latitude errors of uncompensated bias and compensated bias. After the bias was compensated, the latitude error was significantly reduced. Especially for the geometric LD method, it is obvious that after compensating the estimated bias, the LD accuracy was greatly improved.

The Latitude Determination Test at Different Positions
The LD error equation indicates that the LD error is affected by the latitude where SINS is located. Suppose SINS is located in the northern hemisphere with a latitude range of 0 • to 80 • , and the latitude starts from 0 • with a step size of 1 • . The other simulation conditions are the same as Section 5.1 excepting the location of the SINS. At each simulated latitude, 500 Monte Carlo LD tests are performed, and the average latitude estimation errors are plotted in Figures 6-9. The latitude error computed from the Equations (25), (30) and (33) are ploted in Figures 7-9. From the Figures 6-8, we can conclude that the simulation results are consistent with the error analysis results. Figure 9 only shows the latitude error results of analytical method 2, which provides a numerical reference for us to determine the latitude. The latitude error of the magnitude method and the analytical 1 method is mainly affected by U /Ω cos L, so the estimated latitude error will increase rapidly as the location of the SINS increases. However, the magnitude method is affected by g p ∇ U tan L/ḡ, so its latitude error increases much more than the others method. The latitude estimation error of the geometric method is between −2 arcmin and 2 arcmin. The latitude error of geometric method decreases with increasing latitude, and the estimated latitude error is zero at a certain latitude point. Analytical method 2 is less affected by latitude than the others method, its latitude error is smaller than other methods when SINS is in a high latitude area. The simulation results in the southern hemisphere have symmetrical transparent feature.

The Latitude Determination Test by Different Level of IMU
To verify the applicability of the LD methods, the simulation experiments on IMUs of different accuracy levels were conducted. The constant bias of the gyro was 1 • /h, 0.5 • /h, 0.02 • /h, 0.001 • /h, respectively, and the corresponding constant bias of the accelerometer was 1000 µg, 500 µg, 100 µg, 10 µg, respectively. The simulation location was the same as Section 5.1, 500 Monte Carlo LD tests were performed, and the average latitude estimation errors are listed in Table 4. Obviously, error of LD will shapely decline with performance of IMU incline, the simulation LD error is also close to the error computed by the LD error equations.

LD Error(arc-min) Simulation Error Equation Error
Biases of IMU

Experimental Results
To further verify the LD methods in practical applications, the experiments was carried out on an SINS manufactured in our laboratory. The SINS was developed using three dithered ring laser gyroscopes and three quartz accelerometers, and the detail specifications of the SINS are shown in Table 5. The experiments mainly included two aspects, one was multiple LD tests at the same location and the other was LD test at multiple locations.

Multiple Latitude Determination Experiments at the Same Position
The multiple LD testes were implemented in the (Inertial Technology and Integrated Navigation Laboratory) New Main Building of BeiHang University in Haidian District Beijing, China. The experiment equipment mainly included an SINS, a 24 V DC power supply, a marble platform and a laptop for raw data sampling at 200 Hz, as shown in Figure 10. The geographic latitude of the site were 39.977200 • (N). The system needed to be warmed up for 30 min before the start of the test. Then the LD test was implemented by collecting 5 min static raw data multiple times.
According to the LD algorithm and its error analysis, the geometric method, the analytical method 1 and the analytical method 2 were chosen to compute latitude. The estimated geographic latitude of 7 times is shown in Table 6 and one of the curves of three LD methods with time is shown in Figure 11, respectively. In the process of calculating latitude, the CA-CBE algorithm was used to estimate the bias of northward gyro, the bias of upward gyro and the upward bias of the accelerometer, the estimated biases are shown in Table 7.
Then, the LD process was performed again after compensating the estimated biases. The statistical latitude estimation errors for uncompensated and compensated biases are listed in Table 8. As is evident from Table 8, the average latitude errors of the proposed methods were −0.83 arcmin, −1.23 arcmin and −1.02 arcmin, which is consistent with the previous error analysis results. After compensating the biases, the estimated latitude accuracy was significantly improved, which also proves the effectiveness of the bias estimation method. The mean values of latitude errors after compensating the biases were reduced to 0.03 arcmin, 0.21 arcmin and 0.29 arcmin, falling by 96.37%, 82.88% and 71.48% correspond to the geometric method, the analytical 1 method and the analytical 2, respectively.

Latitude Determination Experiments at the Multiple Positions
To further validate the capability, the same RLG SINS is applied in LD experiments at the multiple positions. We chose eight different positions in China for LD experiments. The latitudes of the multiple positions ranged from 26.0 • to 40.7 • . These eight positions were located in Fuzhou City, Xi'an City (2 points), Lianyungang City, Zhengzhou City, Beijing City (2 points) and Huludao City, respectively. The estimated geographic latitudes of the eight different positions are shown in Table 9. Figure 12 illustrates the latitude error computed for the three LD methods, and the x-axis presents the eight positions increasing with latitude meanwhile the y-axis is the latitude error. It is indicated from Figure 12 that the variation of latitude estimation error with increased latitude occurred in the same manner as in . This is consistent with the latitude error equation. It is consistent with the previous LD error equation.

Conclusions
This paper mainly solves the problem that the INS cannot be started in some unknown positions (especially the latitude). This technology can be used to determine the latitude information of the INS, and the LD is performed only based on the output of the static inertial navigation system itself. Four algorithms of LD are comprehensively introduced, and the relationship between the LD error, IMU bias and INS position is emphatically deduced. Simulation and experiments verify the effectiveness of the LD algorithm, as well as the derived error equations. For gyros and accelerometers with bias of 0.01 • /h and 100 µg, the LD error of 500 monte carlo simulation of geometric method is less than 1 arc-min, while for gyros and accelerometers with bias of 1 • /h and 1 mg, the LD error is less than 0.5 • , indicating that the feasibility of extending the method to lower grade IMU. The "CA-CBE" bias estimation algorithm is introduced into the LD stage, which can effectively estimate the partial bias of IMU. After the estimated IMU bias is compensated, the accuracy of the LD is further improved. The LD error of actual SINS is less than 1.5 arc-min when the bias of IMU is uncompensated, and the LD error is less than 0.3 arc-min after compensated bias of IMU. Finally, the LD error is affected by the position of the inertial navigation system. The typical navigation-level SINS is used for imulation and experimental verification the different latitudes. The results are consistent with the derived error formula. The application of this technology can further improve the autonomy of INS.
Author Contributions: S.W. and W.C. proposed the initial idea and wrote the paper. G.Y. and S.W. conceived and performed the experiments. W.C. and L.W. reviewed and edited the manuscript. All authors have read and agreed to the published version of the manuscript. Acknowledgments: In this section you can acknowledge any support given which is not covered by the author contribution or funding sections. This may include administrative and technical support, or donations in kind (e.g., materials used for experiments).

Conflicts of Interest:
The authors declare no conflict of interest.