MITIGATING SYSTEMATIC ERROR IN GYROCOMPASSING

Abstract

A MEMS gyrocompass and method are provided to mitigate systematic error in determination of a north angle. The MEMS gyrocompass includes one or more MEMS gyroscopes having a sense axis within a reference plane. Samples from an output of the MEMS gyroscope are obtained in at least two angles of rotation about an axis perpendicular to the reference plane. First fit coefficients are determined by fitting samples with first fitting functions determined as function of time. Second fit coefficients are determined by fitting components of earth rotation rate projected on the reference plane based on samples obtained by all the MEMS gyroscopes, which fitting is performed with a second fitting function determined as a function of rotation angle of the MEMS gyroscope with respect to a reference angle. The north angle is determined as an angle between the reference angle and true north based on the second fit coefficients.

Claims

1. A method for mitigating systematic error in north angle determined by a MEMS gyrocompass that includes at least one MEMS gyroscope having a sense axis within a reference plane, the method comprising: obtaining samples from an output of the at least one MEMS gyroscope in at least two different angles of rotation about an axis perpendicular to the reference plane; fitting an Earth rate response of the output as a function of a rotation angle simultaneously with individually fitting rate drift of each MEMS gyroscope of the at least one MEMS gyroscope as a function of time, wherein the fittings is performed by: determining one or more first fit coefficients by fitting samples obtained from each individual MEMS gyroscope with first fitting functions determined as the function of time, and determining second fit coefficients by fitting components of an Earth rotation rate projected on the reference plane based on samples obtained from all of the at least one MEMS gyroscope with a second fitting function determined as the function of rotation angle of the at least one MEMS gyroscope with respect to a reference angle; and determining the north angle as an angle between the reference angle and true north based on the second fit coefficients.

2. The method according to claim 1, wherein the first fitting function is a rate drift fitting function.

3. The method according to claim 2, wherein the rate drift fitting function is one of a first order function, a second order function, an exponential function and a sum of exponential functions.

4. The method according to claim 3, further comprising determining the first fit coefficients individually to each MEMS gyroscope.

5. The method according to claim 1, wherein the second fit coefficients are a cosine coefficient (c) and a sine coefficient (s), and wherein the north angle is calculated by arctan (s/c).

6. The method according to claim 1, wherein the at least one MEMS gyroscope comprises two or more MEMS gyroscopes that are arranged such that a vector sum of sense axes of the respective two or more MEMS gyroscopes is zero.

7. The method according to claim 1, wherein the determining of the first fit coefficients and the second fit coefficients are performed simultaneously.

8. A MEMS gyrocompass for mitigating systematic error in determination of a north angle, the MEMS gyrocompass comprising: at least one MEMS gyroscope having a sense axis within a reference plane; a memory; and a processing device communicatively coupled to the memory and that is configured to execute code on the memory that configures the processing device to: obtain samples from an output of the at least one MEMS gyroscope in at least two different angles of rotation about an axis perpendicular to the reference plane; fit an Earth rate response of the output as a function of a rotation angle simultaneously with individually fitting rate drift of each MEMS gyroscope of the at least one MEMS gyroscope as a function of time, wherein the processing device further determines the fittings by: determining one or more first fit coefficients by fitting samples obtained from each individual MEMS gyroscope with first fitting functions determined as the function of time, and determining second fit coefficients by fitting components of an Earth rotation rate projected on the reference plane based on samples obtained from all of the at least one MEMS gyroscope with a second fitting function determined as the function of rotation angle of the at least one MEMS gyroscope with respect to a reference angle; and wherein the processing device of the MEMS gyrocompass is further configured to determine the north angle as an angle between the reference angle and true north based on the second fit coefficients.

9. The MEMS gyrocompass according to claim 8, wherein the first fitting function is a rate drift fitting function.

10. The MEMS gyrocompass according to claim 9, wherein the rate drift fitting function is one of a first order function, a second order function, an exponential function and a sum of exponential functions.

11. The MEMS gyrocompass according to claim 10, wherein the processing device is further configured to determine the first fit coefficients individually to each MEMS gyroscope.

12. The MEMS gyrocompass according to claim 8, wherein the second fit coefficients are a cosine coefficient (c) and a sine coefficient (s), and wherein the north angle is calculated by arctan (s/c).

13. The MEMS gyrocompass according to claim 8, wherein the at least one MEMS gyroscope comprises two or more MEMS gyroscopes that are arranged such that a vector sum of sense axes of the respective two or more MEMS gyroscopes is zero.

14. The MEMS gyrocompass according to claim 8, wherein the processing device is further configured to determine the first fit coefficients and the second fit coefficients simultaneously.

15. A computer program product storing code with instructions executable by a processing device for mitigating systematic error in north angle determined by a MEMS gyrocompass that includes at least one MEMS gyroscope having a sense axis within a reference plane, comprising: code for obtaining samples from an output of the at least one MEMS gyroscope in at least two different angles of rotation about an axis perpendicular to the reference plane; code for fitting an Earth rate response of the output as a function of a rotation angle simultaneously with individually fitting rate drift of each MEMS gyroscope of the at least one MEMS gyroscope as a function of time, wherein the fittings is performed by: determining one or more first fit coefficients by fitting samples obtained from each individual MEMS gyroscope with first fitting functions determined as the function of time, and determining second fit coefficients by fitting components of an Earth rotation rate projected on the reference plane based on samples obtained from all of the at least one MEMS gyroscope with a second fitting function determined as the function of rotation angle of the at least one MEMS gyroscope with respect to a reference angle; and code for determining the north angle as an angle between the reference angle and true north based on the second fit coefficients.

16. The computer program product according to claim 15, wherein the first fitting function is a rate drift fitting function, which is one of a first order function, a second order function, an exponential function and a sum of exponential functions.

17. The computer program product according to claim 16, further comprising code for determining the first fit coefficients individually to each MEMS gyroscope.

18. The computer program product according to claim 16, wherein the second fit coefficients are a cosine coefficient (c) and a sine coefficient (s), and wherein the north angle is calculated by arctan (s/c).

19. The computer program product according to claim 16, wherein the at least one MEMS gyroscope comprises two or more MEMS gyroscopes that are arranged such that a vector sum of sense axes of the respective two or more MEMS gyroscopes is zero.

20. The computer program product according to claim 16, further comprising code for the determining the first fit coefficients and the second fit coefficients simultaneously.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0038] In the following the invention will be described in greater detail, in connection with exemplary embodiments, with reference to the attached drawings, in which:

[0039] FIG. 1 is an Allan deviation plot of a gyroscope.

[0040] FIG. 2 is another Allan deviation plot of a gyroscope.

[0041] FIG. 3 illustrates angle error of one gyroscope as a function of north angle.

[0042] FIG. 4 illustrates angle error of one gyroscope as a function of north angle.

[0043] FIG. 5 illustrates north angle error received from two gyroscopes.

[0044] FIG. 6 illustrates stepwise rotation of a turn table.

[0045] FIG. 7 illustrates linear rate drift.

[0046] FIG. 8 illustrates quadratic rate drift.

[0047] FIG. 9A illustrates step-wise rotation of a turntable.

[0048] FIG. 9B illustrates output of a single gyroscope without rate drift .

[0049] FIG. 9C illustrates gyroscope output as a function of time when no rate drift is present.

[0050] FIG. 10A illustrates step-wise rotation of a turntable.

[0051] FIG. 10B illustrates different output values of a gyroscope caused by rate drift.

[0052] FIG. 10C illustrates gyroscope output with systematic error caused by lack of rate drift fitting.

[0053] FIG. 11A illustrates step-wise rotation of a turntable.

[0054] FIG. 11B illustrates different output values of a gyroscope caused by rate drift.

[0055] FIG. 11C illustrates gyroscope output with systematic error correction by rate drift fitting.

[0056] FIG. 12 illustrates a computer simulation platform with eight single-axis gyroscopes.

[0057] FIG. 13 is a schematic presentation of functional elements of a MEMS gyrocompass.

DETAILED DESCRIPTION

[0058] FIG. 6 illustrates an exemplary stepwise rotation of the turn table that is configured to perform a carouselling operation for a gyroscope. The turn table performs 30 steps of 24? over a duration of one minute (60 s), thus going two full turns (i.e., 720?). It is noted that although technical effects of the invention have been simulated using a turntable setup, the exemplary aspects described herein are applicable with any kind of rotation schemes, which do not have to be known in advance. Only assumption is that rotation occurs with respect to an axis that is orthogonal to the reference plane such that at least two samples are taken at different rotation angles. For example, rotation may be implemented as slow turning of a vehicle in which the MEMS gyrocompass is installed, or the MEMS gyrocompass may be rotated by manually rotating a handheld device comprising the MEMS gyrocompass by an arbitrary, non-zero amount. Moreover, it is noted that although the exemplary aspects of the invention utilize measurement of rotation an axis perpendicular to the reference plane, it is understood that rotation motion may comprise components also deviating from rotation in the reference plane. Effects of such non-useful additional rotation components are omitted by the algorithm described herein.

[0059] The FIG. 7 illustrates linear rate drift (41) of a single gyroscope. Linear rate drift can be approximated with a first-order equation. When the turn table is rotated, output samples (40) of the gyroscope change at every angle. This can be seen as a sine-like pattern of samples (40) taken of the output signal, here shown over two full turns of the turn table carouselling scheme illustrated in the FIG. 6. Without rate drift, the zero level, minimum and maximum of the intermittent sine pattern formed by the samples (40) would not change over time. However, in this example, there is linear rate drift in phase, which causes the output levels into a constant change.

[0060] In this case, rate drift (41) of the gyroscope output can be fitted on a straight line that can be represented with a first-order polynomic function.

[0061] FIG. 7 illustrates quadratic rate drift (41) of a single gyroscope, which can be represented with a second-order polynomic function.

[0062] Rate drift can also be exponential, in other words it can be best approximated with an exponential function. This may be a single exponential function or a sum of exponential functions. For example, there may be two or more different exponentially behaving rate drifts occurring in the gyroscope, which are characterized with different time constants. With short time intervals, linear fitting likely provides best results. Measured angular rates can in some cases even have no rate drift, even though some bias error still affects the output.

[0063] Exemplary rate drift functions can be formulated as follows:

linear: a.sub.0+a.sub.1t

exponential (1): a.sub.0+?.sub.1exp(?t/?.sub.1)

exponential (2): a.sub.0+a.sub.1exp(?t/?.sub.1)+a.sub.2exp(?t/?.sub.2)

quadratic: a.sub.0+a.sub.1t+a.sub.2t.sup.2

[0064] To avoid overfitting, simple equations such as linear and exponential (1) with single time constant ?.sub.1 are preferred over the more complex quadratic and exponential (2) with two time constants ?.sub.1 and ?.sub.2.

[0065] In general, it is noted that the basic algorithm for removing systematic error from output of a gyroscope array with n gyroscopes can be expressed mathematically as a matrix equation. When a polynomic fitting function is used, size of the matrix A has size s?s, where s=2+n(o+1), n refers to number of gyroscopes in the gyroscope array, o is order of rate drift function: o=2 refers to quadratic rate drift function, o=1 refers to linear rate drift function, 0=0 is a constant (no drift). This basic rule of size of the matrix for polynomic fitting is not applicable with exponential fitting functions. Vector y has length s.

[0066] An example of the matrix equation applied for a gyrocompass with two gyroscopes and quadratic rate drift function, in other words a rate drift function of an individual gyroscope follows a second order polynomial equation having formula a.sub.0,i+a.sub.1,it+a.sub.2,it.sup.2, can be represented as:

[00001]$\begin{array}{cc}\mathrm{Ay}=\left(\begin{array}{cccccccc}.Math.& .Math.& .Math.& .Math.& .Math.& .Math.& .Math.& .Math.\\ \cos \left({?}_i+{?}_1\right)& \sin \left({?}_i+{?}_1\right)& 1& t_i& t_i^2& 0& 0& 0\\ \cos \left({?}_j+{?}_2\right)& \sin \left({?}_j+{?}_2\right)& 0& 0& 0& 1& t_j& t_j^2\\ .Math.& .Math.& .Math.& .Math.& .Math.& .Math.& .Math.& .Math.\end{array}\right)\left(\begin{array}{c}c\\ s\\ a_{0,1}\\ a_{1,1}\\ a_{2,1}\\ a_{0,2}\\ a_{1,2}\\ a_{2,2}\end{array}\right)=\left(\begin{array}{c}.Math.\\ y_{1,i}\\ y_{2,j}\\ .Math.\end{array}\right)& \left(2\right)\end{array}$

[0067] In this equation, t.sub.i refers to a time instance, ?.sub.1, ?.sub.2, refer to sense axis angles of gyroscopes on the board with respect to a reference axis within the reference plane, ?.sub.i=?(t.sub.i) refers to board orientation at time instance t.sub.i. C and s refer to cosine and sine coefficients that are used to determine the north angle as an angle between the reference axis and the true north, and a.sub.0,1, a.sub.1,1, a.sub.2,1 refer to rate drift coefficients of gyroscope 1. In this example, rate drift function is quadratic, in other words rate drift function of gyroscope 1 follows an equation having formula a.sub.0,1+a.sub.1,1t+a.sub.2,1t.sup.2. y.sub.1,i refers to gyroscope 1 output at time instance i and y.sub.2,j refers to gyroscope 2 output at time instance j.

[0068] Rate drift fit coefficients a.sub.0,1, a.sub.1,1, a.sub.2,1 are obtained from x=A.sup.?1y.

[0069] The north angle is calculated from the cosine and sine coefficients c and s by equation:

[00002]$\begin{array}{cc}\arctan \frac{s}{c}& \left(3\right)\end{array}$

[0070] In-plane component of the Earth's rotation rate in the reference plane is calculated as:

?{square root over (s.sup.2+c.sup.2)} (4)

[0071] An example of the matrix equation applied for a gyrocompass with two gyroscopes and exponential rate drift function, in other words a rate drift function of an individual gyroscope follows an exponential equation having formula a.sub.0,k+a.sub.1,ke.sup.?.sup.k.sup.t, can be represented as follows:

[00003]$\begin{array}{cc}\left(\begin{array}{ccccccc}\cos {?}_i& \sin {?}_i& 1& t_i& {?}_0^{t_i}\cos ?\left(t\right)\mathrm{dt}& {?}_0^{t_i}\sin ?\left(t\right)\mathrm{dt}& {?}_0^{t_i}{y}_k\left(t\right)\mathrm{dt}\\ \cos {?}_j& \sin {?}_i& 1& t_j& {?}_0^{t_j}\cos ?\left(t\right)\mathrm{dt}& {?}_0^{t_j}\sin ?\left(t\right)\mathrm{dt}& {?}_0^{t_j}{y}_k\left(t\right)\mathrm{dt}\\ .Math.& .Math.& .Math.& .Math.& .Math.& .Math.& .Math.\end{array}\right)\left(\begin{array}{c}{d}_1\\ d_2\\ d_3\\ d_4\\ d_5\\ d_6\\ {?}_k\end{array}\right)=\left(\begin{array}{c}{y}_{k,i}\\ y_{k,j}\\ .Math.\end{array}\right)& \left(5\right)\end{array}$

wherein t.sub.i is time instance, ?.sub.k is orientation of the kth gyro on the board, ?(t.sub.i) is board orientation at time instance t.sub.i, y.sub.k,i is gyroscope k output at time instance i. ? . . . dt is a time integral, evaluated numerically using for example trapezoidal rule or Simpson's rule if the time steps are equal.

[0072] Fit coefficients ?.sub.k related to time constants ?.sub.k=?1/?.sub.k are determined using the matrix equation (6) below for every k gyroscope axis, solving it in the least-squares sense. This can be performed by writing equation (5) as Ax=b, least-squares sense means to minimize ?Ax?b?.sub.2, where ?.Math.?.sub.2 is a vector norm. The fit coefficient ?.sub.k is utilized in further analysis, whereas other coefficients d.sub.1, d.sub.2, . . . , d.sub.6 are not.

[00004]$\begin{array}{cc}\left(\begin{array}{ccccccc}.Math.& .Math.& .Math.& .Math.& .Math.& .Math.& \\ \cos \left({?}_i+{?}_1\right)& \sin \left({?}_i+{?}_1\right)& 1& e^{{?}_1{t}_i}& 0& 0& .Math.\\ \cos \left({?}_j+{?}_2\right)& \sin \left({?}_j+{?}_2\right)& 0& 0& 1& e^{{?}_2{t}_j}& .Math.\\ .Math.& .Math.& .Math.& .Math.& .Math.& .Math.& \end{array}\right)\left(\begin{array}{c}c\\ s\\ a_{0,1}\\ a_{1,1}\\ a_{0,2}\\ a_{1,2}\end{array}\right)=\left(\begin{array}{c}.Math.\\ y_{1,i}\\ y_{2,j}\\ .Math.\end{array}\right)& \left(6\right)\end{array}$

[0073] Cosine coefficient c and sine coefficient s are solved in least-squares sense, and north angle and in-plane component of the Earth rate can be calculated as given above in equations (3) and (4).

[0074] FIGS. 9A to 11C illustrate effects of the invented algorithm on accuracy of measurements.

[0075] In particular, FIG. 9A illustrates another example of stepwise rotation of the turn table. The turn table performs 4 steps of 90? over a duration of twenty seconds (20 s), thus going three full turns (1080?) in one minute.

[0076] The FIG. 9B illustrates output of a single gyroscope rotated by the turn table as a function angle. In this example, there is no rate drift nor rate drift fitting. Output values (62) of the gyroscope received at the same rotation angle but at different times, in other words during different rotation rounds of the turn table, are the same. In this non-limiting example, each dot representing an output value (62) represents three different measurements. The output values (62) can be easily fitted on a sine curve (60). No systematic error is present and north angle can be accurately determined.

[0077] Finally, the FIG. 9C illustrates gyroscope output as a function of time, which gives a good view on the total amplitude of the output in mdps (milli-degrees per second).

[0078] Tested device in the case 9A to 9C gives a north angle or 31.234 degrees with signal amplitude of 4.178 mdps, which corresponds with theoretic Earth's rotation rate with high accuracy. These values represent the correct reference values for further tests.

[0079] FIGS. 10A to 10C illustrate the same device as in the FIGS. 9A to 9C. This time there is linear rate drift in place, but no rate drift fitting in time dimension taking place and thus no systematic error correction. For illustration purposes, a strong linear rate drift was applied.

[0080] In particular, FIG. 10A shows the stepwise rotation of the turn table, which is the same as shown in the FIG. 9A.

[0081] FIGS. 10B and 10C show effects of rate drift in the output of the gyroscope. Because of rate drift, output values (62) of the gyroscope, received at the same rotation angle but at different times, in other words during different rotation rounds of the turn table, are different.

[0082] Instead of having a single dot covering all three measurements, the output shown in the FIG. 10B shows three different output values for each turn table angle. In the FIG. 10C output values show three groups on three different levels.

[0083] Traditionally fitting is performed by fitting a sine or cosine function to the received rotation rates as shown in the FIG. 10B. The sine or cosine function may be of the form C*sin(?.sub.j+?) or C*cos(?.sub.j+?) respectively, where C is the amplitude of the sine or cosine function, which depends at least on the latitude of the device, ?.sub.j is the offset angle of the jth gyroscope's sense axis and ? is a phase offset, which has a fixed difference with respect to true north, and the heading of the device may be determined by varying the phase offset ? to find the value with the minimum error in the fit of the sine or cosine function to the received rotation rates. The sine or cosine function may be fit to the received rotation rates using a least-squares mean fitting method or other fitting method.

[0084] When output values are averaged over multiple samples (in this example, three samples taken on three intermittently repeated similar rotation table positions), the fitted sine or cosine function as function of angle (60) is way above the zero level. Consequently, the fitted curve as function of time (61) is offset. Using these gyroscope output values in north angle and rotation rate calculations causes a significant systematic error in the result. In this example, with heavy rate drift used for visualization, the north angle value received was 17.77 degrees (compared to the correct 31.234 degrees), thus with significant error of 13.46 degrees, and output amplitude was 2.169 mdps, which has a significant 2 mdps error from the correct value of 4.178 mdps, almost 50%. Such error magnitudes are not acceptable.

[0085] FIGS. 11A to 11C illustrate the same device as in the FIGS. 9A to 10C, but this time there is linear rate drift. Difference to previously disclosed examples is that rate drift fitting is taking place and thus correction of systematic error is implemented. A linear, first order rate fitting function was applied in this non-limiting example.

[0086] In particular, FIG. 11A shows the stepwise rotation of the turn table, which is the same as shown in the FIGS. 9A and 10A.

[0087] With proper rate drift fitting taking place, systematic error is removed, and the gyroscope output is brought back zero level as expected.

[0088] It should be appreciated that the benefit of performing rate drift fitting as function of time becomes obvious from the FIGS. 11B and 11C. When the linear rate drift of the gyroscope is fitted with a linear rate drift function in time domain, rate drift fitting fully removes effect of rate drift and thus correct output can be obtained. As a result, resulting north angle calculated based on the fitted output as function of time gives a north angle of 31.234 degrees and amplitude (Earth's rotation rate) of 4.169 mdps, both equal to the values received in the non-drift example shown in the FIGS. 9A to 9C.

[0089] The angles of the sense axes of the one or more MEMS gyroscopes are preferably angled at regular intervals of 360?/N, where N is the number of sense axes. In case of using single-axis gyroscopes, N is also number of gyroscopes. By applying regular intervals of 360?/N between the sense axes of a plurality of mutually similar gyroscopes, vector sum of sense axes of the gyroscopes is zero.

[0090] N may be an even number such that the sense axes of the one or more MEMS gyroscopes are arranged into pairs, wherein the sense axes of each pair of gyroscopes are offset by 180 degrees, thus pairwise zeroing the vector sums of sense axes. N may be an odd number. For example, three gyroscopes with 120? angle between their sense axes produces a zero vector sum.

[0091] The algorithm was simulated with a computer simulation platform with eight single-axis gyroscopes (10-1, . . . , 10-8), placed at sense axis angles ?.sub.i={0?, 180?, 45?, 225?, 90?, 270?, 135?, 315?} on a rotatable board (12) as schematically illustrated in the FIG. 12. However, computer simulations indicate that three or four gyroscopes and three or four axes already provide good results, while adding more than four gyroscopes and axis provides less significant statistical 1/?N-improvement. Sense axis of the gyroscopes (10) are illustrated with the direction of the arrow. Sense axis angles ?.sub.j are determined in relation to a reference axis (15). Only sense axis angles ?.sub.2=180 and ?.sub.3=45? are shown in the FIG. 12 for clarity. The board (12) was rotated by a turn table that was rotated stepwise to implement carouselling, as illustrated in the FIG. 6. Total sampling time used in the test was 60 seconds, with two full turns (720 degrees). The same algorithm is also applicable with continuous carouselling, with samples taken intermittently at regular time intervals. Although the example in FIG. 12 illustrates gyroscopes (10) arranged radially in a circular pattern, any placement of gyroscopes on the plane of the board is applicable as soon as sense axis angles with respect to a reference axis within the plane of the board (12) remain. Number of gyroscopes may also vary. For example, a single gyroscope may have two orthogonal sense axes, in which case just one gyroscope is needed for sensing two different sense axes. Results of computer simulations are shown in the table 2 below.

TABLE-US-00002 TABLE 2 Fit time Time dependency dependency mean std # ARW *100e?6 dps/s order (deg) (deg) 1 0 Linear 0 0.000 0.00 [1, 1, 1, 1, 1, 1, 1, 1] 2 0 Linear 0 1.38 0.00 [1, 0, 0, 0, 0, 0, 0, 0] 3 0 Linear 0 2.23 0.00 [1, 0, ?1, 0, 0, 0, 0, 0] 4 0 Linear 0 3.59 0.00 [1, 0, ?1, 0, ?1, 0, 0, 0] 5 0 Linear 1 0.000 0.00 [1, 0, ?1, 0, ?1, 0, 0, 0] 6 0 Quadratic 0 0.000 0.00 [1, 1, 1, 1, 1, 1, 1, 1] 7 0 Quadratic 0 0.75 0.00 [1, 0, 0, 0, 0, 0, 0, 0] 8 0 Quadratic 1 ?0.084 0.00 [1, 0, 0, 0, 0, 0, 0, 0] 9 0 Quadratic 2 0.000 0.00 [1, 0, 0, 0, 0, 0, 0, 0] 10 0 Quadratic 0 4.78 0.00 [1, 1, ?1, ?1, ?1, ?1, 1, 1] 11 0 Quadratic 1 0.87 0.00 [?1, ?1, ?1, ?1, 1, 1, 1, 1] 12 0 Quadratic 2 0.000 0.00 [?1, ?1, ?1, ?1, 1, 1, 1, 1] 13 0.05 2 ?0.003 0.76 deg/?hr 14 0.1 2 ?0.065 1.55 deg/?hr 15 0.1 Linear 1 ?0.006 1.54 deg/?hr [1, 0, 0, 0, 0, 0, 0, 0]

[0092] It is noted that Table 2 illustrates results of computer simulations performed with the invented fitting method using a testbed with eight gyroscopes.

[0093] Angular random walk (ARW) was set to zero in the first 12 computer simulations.

[0094] Computer simulations #1 to #5 represent examples in which time dependency of rate drift of any individual gyroscope is linear. Nature of time dependency is determined by multipliers 0, 1 and ?1, wherein 0 represents no rate drift, and 1 and ?1 refer to rate drift in opposite directions.

[0095] Computer simulation #1 is a reference case in which all eight gyroscopes have similar linear rate drift. A zero-order fitting function, in other words a constant value is applied as fitting function and averaging samples over the test period removes all errors, mean error in the angle is zero degrees and standard deviation (std) of the error is also zero. Thus, north angle determined based on the measurements is accurate. Although not shown in the table, value of Earth's rotation rate calculated is also accurate.

[0096] Computer simulation #2 illustrates a case in which one of the eight gyroscopes has a rate drift pattern that deviates from rate drift of other seven gyroscopes. If a zero order, constant function is used for fitting, mean value of the systematic error in the output is 1.38 degrees. Computer simulations #3 and #4 add further gyroscopes with non-zero linear rate drift characteristics, which further increases the systematic error in the mean value of the output, if a zero-order fitting function is applied. This is apparently not doing a good job for fitting.

[0097] Computer simulation #5 applies linear fitting according to the invention in the exemplary case otherwise equal to computer simulation #4. This fitting effectively removes the systematic error in its entirety, and north angle is accurate.

[0098] Computer simulations #6 to #12 represent examples in which time dependency of rate drift of an individual gyroscope is quadratic. In other words, the rate drift can be represented by a second-order function.

[0099] Computer simulation #6 is a reference in which rate drifts of all gyroscopes are similar. Like in the computer simulation #1, a constant fitting function can be applied and averaging samples over the test period effectively removes all systematic error, mean error in the angle is zero degrees and standard deviation (std) of the systematic error is also zero.

[0100] Computer simulations #7, #8 and #9 represent a case in which one of the gyroscopes has a quadratic rate drift while other seven gyroscopes are driftless. In the computer simulation #7, constant value fitting is applied. In the computer simulation #8, linear fitting is applied and in the computer simulation #8 a quadratic fitting is applied. Even the linear fitting (#7) gives fairly good systematic error correction, while using quadratic fitting (#8) in time domain enables fully removing systematic error in the output.

[0101] Computer simulations #10, #11 and #12 represent a case in which all gyroscopes have quadratic rate drifts, but in different directions. Constant fitting (#10) leaves a significant systematic error in the output. Use of linear fitting (#11) reduces output systematic error to less than one degree, while use of quadratic fitting (#12) effectively removes all systematic error in the output and also has zero standard deviation.

[0102] Computer simulations #13, #14 and #15 represent a situation with angular random walk (ARW) is present. Computer simulations #13 and #14 indicate that with second order rate fitting, also effects of white noise can effectively be mitigated, although some standard deviation remains in the result due to sporadic characteristic of the ARW. Computer simulation #15 indicates that even if there is a gyroscope that has linear rate drift in addition to presence of ARW, linear fitting can effectively remove all systematic error from the output.

[0103] It is noted that the exemplary fitting algorithm as described herein not only enables removing systematic error. In particular, the algorithm also enables shortening integration times required for removing effects of angular random walk and using longer integration times, since effects of rate drifts do not dominate performance even if integration time is long. Thus, waiting time for receiving relatively accurate north angle information during warm-up time of the device can be shortened.

[0104] The algorithm was found to be effective even with just two gyroscopes, especially when there is a 180-degree offset angle between their sense axes. However, when designing an actual functional gyrocompass, other error causing phenomena should be taken into account. For improving effect of angular random walk (ARW), it is preferred to have more than two gyroscopes, preferably arranged such that vector sum of sense axes is zero. An embodiment with three gyroscopes arranged at 120-degree angles between their sense axes has a vector sum equal to zero always has at least two gyroscopes providing a non-zero signal. Another embodiment is to arrange gyroscopes pairwise such that sense axes in each pair are offset by 180 degrees, which provides strongest possible output signal and thus has the best signal-to-noise ratio.

[0105] FIG. 13 illustrates schematically main functional elements of an exemplary MEMS gyrocompass according to an exemplary aspect. The MEMS gyrocompass comprises MEMS gyroscopes (10), at least one processing device (20) arranged in communication with the MEMS gyroscopes (10) and configured to obtain samples of outputs of the MEMS gyroscopes (10). Number of MEMS gyroscopes (10) is a design option according to exemplary aspects, but specific arrangements of the MEMS gyroscopes (10) have been recognized to have different technical benefits. The processing device (20) is communicatively coupled to a memory (21) comprising code which, when executed by the processing device (20) processes samples obtained from the outputs of the MEMS gyroscopes (10) to perform the exemplary algorithms described herein and thus provide an indication of north angle at an output (22).

[0106] In general, it is noted that the exemplary embodiments described above are intended to facilitate the understanding of the present invention and are not intended to limit the interpretation of the present invention. The present invention may be modified and/or improved without departing from the spirit and scope thereof, and equivalents thereof are also included in the present invention. That is, exemplary embodiments obtained by those skilled in the art applying design change as appropriate on the embodiments are also included in the scope of the present invention as long as the obtained embodiments have the features of the present invention. For example, each of the elements included in each of the embodiments, and arrangement, materials, conditions, shapes, sizes, and the like thereof are not limited to those exemplified above and may be modified as appropriate. It is to be understood that the exemplary embodiments are merely illustrative, partial substitutions or combinations of the configurations described in the different embodiments are possible to be made, and configurations obtained by such substitutions or combinations are also included in the scope of the present invention as long as they have the features of the present invention.