Angular Misalignment Calibration for Dual-Antenna GNSS/IMU Navigation Sensor

We address the angular misalignment calibration problem, which arises when a multi-antenna GNSS serves as a source of aiding information for inertial sensors in an integrated navigation system. Antennas usually occupy some outside structure of the moving carrier object, whilst an inertial measurement unit typically remains inside. Especially when using low- or mid-grade MEMS gyroscopes and accelerometers, it is either impossible or impractical to physically align IMU-sensitive axes and GNSS antenna baselines within some 1–3 degrees due to the micromechanical nature of the inertial sensors: they are just too small to have any physical reference features to align to. However, in some applications, it is desirable to line up all sensors within a fraction-of-a-degree level of accuracy. One may imagine solving this problem via the long-term averaging of sensor signals in different positions to ensure observability and then using angle differences for analytical compensation. We suggest faster calibration in special rotations using sensor fusion. Apart from quicker convergence, this method also accounts for run-to-run inertial sensor bias instability. In addition, it allows further on-the-fly finer calibration in the background when the navigation system performs its regular operation, and carrier objects may undergo gradual deformations of its structure over the years.


Introduction
One of the key problems in using low-grade inertial measurement units (IMU) is their inability to perceive azimuth without external aids. High-grade strapdown inertial navigation systems (INS) do this by measuring the Earth's rotation rate components in their instrumental axes. Apart from direct vector matching, there exist a number of approaches [1], but the crux remains in the accuracy of angular rate sensors (gyroscopes) being well within a small fraction of the Earth's angular rate magnitude. For microelectromechanical sensors (MEMS), run-to-run bias instability typically exceeds this requirement by 1-3 orders of magnitude, making conventional azimuth perception virtually impossible. For sensors having better in-run stability, special gyrocompassing methods have been developed [2], which require some special rotation of the IMU and a decent time span. In some cases, these methods highly rely on the physical stability of the base, so that the slightest mechanical perturbation may entirely ruin the solution. Even conventional sensor fusion with a single-antenna global navigation satellite system (GNSS) has its capabilities quite limited in estimating azimuth attitude error for low-grade gyroscopes.
To tackle this issue, a dual-antenna GNSS setup has become rather popular in a variety of applications [3][4][5][6][7][8] since cheaper and more compact GNSS hardware had emerged on the market. Its baseline vector, i.e., the vector connecting phase centers of the antennas (see Figure 1), being both known in the body reference frame of a vehicle and measured in a navigation reference (say, East-North-Up axes) using GNSS carrier phase observables and an RTK (real-time kinematics) approach [9], yields two attitude angles almost instantly as compared to MEMS gyrocompassing with no special maneuvering necessary. If the baseline lies sideways relative to the carrier, it provides yaw and roll angles of the body, with pitch, i.e., the rotation angle around the baseline itself, remaining unknown. Still, the inertial sensor provides pitch orientation, so that a full attitude solution becomes available in this integrated sensor system. However, some use cases require a sub-degree level of attitude accuracy, which implies both dual-antenna GNSS and inertial systems well aligned within the vehicle's frame of reference. For GNSS antennas with proper phase center calibration, their locations usually match the intended positions in technical drawings. The latter generally have millimeterlevel precision, so that corresponding angular deviations on a 1 m baseline become small fractions of a degree. They do less well with the inertial sensors. Due to the size of MEMS gyroscopes and accelerometers being extremely small, they barely have physical features to align them better than within some 3°. Therefore, after installation, the IMU module has some unknown angular misalignment with respect to GNSS antennas.
Due to its rather niche application, only a few works address the above issue [4,6,10]. We have been unable to find published methods, which still may exist in the Web. From personal communication, it appears that angular misalignment calibration either requires a special static experiment over a long period of time when both attitude solutions are averaged, or it emerges from the manual analysis of deviations in the integrated attitude solution. While the first solution relies upon inertial sensor biases being stable enough and takes a lot of time, the drawbacks of the second approach are self-evident.
In this paper, we suggest a solution to the misalignment problem via its calibration based on sensor fusion algorithms in a special experiment. Apart from faster convergence, this method accounts for run-to-run inertial sensor bias instability. In addition, it allows further on-the fly finer calibration in the background when a navigation system performs its regular operation, and the carrier object may undergo gradual deformations of its structure over the years.
While our calibration method is based on conventional Kalman filtering and INS error equations, we have identified four key issues that appear to be essential to solving the problem. They have not been self-apparent prior to approaching the problem, and therefore, we would like to emphasize these issues for those who may be working on similar problems:

•
Ensuring observability and good estimability in the underlying estimation problem requires specific rather intensive rotations; • For those above-mentioned rotations, attitude integration algorithms must provide a numerical solution which conforms to the INS error equations closely enough; for example, simplest single-step Euler integration does not do so; • Timing delays between GNSS and IMU measurements should be a part of the estimation problem, since their magnitude even of a fraction of an IMU time step is significant enough; • In addition, gyroscope measurements should be extrapolated accordingly when transforming GNSS antenna velocity into the velocity of the IMU.
In the following sections, we formulate misalignment calibration as an optimal estimation problem for a dynamic system with measurements. Its mechanization equations are based on INS error equations, with additional parameters being modeled as constants. For the complete rationale and derivation, please see the sections below. We consider the following sensor setup shown in Figure 2 with the notation listed in Table 1. Figure 2. Sensor setup. Two antennas (A 1 and A 2 ) of a dual-antenna GNSS unit reside at locations given by lever arms¯ 1 and¯ 2 relative to the inertial measurement unit with its reference point M. IMU instrumental axes z 1 , z 2 and z 3 are fixed to inertial sensor array, while lever arm vectors are known in a carrier body reference frame. Three components κ 1 , κ 2 and κ 3 of the Euler rotation vector define the slight misalignment between the two reference frames. Table 1. General notation.

Symbol Meaning
A 1 , A 2 , M GNSS antenna locations and IMU reference point, respectivelȳ 1 ,¯ 2 antenna lever arm vectors, i.e.,¯ k = MA k b 1 , b 2 , b 3 carrier body reference frame: b 1 -roll axis, b 2 -normal axis, b 3 -transverse axis z 1 , z 2 , z 3 respective IMU instrumental axes κ 1 , κ 2 , κ 3 components of the Euler rotation vector transforming from Mb frame to Mz x 1 , x 2 , x 3 navigation frame: IMU attitude matrix, transforming from Mx frame to Mz γ, θ, ψ corresponding roll, pitch and true heading angles, respectively ω z , f z angular rate and specific force vectors as projected onto the instrumental axes ω z , f z their components as measured by respective inertial sensors g,ū the Earth's local gravity and rotation velocity vectors ϕ, λ, h geodetic coordinates (latitude, longitude, altitude) V IMU velocity vector t time We assume every quantity as a function of time, so that t may appear as its argument, if necessary. However, in most cases, it is omitted for brevity.
Each vector symbol may have a subscript letter denoting either a reference frame in which its components are considered (e.g., x, b, z), or a specific axis (E, N, U, etc.) The dot above any quantity[ ] represents its derivative over time. For the vector product with any vectorv, we introduce a linear operator v so that in its coordinate form, the vector product becomes multiplication by a matrix: Let I be the identity matrix of the appropriate size whenever it is being used. We consider all reference frames to be right-handed orthogonal, so that the IMU has been calibrated with sufficient accuracy. The relation between attitude angles and attitude matrix is as follows: cosθ cosψ sinθ −sinθ sinψ cosγ+cosψ sinγ −sinθ cosψ cosγ−sinψ sinγ cosθ cosγ sinθ sinψ sinγ+cosψ cosγ sinθ cosψ sinγ−sinψ cosγ −cosθ sinγ   , (1)

Inertial Solution
To formulate the INS error equations, one must obtain the inertial solution first. It uses measurements from inertial sensors-angular rate sensors (gyroscopes) and accelerometersfor integrating equations of motion and yields attitude, position and velocity solution over time. The equations may appear in different form, and here, we use the attitude matrix and geodetic coordinates to use the same equations not only for calibration but in INS regular operation as well. They are as follows: where g x includes the centrifugal specific force component, and Ω is an angular velocity of the navigation reference frame Mx relative to the Earth, with its components being: with a and e being the Earth's ellipsoid semimajor axis and eccentricity, respectively. The first three equations in (2) actually describe a radius vector r x in some Earth-centered Cartesian reference frame. These two variants may be used interchangeably.
It is common to omit some of the terms such as Coriolis acceleration or even the Earth's rotation in the above equations for lower grade inertial units. In the following derivations we will, however, keep these terms for the pure sake of mathematical rigor. In real-world applications, the corresponding algorithms may happen to be previously implemented and tested in navigation software or libraries. Given that for modern processors, the additional computational burden often appears neglectable, we find it appropriate to leave for each reader the decision of whether to simplify the equations or not.
In addition, in real navigation systems, the equations for the altitude and vertical velocity component in (2) introduce well-known exponential instability [1]. So, instead of integrating them, the system uses an external source for altitude. In our misalignment calibration, however, we may benefit from using vertical velocity measurements from GNSS, like from using horizontal ones. We therefore keep these equations from being integrated and use them to form INS error equations as well.
The starting position and velocity in (2) are trivial to specify, with the coordinates of the calibration experiment known and velocity being zero. Attitude determination, however, requires a special procedure called initial alignment.

Initial Alignment Procedure
In this section, we describe the procedure required to obtain initial estimates for the attitude matrix L 0 . For doing that, the IMU remains at rest on the ground for some time t 0 with its normal axis pointing approximately upwards. While being stationary, its accelerometers measure the ground reaction force opposite to gravity acceleration, so that From the above, after averaging accelerometer outputs over the time period t ∈ [0, t 0 ], one may estimate IMU roll and pitch angles, respectively, as with the prime [ ] symbol meaning a value derived from measurements, and angle brackets · for averaging over t ∈ [0, t 0 ]. The two-antenna GNSS solution then provides an estimate ψ (t 0 ) for the azimuth angle up to some misalignment and other errors. According to (1), the estimated initial attitude matrix becomes: In addition to obtaining the initial attitude matrix, it is usually makes sense for MEMS gyroscopes to obtain rough estimates for their in-run biasesν 0 z1 ,ν 0 z2 ,ν 0 z3 , since for most devices, they exceed tens of degrees per second, being greater than or comparable to the Earth's angular rate of 15°/hr. Given the IMU is stationary during the initial alignment, the simplest form for those estimates is: According to Section 2.1, there is one more feature to the initial alignment in our misalignment calibration experiment. Namely, we are going to integrate the equations of motion (2) along the vertical axis. To avoid introducing exponential instability into the solution, we use a constant gravity model for our misalignment experiment with a gravity acceleration value of The model (8) certainly has some constant bias ∆g U produced by accelerometer errors. However, this bias appears to introduce no error into calibration, being properly estimated along with other parameters (see Section 2.4 for details).

Attitude Integration
For low-grade strapdown inertial systems such as MEMS-based IMUs, one usually implements simpler versions of attitude integration algorithms such as the Euler method for quaternions [11]. However, our simulation has shown (see Section 3.4) that motion patterns which provide better estimability properties of the misalignment calibration should include some kind of conical rotation. Under such motion, those simpler methods tend to introduce significant numerical errors, which do not obey INS error equations. Being systematic, they in turn produce biased estimates in calibration. With that knowledge, we have decided to use a more accurate version of the attitude integration algorithm based on the Bortz kinematic equation [12] for a Euler rotation vectorφ:φ For (9), we use an approximation of the 4-th order Runge-Kutta integration method. Although it may seem excessive to use it for a low-grade IMU, one should keep in mind that its errors are either systematic and closely conformant to INS error equations or stochastic with a nearly zero mean cumulative effect. Moreover, since we have used this algorithm for processing simulated data, it seems consistent to use it for real experiments as well. For the transition between two time instants t i−1 and t i with time step ∆t between them, we have: whereω(t i−1 ) andω(t i ) approximate the instant rotation rate vector using gyroscope measurements, which are the average angular rate components over the respective time step. The function F(ω z , φ z ) is the fourth-order Taylor expansion of the right-hand part of the Bortz equation with φ z 1. The calculated Euler vector increment ∆φ z yields a transition matrix C via Euler-Rodrigues' rotation formula [12] as follows: Together with the transition matrix for the navigation frame using the regular Euler method, we perform mechanization for the attitude matrix L from a time instant t i−1 to t i : x , Using (12), we obtain a calculated attitude matrix L over time, starting with L (t 0 ) from the initial alignment procedure.

Position and Velocity Integration
For the position and velocity, the conventional modified Euler integration has proven to work well, so that according to (2): with an appropriate gravity model for g x , andL(t i− 1 /2 ) being an estimate for the mid-step attitude matrix. Our calibration experiment does not include active linear motion, so (13) may be simplified. The reason for not ignoring here the Coriolis term and the rotation of navigation frame (Ω x ) is our future plan to use the same equations and models for in-run system calibration in its regular operation.
As for the gravity model g x , we use a constant value obtained in the initial alignment, as per Section 2.1.1, for our nearly static calibration experiment to avoid exponential instability. In other cases, appropriate gravity models may be used for integration, which are provided with an external altitude information.

INS Instrumental Errors Model
In general, the choice of mathematical model of INS instrumental errors heavily relies on accuracy class of the INS. For us, the subject is a low-or mid-grade MEMS IMU. We assume that before calibrating the angular misalignment, the inertial sensors themselves are pre-calibrated, so that standard parameters of an INS instrumental errors model, i.e., constant biases, scaling coefficients, etc., are compensated using one of the known methods [13][14][15]. In addition, temperature variations of inertial sensor measurements are not considered in this research. We assume that IMU thermal calibration can be carried out in advance [16,17], and residual errors are stochastic. Otherwise, we may suggest performing the misalignment calibration at a constant temperature. Let us define accelerometers and gyroscopes instrumental errors as where the minus sign in the second expression of (14) originates from a tradition for gimbaled INS. We accept the following model for instrumental errors of accelerometers and gyroscopes where ∆ f 0 z , ν 0 z are null biases of accelerometers and gyroscopes, respectively, ∆ f s z , ν s z are stochastic terms of the measurement error with known a priori moments. We include null biases into the estimation process since they generally happen to be different and not very stable in each INS run as opposed to scaling factors and other parameters. For their instability, we accept Wiener processes: where q ∆ f 0 , q ν 0 represent independent white noise processes with known a priori moments.

INS Error Equations
The INS errors in the geodetic navigation frame are as follows: We further consider the behaviour of INS errors over time up to linear terms. We use INS error equations [18] in the computed geodetic navigation frame y (My 1 y 2 y 3 ). The INS error equations will serve as a dynamic model in the linear estimation problem. We adapt these equations for the case of specific calibration motions, the choice of which will be explained in Section 3.4. The modification consists of replacing the term 2ω 2 0 ∆r U with ω 0 being the Schuler frequency, by the constant error of local gravity force ∆g U resulting from accelerometer errors in initial alignment as described in Section 2.1.1. Thus, we have: . . .
To make INS equations less cumbersome, we omit primes in all coefficients of the INS errors because the equations still hold true to within linear approximation. In addition, we consider the continuous-time version of a linear dynamic system for the sake of notation's simplicity.

Measurements and the Estimation Problem
We use a commonly accepted loosely coupled GNSS-INS integration scheme with the feedback into inertial solution, the reason behind being its equivalence with the tightly coupled integration under a sufficient number of GNSS measurements. We assume that the misalignment calibration experiment is carried out in a favorable GNSS environment. In the estimation, we use a GNSS-derived position and velocity solution obtained from Doppler measurements. For them, timing errors play a significant role, so they are described separately below.

Time Synchronization Errors between INS and GNSS
The benefit and methods of estimating the time synchronization errors between INS and GNSS are shown both by numerical simulation of low-cost GNSS-aided INS integration with feedback [19] and by the processing of real data from aircraft flights with a strapdown INS [20]. For our research, we have also performed the numerical simulation, which supports the above results (see Section 3.4).
Let r x (t) and V x (t) be the calculated INS position and velocity, respectively, computed at some time t. We then define the time synchronization errors τ pos, k , τ vel, k for two antennas (k = 1, 2) between INS and GNSS solutions as follows: V k x (t − τ vel, k ) is a GNSS-derived velocity of the k-th GNSS antenna computed at time t.
We assume magnitudes of time synchronization errors to lie typically within 0.1-0.2 s, i.e., within a few GNSS time steps. Hence, we accept the following relations: In (18) and (19) and further on, we specify only time instants different from t.

Angular Misalignment between Instrumental and Body Frames
In this section, we mathematically formulate the problem of the angular misalignment between INS and dual-antenna GNSS. Let us recall the underlying assumptions: 1.
Origins of the instrumental and body frames are the same; 2.
Constant lever arms of two GNSS antennas A 1 , A 2 in the body frame 1 b , 2 b are known; 3.
Instrumental and body frames slightly differ.
The first proposition follows from the definition of the body frame (see Table 1). As for the second assumption, the components 1 b , 2 b are found or calculated from the technical documentation of the carrier object. If the angular misalignment between reference frames is large, one can deduce its approximate magnitude from the same technical documentation, thus reducing the problem to small angles, as the third assumption states.
From the assumptions above, for¯ 1 ,¯ 2 we have i.e., the attitude of the IMU instrumental frame z relative to the body frame b is determined by Euler rotation vector κ b = (κ 1 , κ 2 , κ 3 ) T . The constancy of κ b fully depends on the carrier object being rigid and stiff enough for GNSS antennas and IMU spatial separation to stay the same, so we believe thatκ

INS Attitude Errors
The INS error equations given in Section 2.3 imply several reference frames as per Table 2 below. Axes Description computed navigation frame as a result of INS coordinate errors z 1 , z 2 , z 3 computed instrumental frame as result of applying operator L to x axes y 1 , y 2 , y 3 computed navigation frame as result of applying operator L T to z axes The origin of true frames x, z is the IMU reference point M. The origin of computed frames x , y, z is a computed IMU position M . Euler rotation vectors α 0 x , β 0 x and γ 0 x represent transformations between close reference frames x, x , and y, so that for the components of any vector v, we have the following: with the linear relation between Euler rotation vectors α 0 x , β 0 x , γ 0 x and INS errors represented as follows: Thus, introducing INS attitude errors, we follow the ideology of [18] with slightly different notation. The complexity may seem excessive, but it keeps mathematical rigor in our derivations.

Position Measurements
Letr k be the radius vector for the k-th GNSS antenna (k = 1, 2), as derived from GNSS pseudoranges [9] and converted to the Earth-centered geodetic navigation frame x. The residual position measurement for our linear estimation problem then becomes: Adding and subtracting r x from the right part of (26), substituting (18), (21), (24), (25) into (26), and using the relation r x + k x = r k x , yield a linearized model for the residual position measurement at the GNSS epoch t for the k-th antenna: where ρ pos, k is a stochastic error with a priori known moments.

Velocity Measurements
Let Θ z be the angular velocity of the frame Mz relative to the Earth so that Θ z = ω z − Lu x . Having the GNSS velocity solution derived from Doppler measurements [21], we introduce a residual velocity measurement: Similarly to the derivation of (27), we add and subtract V x from the right part of (28), substitute (19), (21), (24) into (28), and use the relation for linear velocities of two points (A k and M) of a rigid body V k x = V x − L T Θ z z . The linearized model for the residual velocity measurement at GNSS epoch t for the k-th antenna then becomes: where ρ vel, k contains both GNSS measurement noise and the gyroscope stochastic term ν s z , whose moments are known. The coefficient of τ vel, k in (29), taking (2) into account, can be expressed up to linear terms aṡ Note that (30) may be further simplified if the IMU is stationary, so that V x ≈ Ω x ≈ 0.

Estimation Problem
Thus, we have reduced the problem of calibration of angular misalignment between dual-antenna GNSS and IMU to a linear stochastic estimation problem with the following 23-dimensional state space vector: The dynamic model for (31) consists of (16), (17), (20), and (22). Equations (27) and (29) form the measurement model. For processing, we use the discrete-time equivalent of these equations. The initial estimate of the state space vector is zero, and the initial covariance matrix of the estimation error is known a priori.
There are many methods based on the Kalman filtering approach [22], which provide estimates for the state vector components (31) either in real time or in post-processing. We use the Potter square root filter version based on Cholesky covariance factorization [23].

Motivation
After the calibration has been formulated as an estimation problem in the above Section 2.4, we need to ensure its good estimability properties in terms of converging error covariance [24]. As it is usually the case in INS sensor fusion algorithms, its properties mostly depend on the motion of the IMU. Time-varying coefficients ω z (t), f z (t), and L(t) in (27) and (29) define how well components of the state vector will separate from each other in the estimation process, or, in other words, how differently these components will manifest themselves in INS errors. For higher-order time-varying systems such as that being under consideration, predicting its properties analytically from the equations alone is hardly a solvable task in general. In fact, only numerical analysis may provide practical insights for most systems of this kind. Therefore, testing a range of scenarios for the calibration experiment, and choosing the right one to be actually executed, become primary reasons to perform numerical simulation. The secondary reason being mere verification of the consistency between our models and algorithms is also important to support future conclusions.
In our case, it was the simulation that has forced us to take into account effects which a priori seemed quite neglectable even to experts in the field. All in all, we have decided to describe the numerical simulation as an inseparable part of our research.

Inertial Sensors
As a matter of fact, rotational motion is crucial to calibrating angular misalignment. To simulate such motion, we have developed a virtual three-axis turntable. Each axis, being controlled and simulated individually, can perform a number of commands. Their list includes moving into an arbitrary pre-selected position, uniform rotation at a given rate, harmonic oscillations and stopping the rotation.
For each axis of our virtual turntable, according to a pre-defined set of specific commands, we create a twice continuously differentiable analytical function, representing the angle of rotation around this particular axis over time. Combining rotations for three axes then allows us to simulate a wide range of complex motion patterns and calibration scenarios. Using analytic functions for rotation angles, we derive the absolute angular rate ω z (t) and the specific force f z (t) vectors as projected onto the IMU instrumental reference frame. We then calculate their components on a discrete time grid at a high rate. Downsampling to a desired IMU sampling rate of, say, 250 Hz, using arithmetic average, completes the simulation process.
The higher internal frequency of the simulation allows us to properly reproduce the integration (or averaging) which occurs in real inertial sensors. Internal frequency may be set as high as it is required for a given rotation pattern. In practice, one should try larger and larger values until the change in navigation solution becomes negligible. In our case, for a 250 Hz IMU sampling rate, a 256 times higher simulation frequency of 64 kHz has happened to be enough.
For now, the IMU reference point M remains stationary in our simulation. Inertial sensor errors satisfy the model (15), being added when appropriate.

GNSS Measurements
Having the IMU rotation ready, so that we may assume angular rate vector components ω z (t i ) and attitude matrix L(t i ) to be known at a discrete time grid, the GNSS position and velocity of two antennas A 1 and A 2 need to be simulated (see Figure 2). For them, each lever arm vector 1 b and 2 b is known component-wise in some carrier body reference frame b. Angular misalignment angles κ 1 , κ 2 and κ 3 define an Euler rotation vector κ with the corresponding rotation matrix D according to (11): Using the above matrix, the antenna coordinates in geodetic Cartesian axes become with, longitude, latitude and altitude as well as radius vector appropriately calculated. Velocity vector derivation uses Euler's rotation formula, so that given the IMU does not perform linear motion, for each antenna, we have where the instant rotation rate components ω z are either taken directly from the simulation of inertial sensors before averaging or derived from the gyroscope output similar to (10). When necessary, we add GNSS solution errors to the antenna position and velocity as well. The GNSS position stochastic errors, although having a rather complicated nature in practice, happen to have a minor effect on the estimation of angular misalignment. For GNSS velocity derived from Doppler measurements, their stochastic errors appear to be quite close to white noise.

Simulation Results
To demonstrate the rationale behind certain decisions accepted in our calibration method, we have simulated the following effects listed in Table 3 below. Table 3. Simulated effects.

Description Parameters Typical Values
Different spatial configurations of sensor setup k b ∼ 1 m, Angular misalignment between b and z frames κ b ∼ 3°I nitial alignment errors ∆ψ(t 0 ) = β U ∼1°I nertial sensor errors ∆ f z , ν z ∼1 cm/s 2 , 10°/h GNSS velocity solution bias 1 cm/s Time delays of GNSS solution τ pos, k , τ vel, k ∼10 ms Before proceeding to choose the class of rotational motion for the calibration, one should note that in measurement model (27), (29) the coefficients of the desired κ b parameters contain constant factors k b , and for them, we have rank 1 b 2 b ≤ 3, so it can be potentially less than the dimension of κ b . Therefore, the choice of two GNSS antenna locations ( 1 b , 2 b ) with respect to the IMU has a direct impact on the estimability of all three parameters κ 1 , κ 2 , κ 3 . The necessary conditions for them to be estimable is 1 b = c 2 b for any non-zero constant c. In the simulation, this condition is satisfied at all times, unlike the next argument in Section 4 dedicated to real data processing, where the actual sensor setup did not allow for that.
For the misalignment calibration problem, we suggest using so-called conical motions. They comprise simultaneous rotations around two perpendicular axes. For them to be carried out on a turntable, two of its axes perform harmonic oscillations out of phase by a quarter of a full period with each other. The term conical arises from one of the instrumental axes moving along the generatrix of a certain circular or elliptical cone (depending on amplitudes of the above harmonic oscillations). Figure 3 demonstrates the difference in estimation process for two different rotation types in terms of κ b . The first type of motion (see the left inset) is a sequence of four rotations by 90°around vertical axis x 3 , which approximately coincides with the second instrumental one z 2 . After each rotation by 90°, there is a static position. Axes z 1 and z 3 lie near the horizontal plane x 1 x 2 . As a result, two GNSS antennas move along horizontal circles. Rotations such as these, so that the IMU has different heading angles with a roughly 90°increment, are similar to maytagging-a conventional technique used for gyrocompassing using low-grade inertial sensors [2]. The second type of rotation (see the right inset of Figure 3) is the conical motion described above. The figure clearly shows that conical rotation provides a strictly monotonic decrease of the estimate error covariance converging to zero over time, while the maytagging does not do so.
The second point which the numerical simulation has given us an insight into refers to the timing errors between inertial and GNSS information while performing the angular misalignment calibration. Although it is almost self-evident that these should be taken into consideration, the actual figures have become a surprise even for experienced engineers working in the field. Figure 4 indicates that even 12-millisecond timing errors τ vel, k , (k = 1, 2), which correspond to only 3 IMU samples at 250 Hz, can be critical. The dashed lines stand for the errors in estimating τ vel, k (left plot) and κ b (right plot) when calibration models include the timing skew. Solid lines represent the case when τ vel, k is omitted under simulated 12-millisecond delay in GNSS measurements. It apparently becomes an issue for estimating κ 2 , κ 3 , with their estimates swaying away from reference values. In addition to Figure 3, the plots below confirm that once the conical rotation starts, the estimated misalignment errors immediately begin to converge. Timing errors τ vel, k appear to have good estimability right away from the very first rotation.  From the same assessment, it follows that even phase delays of a fraction of the inertial time step in angular rate measurements should be accounted for in both simulation (34) and estimation (30), since they produce significant calibration errors, albeit not as large as in the example that Figure 4 illustrates.
Another issue that the numerical simulation has revealed appeared to be a substantial difference between attitude integration methods. Although simpler integration methods are usually considered sufficient for MEMS gyroscopes, Figure 5 indicates that replacing the method alone changes the estimation error completely once conical rotation starts. For conventional Euler integration, which accounts only for the first term in the Bortz kinematic Equation (9), estimation errors of κ 2 and β U shown in green and yellow, respectively, do not converge to zero over time. We believe that this is due to the fact that a more accurate attitude integration algorithm provides errors much more closely conforming to the INS error Equation (15) for their systematic parts.
In this section, we deliberately do not show any results with sensor errors containing stochastic terms, because qualitative analysis does not depend on them. In fact, such simulations with Gaussian noises were carried out as well to test our processing software. Overall, our experience has shown that estimation problems of this kind are barely solvable in practice without proper simulation.

Experimental Results
For preliminary validation, two similar calibration experiments have been performed using the setup shown in Figure 6. In these experiments, we used high-precision GNSS equipment, namely Javad™ Prego ® receivers and AirAnt ® antennas. For inertial sensors, we took an iSense™ AIST-350 thermally stabilized MEMS IMU based on LPY510 gyroscopes and ADXL326 accelerometers by ST Microelectronics™ and Analog Devices™, respectively. GNSS receivers were operating at 10 Hz, while IMU records had a 250 Hz sampling rate. Both experiments, after the 4 min initial alignment phase, included three types of rotation for approximately 6 min each, with 4 min static intervals between them:

1.
Oscillations around the roll axis (referred to as γ(t) oscillations below) with amplitude up to 10°and nearly constant yaw ψ(t) and pitch θ(t) angles; 2.
Oscillations around the yaw axis (referred to as ψ(t) oscillations below) with amplitude around 15°with nearly constant roll γ(t) and pitch θ(t) angles, both close to zero; 3.
Conical rotations, in which IMU performed minor linear motion, and GNSS antennas were traveling along circles in anti-phase with each other, antenna lever arm vectors sweeping along generatrix lines with 40-50°aperture, at a period of around 2-4 s. These motions may be considered as a composition of the previous two.
All rotations were performed manually. Lever arms for both antennas had their lengths around 0.8 m. Note that in our experiment, lever arm vectors¯ 1 and¯ 2 happen to be collinear, so that the IMU reference point M lies on their baseline A 1 A 2 . Hence, parameter κ 3 is not estimable under the given geometry. From this point on, we consider the estimability of only two parameters κ 1 and κ 2 for the angular misalignment between IMU and dual-antenna GNSS reference frames. Their true values, as described in the Introduction, remain unknown, and there exist no reasonable means of measuring them directly. However, since between the two experiments, our instrumental setup has not changed, we expect estimates for κ 1 and κ 2 to repeat. Figure 7 illustrates the estimation process over time for both experiments throughout different types of rotation. Estimates for κ 1 and κ 2 appear as solid lines along with their predicted 2-σ confidence intervals. The latter seem to well overlap by the end of the calibration, indicating that the estimates are consistent with each other in two experiments at a sub-degree precision. Experimental data used in this Section are available in Supplementary materials for processing. Performing the first two types of rotation in our experiment was motivated by the potential possibility of replacing conical motion, being more difficult to implement in practice, with two of its constituents, namely rotations around each axis individually. To address this, let us refer to Table 4 with Figure 7 serving for further clarification. he estimated error covariance for κ 1 , κ 2 noticeably decreasing after the execution of conical motion implies that rotations around roll and yaw axes separately do not provide proper convergence. So, from the covariance analysis, we recognize that conical rotation is the preferred type of motion ensuring estimability for the angular misalignment calibration in real experiments and further supporting the results of numerical simulation.

Discussion
We have reduced the problem of angular misalignment calibration between the instrumental reference frame associated with an IMU, and the carrier body reference frame with known locations of two GNSS antennas in it, to a conventional linear stochastic estimation problem. The research is relevant to all applications aimed at tracking orientation using a low-grade IMU and dual-antenna GNSS within a sub-degree level of precision. Appropriate measurement models have been derived to extend conventional loosely coupled GNSS/INS sensor fusion filtering for including parameters required for the misalignment calibration. It uses a GNSS position solution and velocity derived from Doppler observables. Since the algorithm has inertial sensor biases in its state vector subject to estimation, it is expected to be immune to run-to-run bias change inherent to most lower-grade inertial sensors. Presumably, after compensating for the initial misalignment, the estimation may continue running in the background to account for slower structural deformations over time. This, however, requires additional research for confirmation.
The numerical simulation in Section 3.4 has shown that: 1.
The preferred motion pattern for calibrating angular misalignment includes conical rotation; 2.
The following key issues appeared to be essential to successful estimation: • Taking the time synchronization error between IMU and GNSS data into account at the few-millisecond-level; • The above includes phase delay inherent to integrating (or averaging) gyroscopes; • Modifying the attitude integration algorithm to produce errors properly obeying the INS error equations.
The processing of real experimental data has shown the feasibility of the proposed calibration method, and it produced consistent results in agreement with the numerical simulation. The final predicted standard deviation of the misalignment error does not exceed 0.25°.
However, more extensive validation is planned for the future, since the real misalignment angles, i.e., the "ground truth", seem to be practically unavailable in real applications that use MEMS sensors. The validation may be achieved through a variety of approaches, the first being increasing the number of experiments to a statistically significant amount. The second option is to perform the calibration using a high-grade inertial sensors. A navigation-grade INS is able to produce attitude angles autonomously, so that its instrumental frame may be directly aligned to GNSS antennas within some few arc minutes. After that, MEMS gyroscope and accelerometer data may be used to simulate low-grade IMU output for our calibration algorithm to be applied to. The third approach is an indirect one, the idea behind being to show that the navigation solution becomes more accurate after compensating for the estimated misalignment angles. Having its own importance in itself, this approach will become our primary focus for future research.
In this work, we used the GNSS-derived position and velocity obtained from code pseudoranges and Doppler measurements, respectively, since there exists GNSS equipment not able to record phase measurements. The first two measurement types seem to be available for an external processing in a wider range of GNSS devices rather than phase measurements. For many general applications, lower-accuracy GNSS sensors are preferable due to their smaller cost and size. Still, using phase measurements, and developing the corresponding extension to the estimation problem, including INS-aided carrier phase ambiguity resolution on-the-fly (i.e., deeply-coupled GNSS/INS for MEMS) is the next possible development branch in our research.