Modeling and Implementation of Multi-Position Non-Continuous Rotation Gyroscope North Finder

Even when the Global Positioning System (GPS) signal is blocked, a rate gyroscope (gyro) north finder is capable of providing the required azimuth reference information to a certain extent. In order to measure the azimuth between the observer and the north direction very accurately, we propose a multi-position non-continuous rotation gyro north finding scheme. Our new generalized mathematical model analyzes the elements that affect the azimuth measurement precision and can thus provide high precision azimuth reference information. Based on the gyro’s principle of detecting a projection of the earth rotation rate on its sensitive axis and the proposed north finding scheme, we are able to deduct an accurate mathematical model of the gyro outputs against azimuth with the gyro and shaft misalignments. Combining the gyro outputs model and the theory of propagation of uncertainty, some approaches to optimize north finding are provided, including reducing the gyro bias error, constraining the gyro random error, increasing the number of rotation points, improving rotation angle measurement precision, decreasing the gyro and the shaft misalignment angles. According them, a north finder setup is built and the azimuth uncertainty of 18” is obtained. This paper provides systematic theory for analyzing the details of the gyro north finder scheme from simulation to implementation. The proposed theory can guide both applied researchers in academia and advanced practitioners in industry for designing high precision robust north finder based on different types of rate gyroscopes.


Introduction
The rate gyroscope (gyro) north finder is an orientation instrument that finds the azimuth by measuring a projection of the earth rotation rate onto the gyro's sensitive axis. It has been widely used in e.g., aiming techniques, surveying techniques, and navigation of autonomous ground vehicles [1][2][3]. It is also capable of providing the required azimuth reference information when the GPS signal is blocked, for instance in underground environments. Usually, north finding with an uncertainty of 206"(206" = 1 milliradian) is required in these applications [4]. However, more accurate north finding systems are always popular in practical application, e.g., if the aiming distance is 1 km, the aiming error is about 1 m when the azimuth uncertainty is 206", and the aiming error is about 0.17 m when the azimuth uncertainty is 36" . Thus improving the north finding system precision is valuable to both applied researchers in academia and advanced practitioners in industry. In current gyro north finder systems, different types of gyros are used to measure the angular rate with respect to an inertial frame of reference, including Micro-Electro-Mechanical System (MEMS) gyros, ring laser gyros, fiber optic gyros, resonant optical gyros and dynamically tuned gyros [5]. According to their general usage categories, gyros can be sorted as consumer, tactical, navigational and strategic type. Recently, the structures of different gyros have been comprehensively analyzed to constrain the gyros' bias error, scale-factor error and random error, e.g., MEMS gyros [6][7][8][9], ring laser gyros [10,11], fiber optic gyros [10,12], resonant optical gyros [13][14][15][16] and dynamically tuned gyros [17][18][19]. Other very active areas include constraining the gyro bias error and modeling of temperature drift through different kinds of filter methods [20][21][22][23][24][25][26]. However, only few works present the details of building a complete gyro north finder system that provides high precision azimuth reference information.
It is a challenge to design a high precision robust north finder system and to provide the corresponding azimuth measuring model. Many factors should be taken into consideration. Firstly, the model should take the gyro's characteristics into account, including the gyro bias, the gyro random error and the scale-factor error. Secondly, the relationship of the gyro and shaft misalignment angles against the azimuth model needs to be clarified. The gyro misalignment angle is the mechanical install error between the shaft (rotation axis) and the gyro axis; the shaft misalignment angle is the tilt between the practical shaft and the vertical one. Finally, the azimuth measurement precision is affected by several elements such as the measurement points and encoder precision; thus these elements should be considered by the system model.
Currently, two of the most commonly used approaches to attain high precision robust azimuth results are two-point gyro north finding [27] and multi-point continuous rotation gyro north finding [4,[28][29][30]. In the former case, the gyro is rotated by 180 • . By subtracting the azimuth measurement results at the two positions, the additive bias error and scale-factor error is mitigated. However, it is not robust to the gyro random error. For achieving high azimuth measurement precision, the latter approach is more often used, because it is not only robust to the gyro bias error and the scale-factor error, but also to the random error. In particular, Prikhodko provides a comprehensive analysis about how the gyro misalignment, the gyro bias error, and the number of measurement points affect the north finding precision, but does not provide the azimuth measuring model [4]. Recently, a north finding system consisting of an uniaxial MEMS accelerometer and a fiber optic gyro is proposed, and a translocation method is adopted to eliminate sensor constant drift and scale factor error [31]. Bojja improves this system by replacing the uniaxial MEMS accelerometer with a two-axis MEMS accelerometer, and an indexing method is used to remove the sensor bias [32]. However, these methods still possess three shortcomings. The first one is that the gyro detects not only the earth rotation rate but also the turntable rotation rate due to the gyro misalignment and the shaft misalignment. The second one is that the vibration produced by the continuous rotation corrupts the gyro outputs. Finally, there exists no generalized mathematical model which analyzes the elements that affect the azimuth measurement precision. Hence, proposing an appropriate azimuth measuring scheme and providing the corresponding azimuth measuring mathematical model for designing and evaluating the north finder are urgently needed.
In this paper we propose a multi-position non-continuous rotation gyro north finder scheme to measure the azimuth between the observer and the north direction. We are the first to deduce the mathematical model of the gyro outputs against the azimuth at each rotation position explicitly. Our approach demonstrates the effects of the gyro and the shaft misalignments on the gyro outputs.
Combining the gyro outputs model and the theory of propagation of uncertainty, we thus obtain an optimized azimuth measurement scheme.
The rest of the paper is organized as follows. Section 2 presents the gyro north finder structure and the multi-position non-continuous azimuth measuring scheme. Section 3 gives the gyro outputs model. The optimized north finding scheme is described in Section 4. Experimental results are discussed in Section 5. Finally, conclusions are drawn in Section 6. Figure 1 demonstrates the layout of the north finder whose length, width and height are 420 mm, 420 mm and 416 mm respectively. Inside a pedestal, an inclinometer, a gyro, a platform, a shaft, a motor, and an encoder (see Figure 1a) are placed. The inclinometer and the gyro are mounted on the platform that is fixed to the shaft. The pedestal is shored up by three supporting legs whose heights are adjusted by the three knobs (see Figure 1b). These devices are divided into two categories according to their functions. The first category measures the shaft misalignment angle and aligns the shaft. It includes the inclinometer, the three knobs and the three support legs. The inclinometer is used to measure the angle of the platform with respect to the horizontal plane, and the shaft tilt is adjusted by the three knobs using the inclinometer outputs. Secondly, the gyroscope, the platform, the shaft, the motor, and the encoder are used to measure the earth rotation rate onto the gyro sensitive axis at each rotation position. The gyro detects the projection of earth rotation rate onto its sensitive axis. The motor controls the gyro (the platform and the shaft) to move non-continuously from position 1 to position n. The encoder records the rotation angle of the gyro. The non-continuous rotation scheme samples specify the gyro outputs at appointed positions. As shown in the first graph in Figure 2, the vector OA points towards the observer i.e., position 1; the vector OG denotes the gyro sensitive axis; the vector ON point to north direction and ψ indicates the azimuth of the observer against the north direction. The angle between the observer and the current gyro axis is denoted as γ 1 . As the gyro rotates non-continuously around the shaft by 360 • in n equal angular increments, we obtain the angle γ i (1 ≤ i ≤ n) and the gyro outputs at each position from 1 to n. The azimuth is achieved by fitting the gyro outputs at these positions using the azimuth measurement model. The following section describes our model in detail.

Modeling of the Gyro Outputs
The gyro detects the projection of the earth rotation rate on its sensitive axis and outputs its value. In this section, our mathematical model describing the relationship between the gyro outputs and the azimuth at each appointed position, i.e., position 1 to n, is proposed. Using this model, the azimuth is easily obtained from the gyro outputs. Subsection 3.1 proposes an azimuth measuring model in ideal condition. Since the azimuth measurement result is strongly affected by the gyro misalignment and the shaft misalignment, Subsection 3.2 presents a new measuring model with the former element taken into account, and Subsection 3.3 further improves this model by taking both the affecting elements into consideration.

Gyro Outputs Modeling Under Ideal Conditions
The aim of the north finder is to detect the azimuth ψ between the observer (first position of the gyro) and the north direction. Under ideal conditions (without the gyro and shaft misalignments), the azimuth measurement is performed in a local horizontal plane with the gyro axis kept parallel to the horizontal plane. In Figure 3, the sphere is the earth which rotates around itself with an angular rate ω e . The north finder is placed at point O with latitude ϕ. OXYZ is the geographical coordinate system, where the axes OX, OY, and OZ coincide with east, north, and up direction, respectively. Under ideal conditions, the gyro axis OG is parallel to the OXY horizontal plane at each rotation position i. At position 1, the gyro axis OG coincides with the direction OA. It is clear from Figure 3 that the earth rotation angular velocity ω can be expressed in the geographical coordinate system as a column vector: where ω e is the magnitude of the earth angular velocity, with ω e = 15.041 • /h. In order to detect the azimuth ψ, the earth angular velocity (i.e., the column vector ω) is rotated counterclockwise by ψ around the OZ axis. Using the right hand rule (all the following rotation matrices follow these rule), the rotation matrix is given by It should be noted that in our system, the OXYZ coordinate system is fixed, while the vector ω is rotated. Then the rotated vector ω o is obtained by using the matrix multiplication R Z (ψ)ω: Here, the projection of the vector ω o on the OXY horizontal plane coincides with OA. When the gyro axis coincides with OA, the gyro output at this position is given by where e b and e r denote the gyro bias error and the gyro random error, respectively. At a specific latitude, the azimuth ψ is achieved from Equation (4) if the gyro output ω og is known. In practice, however, since only one gyro output is used, it requires the high bias stability of the gyro. Otherwise the small angular rate ω e cos ϕ cos ψ will be too small compared to the gyro bias error e b . Therefore, it is necessary to sample the gyro outputs at several positions. Hence, a multi-position non-continuous rotation north finding scheme is used in our system. As described in Figure 2, the gyro is rotated counterclockwise by γ i around axis OZ from position 1 to position n. We assume that the vector ω o is rotated along side with the gyro. The corresponding rotation matrix is At position i, the vector ω r is obtained by multiplying R Z (γ i ) by ω o : The projection of ω r on the gyro axis at position i hence becomes Equation (7) is the mathematical model of the gyro outputs against azimuth at position i without taking the gyro and shaft misalignments into consideration.

Gyro Outputs Modeling with Gyro Misalignment
Conversely, in the north finding system, the gyro misalignment [4] and the shaft misalignment [33] usually do exist and corrupt the azimuth measurement result. As shown in Figure 4 (the sphere is in this case not the earth, but used to facilitate the understanding of the misalignments), under ideal conditions, the shaft is aligned vertically to the horizontal plane and is represented here as the vector OR. The gyro sensitive axis OG lies in the horizontal plane and points in the observer direction. When the shaft rotates, the gyro axis always remains in the horizontal plane, i.e., the blue part. If the gyro is misaligned, there exist an angle ε between the gyro sensitive axis OG g and the horizontal plane. It is defined as the gyro misalignment angle. As the shaft OR rotates, the trajectory of OG g forms a circular cone (the green part) and ε remains the same at each rotation position. If the shaft is misaligned, the condition becomes complex. Assuming OR 1 is the misaligned shaft, the orange plane will be the moving trajectory of the gyro sensitive axis OG s as it rotates. Here, η is the angle between OR and OR 1 and is defined as the shaft misalignment angle, and α i is used to denote the angle of the gyro axis to the OXY horizontal plane. Unlike the case when the gyro is misaligned, the angle between the gyro axis OG s and the OXY horizontal plane is different at different rotation positions. Since the characteristics of the gyro misalignment and the shaft misalignment are different, the two misalignments will be discussed separately, and eventually a generalized gyro outputs model will be presented. The gyro axis outputs model with only the gyro misalignment (ε = 0, η = 0) is firstly presented. Then the model with also the shaft misalignment (ε = 0 or ε = 0, η = 0) is proposed. Figure 4. The gyro axis's trajectory with the gyro misalignment or the shaft misalignment when the shaft is rotated counterclockwise around itself for 360 • (the sphere is not the earth, but is used to facilitate understanding of the misalignments).
In this part, η = 0 is assumed, i.e., only the gyro misalignment angle ε is considered. The gyro misalignment angle is obtained by rotating the gyro axis counterclockwise by −ε around the X axis. Likewise, we rotate the vector ω r counterclockwise by −ε around the X axis, and then a new vector ω g is obtained. Now the gyro outputs are equal to the projection of the vector ω g on the gyro axis. The rotation matrix is The vector ω g is obtained by multiplying R X (−ε) by ω r : Thus the earth rotation rate detected by the gyro axis is equal to the projection of ω g on the gyro axis:

Gyro Outputs Modeling with Shaft Misalignment
In this part, the shaft misalignment condition (η = 0) will be explained. The gyro axis is obtained by rotating the gyro axis counterclockwise by −α i around the X axis. Similarly, we rotate the vector ω r counterclockwise by −α i around the X axis, and then a new vector ω s is obtained. Now the gyro outputs are equal to the projection of the vector ω s on the gyro axis. The rotation matrix is The vector ω s is obtained by multiplying R X (−α) by ω r : Thus, the earth rotation rate detected by the gyro axis is equal to the projection of ω s on the gyro axis: As shown in Figure 4, the value of α i varies at different rotation positions i when η = 0. Hence, the calculation of α i is discussed. Based on the value of ε (ε = 0, ε = 0), α i is divided into two cases. As shown in Figure 5, the first case is ε = 0: OG s1 denotes the corresponding gyro axis at position 1, and α 1 = η. The trajectory of gyro axis under this condition (η = 0, ε = 0) is a plane. The second case is ε = 0: the corresponding gyro axis is OG sg1 , the angle G sg1 OG s1 is ε, and α 1 = η −ε. The trajectory of the gyro axis under this condition (η = 0, ε = 0) is a circular cone. Figure 5. The gyro axis's trajectory with the shaft misalignment and the gyro misalignment when the shaft is rotated counterclockwise around the shaft for 360 • (the sphere is not the earth, but used to facilitate understanding of the misalignments).
Here, the second case is discussed. The first case is achieved by defining ε = 0. The gyro axis at position i is OG sgi , and OG sgi is defined as a unit vector, and then point G sg1 is located at (0, cos(η−ε), −sin(η−ε)). Point G sgi 's trajectory is a circle with its center O 1 at (0, sinε sinη, sinεcosη), and point G sgi 's trajectory forms the shaft misalignment plane. The normal vector of the shaft misalignment plane is given by n = (0, sin η, cos η) Based on Equation (14) and point G sg1 , the shaft misalignment plane equation equals In the rotating process, if the shaft O 1 R 1 is rotated γ i counterclockwise around itself, so does O 1 G sg1 . Assuming point G sgi = (x i , y i , z i ), then OG sgi is a unit vector, point G sgi is on the shaft misalignment plane, and the angle between vector O 1 G sg1 and vector O 1 G sgi is γ i . Based on these conditions, the following set of equations about point G sgi is thus obtained: The coordinate of G sgi is found by solving Equation (16): Therefore, cos α i is equal to the cosine of the vectors OG sgi and OA i (A i is in OXY plane, and A i = (−sinγ i , cosγ i , 0)): Substituting Equation (18) into Equation (13), the gyro axis outputs model with the shaft and gyro misalignments is achieved.

Estimate Azimuth Uncertainty
In order to estimate the measurement uncertainty of the north finding system, each critical measurement component should be identified and its uncertainty must be quantified. The parameters in the north finder system include the gyro bias error, the gyro random error, the number of measurement points n, the encoder uncertainty δ γ , the gyro misalignment angle ε, and the shaft misalignment angle η. Combining the gyro outputs model and the theory of propagation of uncertainty, these parameters are analyzed separately, and finally the overall azimuth uncertainty is presented.
Initially, the gyro measurement uncertainty is analyzed. According to Equation (19), the derivative of ω sgi to ψ with η = ε = 0 is ∂ω sgi ∂ψ = −ω e cos ϕ sin(γ i + ψ) The gyro bias is a nonzero output when the input is zero. Usually, the bias error is very small in a short period time when the temperature is stable. The gyro random error is a fluctuating output when the input is the same. Here we take the gyro bias error e b and the gyro random error e r as gyro measurement uncertainty [5]. At position i, the gyro measurement uncertainty's effect on ψ's uncertainty is described as where σ ω is the gyro measurement uncertainty. In one circle measurement from position 1 to position n, the total azimuth uncertainty σ ψωsg caused by the gyro measurement uncertainty is In Figure 2, the rotation angle is equal at every time in one circle, and the number of samples acquired n covers exactly a complete period [34], and then we have Substituting Equations (20) and (23) into Equation (22), the azimuth uncertainty σ ψω sg is In Equation (24), ω e = 15.041 • /h and ϕ is a constant in specified latitude (for instance, ϕ = 43.8 • in Changchun, China). As shown in Figure 6a (arc second is used to denote the unit of azimuth uncertainty, and 1 arc second = 1" = 1/3600 • ), increasing the number of measurement points n and reducing the measurement uncertainty σ ω will lower the azimuth uncertainty. However, it is noted that azimuth uncertainty decreases slowly when the σ ω is quite small (see Figure 6b). In particular, when σ ω = 0.005 • /h, the azimuth uncertainty decreases from 15 arc second (n = 70) to 10 arc second (n = 240); when σ ω = 0.001 • /h, the azimuth uncertainty decreases from 3 arc second (n = 70) to 2 arc second (n = 240). Therefore, a proper number of measurement points should be chosen according to the required azimuth precision and the measurement time. In addition, the latitude affects the azimuth uncertainty too. Based on Equation (24), the azimuth uncertainty gradually decreases as the latitude decreases from high to low.  Subsequently, the encoder measurement uncertainty is analyzed. According to Equation (19), the azimuth uncertainty caused by the encoder measurement uncertainty σ γ at position i is Substituting Equation (23) into Equation (25), in one circle measurement the overall azimuth uncertainty caused by the encoder measurement uncertainty with η = ε = 0 is According to Equation (26), the encoder uncertainty σ γ has a significant effect on the azimuth results, and the encoder uncertainty will propagate to the azimuth uncertainty directly. Therefore, a high precision encoder is required for a high precision north finder.
Next, the gyro misalignment's effect is analyzed. According to Equation (19), gyro outputs with η = 0 equal ω sgi = ω e cos ϕ cos(γ i + ψ) cos ε + ω e sin ϕ sin ε Here ε is a constant for the assembled gyro, ω e = 15.041 • /h, and ϕ is constant in a specified latitude. Thus Equation (27) can be rewritten as ω sgi = C 1 cos(γ i + ψ) + C 2 C 1 =ω e cos ϕ cos ε, C 2 = ω e sin ϕ sin ε (28) Here C 1 and C 2 are constants, hence ε has no effect on the azimuth uncertainty when η = 0. Finally, the shaft misalignment's effect on the azimuth results is analyzed. The azimuth uncertainty caused by the shaft misalignment is Substituting Equation (19) into Equation (29), σ ψη is obtained as follows, where the symbols a 1 , a 2 , a 3 , a 4 , a 5 , a 6 , b 1 , b 2 are introduced to simplify the equation: n cos ϕ 2 (b 1 +b 2 ) a 1 = cos η 2 sin ε 2 (0.125 + 0.25 cos ψ 2 ) a 2 = sin η 2 cos ε 2 (0.0625 + 0.25 cos ψ 2 ) a 3 = sin η 2 cos η 2 (0.375 sin ε 4 + 0.2734 cos ε 4 ) a 4 = sin η 2 cos η(0.0782 cos ε 4 − 0.125 sin ε 2 cos ε 2 ) a 5 = 0.0625 cos ε 2 cos η 2 sin ε 2 + 0.0234 cos ε 4 sin η 2 + 0.3125 sin ε 2 cos ε 2 a 6 = 1 − (0.375 cos ε 2 + 0.25 cos ε 2 cos η + 0.5 sin η 2 sin ε 2 + 0.375 cos η 2 cos ε 2 ) b 1 = sin η 2 sin ε 2 (0.125 + 0.25 sin ψ 2 ) + cos η 2 cos ε 2 (0.0625 + 0.25 sin ψ 2 ) b 2 = 0.25 cos ψ 2 cos ε 2 + 0.0625 cos ε 2 + 0.125 cos ε 2 cos η (30) According to Equation (30), the relationship of the azimuth uncertainty against the shaft misalignment is shown in Figure 7. It is clear from Figure 7a that the azimuth uncertainty σ ψη increases as the shaft misalignment angle η grows; however, as the gyro misalignment angle ε increases, the azimuth uncertainty slightly falls. When the shaft misalignment angle η equals zero, the azimuth uncertainty remains at zero regardless of the value of ε. Differently, the azimuth value ψ shows no effect on azimuth uncertainty (see Figure 7b). Usually, the gyro misalignment angle ε is very small and measuring it is troublesome. Therefore, a much better choice for achieving a lower azimuth uncertainty is to reduce the shaft misalignment angle. Based on the theory of propagation of uncertainty, the total azimuth uncertainty σ ψ is expressed as: According to Equation (31), in order to improve the azimuth measurement precision, the following procedures should be conducted: 1. Constrain the gyro bias error and the gyro random error. Regarding to the former, on one hand, selecting a high precision gyro usually attains a low bias error; on the other hand, using filtering approaches [20,24,32] is an alternative method to reduce the gyro bias error. For the gyro random error, a simple and effective method is to average a large number of gyro outputs at each rotation position. Other filtering approaches can be found in [21,22,24,25] . 2. Increase the number of measurement points (ensure that the north finding time is less than the time scale of the bias instability [4]). 3. Execute the azimuth measurement at a low latitude area. When the north finding experiment is executed at the equator (ϕ = 0 • ), all the earth rotation rate is detected by the gyro sensitive axis. In contrast, the north finding task is impossible to achieve when ϕ = 90 • . Generally, the north finding experiment should be conducted away from the Antarctic Circle and the Arctic Circle. 4. Improve the encoder measurement precision. In the multi-position north finding scheme, a high precise rotary platform and a corresponding high precision encoder are prerequisites to achieve high precision azimuth results. 5. Align the shaft to make the shaft misalignment angle as small as possible. As described in Equation (28), the gyro misalignment angle has no effect on the azimuth results when shaft misalignment angle is 0. However, when a shaft misalignment exists, both the shaft misalignment angle and the gyro misalignment angle corrupt the gyro outputs, as shown in Figure 7a. Therefore, the best choice is to align the shaft.
The criteriafor the selection of the system components, including the rate gyro, the inclinometer and the encoder, are listed as follows: 1. Building a north finding system with azimuth uncertainty of σ ψ (Equation (31)) is the goal. Through rough calculation (distribute azimuth uncertainty), the value of σ ψω sg , σ ψγ and σ ψη (Equation (31)) are determined. 2. Based on the value of σ ψω sg and Equation (24), the gyro bias error and the gyro random error are decided, and then the specified rate gyro can be selected. 3. Combining the value of σ ψγ and Equation (26), the encoder uncertainty is calculated, thus the encoder is selected. 4. According the value of σ ψη and Equation (29), the shaft misalignment angle is calculated, thus the shaft misalignment angle measurement instrument can be chosen. In practice, we usually align this angle to a very small value. The details of shaft misalignment angle calculation and alignment process refer to [35].

Implementation of a High Precision Robust Gyro North Finder
To evaluate the proposed north finding method, a north finder setup is built as shown in Figure 8. Currently, depending on the performance specifications of the gyros, MEMS gyros are used in consumer area, ring laser gyros are used in tactical area, resonant optical gyros, fiber optic gyros and dynamically tuned gyros are used in navigational area [5]. In order to achieve high precision azimuth measurement, a dynamically tuned gyro in navigational grade is used in our setup. The dynamically tuned gyro in this setup has a bias error of 0.003 • /h. Through filtering, the gyro random error is found to be 0.002 • /h. Hence, the gyro measurement uncertainty is 0.005 • /h. At the same latitude (ϕ = 43.8 • ), all experiments are carried out using the scheme in Figure 2. At each position i, the gyro stays for 2 s and is then automatically rotated counterclockwise by 360/n • to position i + 1 by the motor.

Gyro Outputs Modeling with Shaft Misalignment
As described in Figure 7, the shaft misalignment will corrupt the gyro outputs. In order to achieve high precision azimuth result, it is essential to measure the shaft misalignment angle η between the practical shaft OR 1 and the vertical shaft OR, and to adjust the practical shaft to a vertical state. The shaft misalignment angle measurement procedure includes two steps using the inclinometer in our system [35]. Firstly, the inclinometer outputs are obtained at position 1, and θ X1 , θ Y1 indicate the inclinometer X+ axis and the Y+ axis outputs respectively. Then the inclinometer is rotated counterclockwise by 180 • around the shaft to position i, and θ X2 , θ Y2 indicate the inclinometer X+ axis and the Y+ axis outputs respectively. Generally, (θ X2 − θ X1 )/2 and (θ Y2 − θ Y1 )/2 are small angles (less than 1 • ), and then the relationship of η against the inclinometer outputs is The shaft misalignment angle is adjusted by regulating the three knobs. η = 0 is acquired when The resolution of the inclinometer in our experiment is 0.0005 • (NS-5/P2), and the standard deviation of the shaft alignment angle reaches at 0.001 • . Using the approaches in [35], a shaft alignment precision of 0.003 • is achieved. Two comparison experiments are conducted to evaluate the shaft misalignment's effect on azimuth measurement results. The first experiment is conducted at η = 1 • , and the second one is operated at η = 0 • . The rotation points n are 72 in both experiments and their gyro outputs are shown in Figure 9. Ideally, the curve of the gyro outputs from position 1 to position n is a standard sinusoidal. However, when the shaft is misaligned, the practical gyro outputs deviates from the ideal gyro outputs. As shown in Figure 9, the gyro outputs of η = 1 • ) is bigger than those of η = 0 • in the troughs; in contrast, the gyro outputs of η = 1 • ) is smaller than those of η = 0 • in the peaks.
The function of ω sgi = A cos γ i + B sin γ i + C is used to fit the gyro outputs. By subtracting the fitting results and the gyro outputs, the gyro outputs residual of η = 0 • and η = 1 • are obtained as shown in Figure 10. Naturally, the residual is randomly located around the line of x = 0. However, it is clear from Figure 10 that the locations of the residual of η = 1 • does not fit this pattern. This is due to the fact that the gyro outputs of η = 1 • deviate from the sinusoidal gyro outputs using sine wave fitting. Therefore, it is strongly recommended to align the shaft to vertical state before north finding. The following experiments are all conducted in shaft aligned state.

Azimuth Measurement Results
Increasing the number of measurement points in one circle measurement lowers the azimuth uncertainty but takes more time. This is a trade-off of azimuth uncertainty and measurement time. Therefore, a proper number of measurement points should be chosen according to Equation (31) to attain a required azimuth precision. In our north finder system, σ γ = 0.001 • and η = 0.003 • . Based on Equations (26) and (29), σ ψγ = 3.6" and σ ψη = 0.7". Substituting ω e , σ ω and ϕ by 15.041 • , 0.005 • and 43.8 • respectively, and specifying n into Equation (24), δ ψωe4s is achieved. Finally, the overall azimuth uncertainty δ ψ is calculated using Equation (31). Here, two criteria assist us to select the number of measurement points. The first one is to calculate the specified n from the given required azimuth uncertainty according to Equation (31). The second one is to ensure that the number of measurement points n large enough to attain desirable fitting results. Based on these two criteria, the minimum value of n is set at 72 in our system, and the north finding results of n = 72, n = 90, n = 120 and n = 180 are given.
When n is set at 72, 10 repeated azimuth measurement experiments are successively conducted based on the working process in Figure 2. For each run, the gyro outputs and the corresponding encoder outputs at the 72 positions are recorded. ω sgi = A cos γ i + B sin γ i + C is then used to fit the gyro outputs. The ten fitting results are given in Table 1. Averaging the 10 runs' azimuth results in Table 1, the average azimuth is 65.5378 • . The standard deviation is 33". In each run, the maximum of the sum of residual square is less than 0.15 × 0.001 ( • /h) 2 . The first run's gyro outputs is plotted in Figure 11. As expected, the gyro maximum output is 10.83 • /h, pointing north; the minimum output is −10.83 • /h, pointing south. The gyro outputs fit the sine wave quite well (see Figure 11a), and the maximum residual is less than 0.04 • /h.  Table 2 gives the 10 groups of azimuth measurement results when n is set as 90. According to Table 2, the average azimuth is 65.5467 • , and the standard deviation is 25 arc second. The maximum of the sum of residual square in each run is less than 0.9 × 0.0001( • /h) 2 , which is 0.6 × 0.0001 ( • /h) 2 smaller than the former 10 runs when n is set at 72. Similarly, the first run's experiment data in Table 2 are chosen to illustrate the azimuth fitting results, as shown in Figure 12. The maximum residual is less than 0.03 • /h, about 25% lower than that of the first 10 runs. Thirdly, setting n = 120, Table 3 gives the 10 groups of azimuth measurement results. As listed in Table 3, the average azimuth is 65.5504 • , and the standard deviation is 22". The maximum of the sum of residual square in each run is less than 0.9 × 0.0001 ( • /h) 2 . Figure 13 shows the fitting results of the 1-th group experiment data in Table 3. The maximum residual is less than 0.025 • /h, which is 37.5% lower than the first ten 10 runs. Finally, Table 4 gives the 10 groups of azimuth measurement results when n = 180. According to Table 4, the average azimuth is 65.5445 • , and the standard deviation is 18 arc second. The maximum of the sum of residual square in each run is less than 0.82 × 0.0001( • /h) 2 . The first group experiment data in Table 4 are chosen to illustrate the azimuth fitting process, as shown in Figure 14. The maximum residual is less than 0.023 • /h, about 40% smaller than that of the first 10 runs. Based on Equation (31), the simulated azimuth uncertainties are 19", 17", 15" and 13" when n are 72, 90, 120 and 180 respectively. As shown in Tables 1-4, the experimental azimuth uncertainties are 33", 25", 22" and 18" arc second accordingly. The differences between the azimuth uncertainties and the simulated azimuth uncertainties are 14" (n = 72), 8"(n = 90), 7"(n = 120) and 5"(n = 180). Clearly, the experimental azimuth uncertainties are closer and closer to the simulated azimuth uncertainties as n goes up. Equation (31) is obtained according to Equation (19), and Equation (19) is a continuous function. However, in our experiment, the number of measurement points n is a limited number, i.e., the gyro outputs are digitized. Hence, the larger the value of n is, the better the fitting results are. In our experiment, the rotation time from position i to position i + 1 is about 0.2 s, and the sampling time at each measurement position is 2 s, so 2.2 s are taken at each measurement. Therefore, the measuring time are about 158 s, 198 s, 264 s, 396 s when n are 72, 90, 120 and 180 respectively. Increasing the number of measurement points leads to higher precision but longer measuring time.

Conclusions
In this paper we propose a multi-position non-continuous rotation gyro north finding method to measure the azimuth between the observer and north direction. We proposed a mathematical model for gyro outputs at each rotation position. Our model concerns the key factors related to the azimuth results, including the gyro bias error, the gyro random error, scale-factor error, encoder precision, the number of measurement points, the gyro misalignment angle, and the shaft misalignment angle.
Combining the gyro outputs model and the theory of propagation of uncertainty, the optimized azimuth measurement scheme is achieved. According to this, the corresponding approaches to attain high precision azimuth results are provided. Finally, the azimuth uncertainty of 18" (equals to 0.087 mrad) is obtained in our system. Currently, the azimuth uncertainty of the reported north finding system [4] and the product on the market are single digit mrad. It is noted that the specifications of the gyros has a significant effect on the azimuth uncertainty according to Equation (31), therefore, high performance gyros, e.g., navigational grade gyros, achieve better azimuth measurement results than consumer grade gyros. This paper provides systematic theory for analyzing the details of the gyro north finder scheme from simulation to implementation. It is useful and valuable to both applied researchers in academia and advanced practitioners in industry.