On the Recursive Joint Position and Attitude Determination in Multi-Antenna GNSS Platforms

: Global Navigation Satellite Systems’ (GNSS) carrier phase observations are fundamental in the provision of precise navigation for modern applications in intelligent transport systems. Differential precise positioning requires the use of a base station nearby the vehicle location, while attitude determination requires the vehicle to be equipped with a setup of multiple GNSS antennas. In the GNSS context, positioning and attitude determination have been traditionally tackled in a separate manner, thus losing valuable correlated information, and for the latter only in batch form. The main goal of this contribution is to shed some light on the recursive joint estimation of position and attitude in multi-antenna GNSS platforms. We propose a new formulation for the joint positioning and attitude (JPA) determination using quaternion rotations. A Bayesian recursive formulation for JPA is proposed, for which we derive a Kalman ﬁlter-like solution. To support the discussion and assess the performance of the new JPA, the proposed methodology is compared to standard approaches with actual data collected from a dynamic scenario under the inﬂuence of severe multipath effects.


Introduction
As contemporary applications such as driverless cars or autonomous shipping are called to revolutionize intelligent transportation systems (ITS), there is a growing need for the provision of precise and reliable navigation information. Global Navigation Satellite Systems (GNSS) play a fundamental role, becoming the main information supplier of positioning, navigation, and timing (PNT) data. While standard GNSS techniques-based on code observations-provide a decent performance for many applications, they do not comply with the far more stringent precision requirements of modern safety-critical scenarios. That is the reason why the transition to carrier phase-based techniques is required to reach precise navigation. Indeed, carrier phase observations present noise levels two orders of magnitude lower than their code counterpart. However, carrier phase observations are ambiguous, since only their fractional part is measured by the receiver [1]. The unknown number of integer cycles between satellite and receiver, so-called ambiguities, must be estimated to enable high-precision navigation. The ambiguities' estimation process is widely known as ambiguity resolution (AR) [2][3][4][5], which in turn results in a challenging estimation procedure.
(1) The GNSS carrier phase-based positioning and attitude models are revisited, and the connection to the JPA problem is formalized. (2) Application of Lie theory principles to the GNSS-based attitude problem (and consequently, also the JPA). Thus, the attitude is parametrized with the unit-quaternion, and a recursive solution based on an error state Kalman filter (ESKF) is formulated. (3) The performance of the proposed recursive solution for the JPA problem is addressed on a realistic signal-degraded scenario, and the comparison w.r.t separately solving the positioning and attitude problems (i.e., the standard solution) is analyzed.
The rest of the paper is organized as follows. Section 2 presents the notation employed throughout the paper, the coordinate frames involved in the estimation process, and the state definitions of the JPA. Then, Section 3 discusses the carrier phase-based positioning and attitude models. Section 4 relates the GNSS observations models to the JPA problem, and the Kalman Filter-like recursive solution is presented. The experimentation and posterior discussion are presented in Section 5. Finally, Section 7 concludes the paper with an outlook and discussion of future work directions.

Notation
Italics indicate a scalar quantity, as in c; lower case boldface indicates a column vector quantity, as in a; upper case boldface indicates a matrix quantity, as in A. The matrix/vector transpose is indicated by a superscript (·) as in A . [A, B] and [ A B ] denote the matrix resulting from the horizontal and the vertical concatenation of A and B, respectively. I M is the identity matrix of dimension M. 1 M,N is an M × N-dimensional matrix with all components equal to one. 0 M,N is an M × N-dimensional matrix with all components equal to zero. · denotes a Euclidean norm, and b 2 A = b A −1 b is the weighted inner norm. E(·) and D(·) represent the expectation and dispersion operator respectively. vec(·) is the operator that stacks the columns of a matrix A into a single vector. A = diag(·) indicates a diagonal matrix whose entries in the diagonal are given by (·). The Kronecker product is denoted with ⊗. The quaternion multiplication is denoted by •.

Coordinate Frames
In general, a navigation system involves the determination of the position, velocity, and attitude of a moving (rigid) body. The kinematic quantities generally involve two coordinate frames: (i) the frame whose motion is described, the body or vehicle frame B, and (ii) the frame with respect to which that motion is, denoted as the global or inertial frame G. For the problem at hand, the body frame is defined such that the X-, Y-, and Z-axes are aligned with the right, forward, and up directions of the vehicle, respectively, while the global frame corresponds to the Earth-centered Earth-fixed (ECEF) frame. A pictorial example of the aforementioned frames is given in Figure 1 (left). The left subscript on a vector indicates the frame in which that vector is expressed, while the right subscript refers to a specific position (i.e., any of the onboard vehicle antennas or the base station). In multi-antenna GNSS platforms, one generally considers the position of the master antenna as the center of the body frame and the unknown position to be determined G p G p m , while the positions of the remaining N slave antennas G p j (with j = 1, . . . , N) are referred to the former with the respective baseline vectors. For instance, the position in the body frame for the jth slave antenna is as: where the baseline vector B b j,m in the body frame is known a priori. Similarly, the position in the global frame for the jth antenna can be easily formulated upon the knowledge of the master position and the rotation relating the body and global frames.
where R ∈ SO(3) and q ∈ S 3 are the rotation matrix and unit-quaternion to represent the body-to-global rotation. The quaternion is written as q = [cos(φ/2), e sin(φ/2)] with φ a rotation angle and e a unit rotation axis. The quaternion rotation follows the Hamilton convention (real part first, right handedness). Henceforth, either R B u or q • B u • q * are used indistinctly for expressing the body-to-global rotation of a generic vector B u. However, since rotation matrices result in difficult renormalization and computational inefficiency, its use on recursive attitude determination is often discouraged [16][17][18]. Instead, this work focuses on the unit quaternion parametrization for attitude estimation. Additional details on the quaternion group and algebra are provided in Appendix A.

State Definitions
The proposed JPA problem relates to the estimation of the vehicle kinematics: position G p, velocity G v, and attitude q, as well as the GNSS integer ambiguities a for the n + 1 tracked satellites across the onboard and base station antennas, as later defined in Section 3. The state estimate x is expressed as a discrete-time state-space model, described at a given time t as: Optionally, the state vector could be augmented to account for the biases of an inertial sensor, whenever that navigation modality is fused along with GNSS. The state defined in (3) is often applied to describe the motion of vehicles following a constant-velocity non-turning movement [19,20], e.g., automobiles, vessels, or aeroplanes moving in a straight line. To describe a wider range of motions, (3) could be extended with the acceleration and angular rate of the vehicle, or an interacting multiple model (IMM) estimator could be applied [21,22].
Preserving the quaternion unit norm constraint is a challenging task, for which Lie theory provides useful algebraic tools whose application in navigation problems is increasingly used [23][24][25]. While the state x lives in a manifold, its perturbations δx belong to the tangent space of the manifold. The error state Kalman filter (ESKF) allows exploiting the Lie algebra by defining the "true" (unknown) state as the group composition of the total-statex and the error-state δx as x =x ⊕ δx [26][27][28]. At time t, the error state is described by: where the rotation vector δθ can be identified with its associated quaternion δq as follows: First, the rotation vector in the Euclidean space δθ ∈ R 3 is related to the Lie algebra eφ ∈ s 3 through the isomorphism (·) ∧ : R 3 → s 3 . Then, the Lie algebra is identified with the S 3 manifold via exponential mapping δq = exp(eφ) = [cos(φ/2), e sin(φ/2)] . Finally, the composition of the total quaternionq with the error quaternion δq is given by: The reader interested in the basics of Lie theory may refer to [29,30] or to a recent selection of these principles in relation to robotics in [31].

GNSS Carrier Phase-Based Positioning and Attitude Observation Models
Let us consider n + 1 GNSS satellites simultaneously tracked (single frequency) at one antenna (i.e., if considering N + 1 antennas at the vehicle, the observation model remains the same for each antenna) installed on a vehicle and at a base station of known coordinates. At a particular time, the code and phase observables for the ith satellite at the jth antenna (either onboard or base) are: where: ρ, Φ are the code and phase observations (m), Due to the influence of imprecise ephemeris and atmospheric-related errors, high precision cannot be achieved directly exploiting the observations on (7). Instead, the double-difference combination of observations is fundamental to enable precise positioning and attitude determination, as discussed next. For the remainder of the section, superscripts refer to a particular satellite (r for the pivot and i = 1, . . . , n for the others), while subscripts denote a particular receiving antenna (m for master antenna, j = 1, . . . , N for the slave antennas, and b for the base station), as illustrated in Figure 2.

RTK Positioning Model
RTK is a differential positioning procedure, for which the unknown position of a GNSS antenna is determined w.r.t. a stationary station of known coordinates [32]. To achieve centimeter-level precision positioning, the use of observation double-differencing and integer AR (IAR) is required (details in Section 3.2). Double-difference is the combination formed between two satellites and two receivers to eliminate or minimize nuisance parameters (i.e., atmospheric propagation delays, clock offsets, relativistic effects, etc.) from the involved observables. In the positioning model, the particular double-difference (DD) code and phase measurements (i.e., for the ith satellite, w.r.t. the base b and pivot satellite r, at the master antenna m) are given by: where the notation (·) i,r b,m refers to the DD observation conformed by the base station and master antennas, pivot, and ith satellite and (·) 1:n,r b,m is the corresponding n-dimensional vector of DD observations. The resulting positioning observation model is generally expressed in its linearized form as: where λ is the carrier wavelength, x pos = [ G b b,m , a pos ], a pos ∈ Z n is the vector of ambiguities, and G b b,m ∈ R 3 is the baseline vector between the master and base station positions. n pos is a zero-mean noise term with covariance Q y pos . The so-called design matrices are: where G is the geometry matrix composed by the satellite steering line-of-sight vectors w.r.t. the base . Re-arranging some terms in (10), the positioning model can be cast as: Resolving the baseline G b b,m between the base station and the master antenna grants knowing the position of the master antenna as: Notice that in contrast to (7), where the unknown ambiguities N i m ∈ R, because of double-differencing, we now have integer ambiguities a pos ∈ Z n . That is, to solve the estimation problem defined by the linearized Gaussian observation model (10), we have to resort to AR techniques (see Section 3.2), which applied to the previous equations, define the single-antenna RTK solution.

Integer Ambiguity Resolution for Positioning
The system of observation equations in (10) leads to an optimization problem with mixed real and integer parameter estimation. For simplicity, let us rename the variables for the positioning model in (10), i.e., y y pos , Q y Q y pos , a a pos , and b G b b,m . Thus, we want to solve: Due to the integer nature of a, a closed-form solution to (14) is not known. Instead, the work in [33] proposed the decomposition of the quadratic objective function (14) into the sum of three successive LS problems, where the first term on the right-hand side is the so-called float solution, min y − Aâ − Bb where Qb ,â refers to the cross-covariance matrix between the float estimates for position and ambiguities. A relevant remark is that, whenever the estimated integer ambiguities do not match the true ones, the fixed solution will be biased. The precision of the solution improves only when the correct ambiguities are estimated.

A Comment on Multi-Antenna DD Observations and Processing
Notice that the previous standard RTK solution was formulated for a single master antenna. One could think that having multiple (N + 1) antennas, i.e., N times more observations, would certainly improve the overall estimation. The trivial solution is to reformulate the model in (10) taking into account an N + 1 (i.e., master + N slave antennas) multi-antenna GNSS platform. Notice that in that case, we have N + 1 unknown baseline vectors between the base and the different antennas to be estimated alongside the n × (N + 1) integer ambiguities. The multi-antenna (MA) measurement vector is then: where y b,m corresponds to (9) and the DD observations formed by the base station and the slave antennas y b,j , j = 1, . . . , N, can be built as in (8) where the state to be estimated is now However, obviously, this is like considering N + 1 independent single-antenna RTK problems, because this trivial formulation does not take into account that all the antennas belong to the same rigid body, and therefore, there exists a link among antennas, which is disregarded. If we consider known baseline vectors among the different slave antennas and the master in the body frame, B b j,m (master to the jth antenna baseline), we can reformulate the previous model (18) only in terms of the master-to-base baseline, G b b,m , but this requires the platform orientation, i.e., its attitude. Therefore, to exploit multi-antenna RTK solutions correctly, we need to perform a JPA estimation, being the underlying motivation of this contribution.

GNSS Attitude Model
GNSS-based attitude determination can be realized for vehicles equipped with at least two GNSS antennas, whose positions in the local frame of the vehicle have been surveyed, i.e., the baseline among antennas in the body frame is known. Similarly to the RTK positioning introduced in Section 3.1, precise GNSS-based attitude estimation requires IAR and observation double-differencing. In this case, however, the DD combinations performed between the slaves and the master antennas (i.e., for the ith satellite, w.r.t. the jth antenna and pivot satellite r, at the master antenna m), The complete set of observations is gathered in vector y att ∈ R 2n×N , which stacks the observations formed with the combination of slave and master antennas y j,m as: ]. Then, the attitude model can be cast as a function of the quaternion-parametrized attitude operation: , D(y att ) = Q y att (21) where Z ∈ Z n,N is the matrix containing the DD ambiguities, which expressed in vector form is denoted as a att = vec(Z), a att ∈ Z n×N . Differently from the RTK positioning case, the GNSS attitude model presents additional challenges: (i) the double product of the quaternion for the rotation operation leads to a nonlinear equation (hence, the model being expressed as h(q) and not in matrix form); (ii) a nonlinear constraint to assure that the quaternion presents the unit norm (otherwise, it would not represent a proper rotation) is to be imposed. A short discussion on IAR techniques for the attitude model (21) is provided next, with the associated recursive solution being proposed in Section 4.

Integer Ambiguity Resolution for Attitude Determination
As previously discussed, the GNSS attitude model in (21) leads to a nonlinear optimization problem with mixed real and integer parameter estimation, subject to a nonlinear constraint. The optimization problem is formulated as: Notice that, although the constraint function is not made explicit (i.e., (22) is subject to q 2 = 1), such a constraint is made implicit by indicating that the quaternion belongs to the three-sphere q ∈ S 3 . Similarly to the RTK positioning problem, a closed-form solution to (22) is not known, and a decomposition into three sequential LS adjustments is applied: Similarly to the IAR for positioning in Section 3.2, the three consecutive estimates correspond to the float solution, the integer ambiguities, and the fixed solution, respectively, with slight differences residing on the resolution of (23a) and (23c). For non-recursive estimation, (23a) can be solved applying an iterative on-manifold Gauss-Newton method [34,35]. Notice that the aforementioned adjustment constitutes a local search for (23a), and thus, the initial quaternion should be carefully chosen [36]. Likewise, solution fixing (23c) also makes use of the Lie algebra, the fixed solution being updated via exponential mapping, as described in [37].
The multivariate constrained LAMBDA (MC-LAMBDA) constitutes an alternative for the procedure (23), and it can be considered as the most widespread solution to the GNSS-based attitude-only determination problem. MC-LAMBDA generally parametrizes attitude as a rotation matrix, whose nine elements are initially estimated in the float solution disregarding the orthogonality constraints (i.e., R R = I, det(R) = 1). Then, the ILS search is performed in conjunction with the evaluation of the solution fixing-which, in the MC-LAMBDA case, comprises a weighted orthogonal Procrustes problem-the estimated vector of ambiguities being the minimizer of the combined cost function [8,9,38]. Despite the popularity and proven effectiveness of MC-LAMBDA, it remains computationally demanding, and its use on recursive estimation is not as straightforward as the proposed solution based on Lie theory.

Kalman Filtering for the Joint Positioning and Attitude Estimation
Generally, the determination of the positioning and attitude solution is either separately carried out or as part of the same filtering solution, but disregarding the cross-correlation between the attitude and positioning models. This section introduces a recursive solution for the JPA problem, based on the ESKF described in Section 2.3, where the cross-correlation between positioning and attitude models is recognized and exploited for an enhanced performance.
Similarly to the positioning and attitude problems, JPA requires the determination of a vector composed by integer and on-manifold values. Thus, a three-step decomposition (float estimation, IAR, and fixing estimation) is applied like in Section 3.5. Remarkably, the IAR process does not incur modifications-other than the vector of ambiguities comprising both the positioning-and attitude-related ambiguities-and any admissible integer estimator can be applied, with the ILS being an optimal estimator [39,40]. With regards to the solution fixing, JPA proceeds similarly to the attitude case in Section 3.5, applying a WLS. The output of such a WLS is expressed as the error state (with minimal parametrization of the attitude), which is translated to the total state respecting the composition operation defined in Section 2.3. Thus, the major aspect to address in the recursive JPA problem regards the definition of the associated ESKF and the observation models, which are described in the remainder of the section.
The belief function for the state (3) is presumed to be Gaussian and conditioned on the choice for the initial state: where P t is the covariance matrix. Notice that, since the error state comprises the state uncertainty, P t ∈ R 9+n(N+1),9+n(N+1) presents minimal parametrization (i.e., the attitude represented in P t regards the uncertainty in the rotation vector θ ∈ R 3 ). The time evolution of x t is dictated by the process and measurement functions, f(·) and h(·), respectively, also known as prediction and correction models: where the process and observation noise are assumed to follow zero-mean normal distributions w t ∼ N (0, Q t ), n t ∼ N (0, Q y ), and y t is the vector of observations. The time index t is used generically: the prediction step can be performed periodically, and the correction step is triggered whenever GNSS observations are received. Thus, we use t − 1 to refer to the most recent change on the state estimation (either due to a prediction or a correction step). The ESKF operates the nonlinear equations for prediction and observation on the total statex, while the disturbances represented in the error state are assumed as small signals, meaning that second-order products are negligible and the Jacobians are easily derived.

Prediction
Step A typical dynamical model regards a vehicle to move according to the constant-velocity non-turning model [19]. Thus, the process function f(·) leads to the following linear model: where the empty spots for the matrix in (27) correspond to zero values. ∆t corresponds to the time spanned between the time at the latest prediction or correction step and the current time. Then, the covariance predictor is: with Q t = diag σ 2 p 1 3,1 , σ 2 v 1 3,1 , σ 2 θ 1 3,1 , σ 2 a 1 n(N+1),1 the diagonal matrix gathering the process noise for the position σ 2 p , velocity σ 2 v , orientation σ 2 θ , and ambiguities σ 2 a . The process noise variances constitute the tuning parameters for the filter and shall scale with the time ∆t. The Jacobian w.r.t. the error state F x is as: As discussed in Section 2.3, the presented prediction model describes the movement behavior of a non-maneuvering vehicle. A more precise prediction step would integrate inertial measurements, for which the modification of the proposed ESKF can be realized straightforwardly following [28,37].

Correction Step
At this stage, the observations y pos , y att corresponding to the positioning and attitude models are integrated into the recursive JPA estimate. For the sake of convenience, when building the Jacobian matrices, the vector of observations y JPA = y ∈ R 2n(N+1) is re-structured so that DD phase observations appear before the DD code observations such as: where we have introduced the additional notation y Φ/ρ,pos/att to refer to the set of phase and code DD observations derived from the base-master (positioning) or the slaves-master (attitude) combination. The update on the state is done according to the classical KF equations: with the measurement function h(·) relating the observations to the state estimate as explained in Section 3 and H t the Jacobian matrix for the measurement functions. Applying the chain rule on the group composition x =x ⊕ δx, H t is defined as: where Hx is the standard Jacobian of h(·) w.r.t. the nominal statex as: where the partial derivatives in Hx are indicated in a second row, A, B are as in (11), and J q ( B b j,m ) ∂(q • B b j,m • q * )/∂q q t is given in (A21). H δx is the Jacobian of the true state w.r.t. the error state as: with [q t ] L defined in (A9) and the partial derivatives indicated below the matrix. Proper stochastic modeling is key to assure the performance of the proposed JPA problem. Thus, the observations' covariance matrix Q y is defined as follows: where the subscripts used in (37) follow the notation introduced in (30). Thus, Q Φ , Q ρ ∈ R n(N+1),n(N+1) are the covariance matrices for the complete set of DD carrier and code observations. The former matrices Q Φ , Q ρ are then conformed by the set of carrier/code DD for positioning Q Φ,pos , Q ρ,pos ∈ R n,n and for attitude Q Φ,att , Q ρ,att ∈ R n×N,n×N . Matrices Q Φ,pos,att , Q ρ,pos,att ∈ R n,n×N denote the cross-correlation between the code/phase observations for positioning and attitude. Notice that conventional processing estimation schemes, which separately estimate the positioning and attitude problems, disregard the cross-correlation between the two problems (i.e., Q Φ,pos,att , Q ρ,pos,att = 0). However, this cross-correlation does exist, since the DD observations for positioning and attitude have in common the pivot satellite noise perceived at the master antenna. Notice that code and phase observations are generally assumed to be uncorrelated.
Stochastic modeling for the proposed JPA set of observations can realized straightforwardly as: where D JPA ∈ R N+1,N+1 is the differencing matrix, which introduces the base/slaves-master noise correlation, and D ∈ R n+1,n the differencing matrix to relate the pivot to the remaining satellites: Finally, W −1 Φ , W −1 ρ ∈ R n+1 are the diagonal matrices composed by the variance of the original n + 1 code and phase observations, respectively: where each σ r/1:n Φ/ρ indicates the standard deviation of the phase and code observations for the n satellites and the pivot r one. Characterization of the noise present in GNSS observables is generally realized applying satellite elevation-or carrier-to-noise density ratio C/N 0 -based models [41,42]. Recently, stochastic modeling w.r.t. the GNSS baseband signal resolution (i.e., bandwidth, modulation, autocorrelation function) has been proposed [43,44]. When the receiver characteristics and the C/N 0 model are unknown, the following elevation-based model is advised: where el i indicates the elevation of the ith satellite, a = 2 mm, b = 2 mm, and f = 100 [11]. Additionally, when considering medium baseline lengths, -i.e., when the distance to the base station exceeds 20 km-, an additional noise term is added to the positioning covariance matrices Q Φ,pos and Q ρ,pos to account for the atmospheric differences between the base and vehicle antennas [45,46]. To emphasize the importance of the JPA model to exploit the observations' cross-correlation fully, Figure 3 presents a pictorial example. Such an example considers four observations (with variances shown in the color bar) to be tracked by the base station, the master, and N = 2 slave antennas. The code DD covariance matrix Q ρ is shown in Figure 3 (left) for the separate estimation of positioning and attitude and in Figure 3 (right) for the JPA problem. One can observe that separately estimating the positioning and the attitude problems does not exploit the cross-correlation between both models, given by the noise for the pivot satellite at the master antenna.

Experimental Setup
The performance characterization of the proposed ESKF-based JPA problem was addressed for the navigation of a vessel navigation in an inland waterway channel, which is an interesting ITS application. The measurement campaign was conducted in Koblenz (Germany) on 16 May 2017 (DOY 136, UTC 09:00-14:00), the tracked vehicle the being MSBingen, a multi-purpose research vessel of the German Waterways authorities. The equipment setup listed three navXperience R 3G+C GNSS antennas connected to three separate dual frequency Javad R Delta receivers, a fiber optic gyroscope (FOG) IMU iMAR R IMU FCAI [47], and an active reflector under the master antenna. Figure 4   The trajectory followed and the location of the base station and the total station are illustrated in Figure 5, with the distance between the base and the vessel constituting a short baseline ranging between 900 m and 2.5 km. The vessel navigated from the port on the right side of Figure 5 up river and then returned back to the point of departure. Moreover, the vessel performed several turns around the three consecutive bridges nearby the total station, as shown in Figure 6. The multiple bridges were present along the navigation, especially in the surroundings of the total station position. These structures shaded the reception of satellite signals and induced multipath/non-line-of-sight (NLOS), and the navigation solution could be easily affected, leading to multiple cycle slips/loss of track to occur in the proximity of the bridges.
The reference trajectory of the vessel was obtained based on optical technology [48], using a total station on land and an active reflector mounted under the master GNSS antenna for automatic target tracking. This technology assured a positioning accuracy around one centimeter and whose error pattern was independent of GNSS, while its availability was assured even during the maneuvers realized around the bridges. The integration of the angular rates measured by the high-quality IMU was used as the benchmark solution for the attitude estimates.  The ESKF-based JPA solution was as presented in Section 4, while the positioning-and attitude-only estimates were derived from the recursive formulation of the models in Section 3. The navigation solution used GPS observations for the L1 and L2 frequencies, with a cut-off elevation of 15 • and a sampling rate of 1 Hz. Stochastic modeling of the observation was realized using the elevation-based model in (41). The IAR process used the ILS method based on search-and-shrink [3], and the navigation solution was only considered valid when the ratio-test for ambiguity resolution was passed [49,50]. Despite using the ratio-test, failure cases-i.e., when an ILS solution is considered valid, but the ambiguities are incorrectly estimated-could occur. In this experiment, the failure rates for both JPA and separately estimation of positioning and attitude were comparable and rather low (≤1% during the period of bridge passing). Thus, the fix ratio-i.e., the percentage of epochs for which the ratio-test is passed-was considered to evaluate the availability of the precise navigation estimates.

Results and Discussion
First, the performance characterization of the proposed JPA was compared against positioning and attitude determination solutions obtained in a separate manner (i.e., two Kalman filters were processed in parallel, one integrating the positioning observations and the other in charge of the attitude model). Notice that the inertial information was not integrated in any of the aforementioned filters, since the focus of the contribution was to address the gain of JPA against conventional GNSS attitude and positioning models.
The comparison between JPA and the positioning-and attitude-only models is reported in Table 1, in terms of the fix ratio over the course of the measurement campaign. The first column depicts the number of satellites n + 1 and the corresponding percentage of time for which that number was tracked. The following columns correspond to the fix ratio of the assessed positioning and attitude methods: the second column is for the proposed JPA, the third for the positioning-only, the fourth for the attitude-only, and the last column the union of the positioning and attitude (i.e., the simultaneous occurrence of having a fixed solution for the positioning-and attitude-only solution of the third and fourth columns). The first thing to notice is that, generally, the fixing ratio grew with the number of satellites, and the chances of having a fixed solution for five or less satellites was very uncommon. Some conclusions could be drawn from the comparison of JPA against separately estimating positioning and attitude: (i) The attitude model was "stronger" than the positioning one despite the nonlinearities in the observation function. There were two reasons to ground the positioning-attitude difference in performance: on the one hand, small residuals due to the atmospheric propagation delays between the vehicle and the base station might be present despite the short baseline; on the other hand, the data redundancy was higher in the attitude than in the positioning model (i.e., n(N + 1) code observations guided the float estimate of the four-dimensional unknown quaternion, while n code measurements supported the search of the three-dimensional unknown position). (ii) In average, JPA performed better than the union of separately estimating the position and attitude problems. Thus, the former provided a fixed solution for 74.60% of the time, while the latter was limited to 57.11% of time-which was a difference in precise navigation availability of over 45 min. Such increased performance was likely due to exploiting the cross-correlation between the positioning-and attitude-related observations, as discussed in Section 4.
(iii) The standalone attitude problem presented a higher fix ratio compared to the JPA, as it was unaffected by the residual atmospheric propagation delays present in the positioning problem. Thus, a practical application might be interested in executing in parallel the attitude-only and the JPA filters, leading to high availability positioning and attitude estimates and a mechanism to monitor the integrity of the algorithms if discrepancies between the estimates occurred.
Next, the positioning performance was analyzed. Besides the higher availability for the JPA against positioning-only estimates, the fixed solutions were very much equivalent up to the mm-regime. As depicted in Figure 6, a GNSS-independent localization ground truth was available, obtained using optical (laser-based) technology. For that, a total station was located on the small island in the center of the river, and automatic tracking of the active reflector located below the master antenna was performed. This area was most interesting, since the three bridges induced multipath biases and the track of satellites was often lost. As expected, the chance for having a fixed navigation solution under the bridges resulted in being null, although the navigation fix was rapidly recovered. The standard deviation for the fixed position solution was very close to a centimeter, although a few wrong fixes (with horizontal errors of up to 15 cm) occurred as well. To counteract or minimize the likelihood of wrong fixes to occur, the use of a second constellation, the integration of IMU observations, or the application of GNSS-robust techniques would be of great benefit [51,52].
Finally, the JPA-based attitude determination performance was examined. As with the positioning case, the fixed attitude solutions were equivalent between attitude-only and JPA estimates. Figure 7 depicts the estimated orientation using the Euler angles-pitch, roll, and heading-defined here as the rotation from the body to the global frames. The highest degree of similarity was achieved on the heading estimation, while the pitch and especially the roll were characterized with higher levels of noise. The better heading performance could be explained in relation to the GNSS satellite geometry-offering a higher accuracy on the horizontal plane-and the antenna configuration onboard the vehicle, as the two front antennas were coplanar with the local horizontal plane of the vehicle. It was noticeable that the IMU allowed to track subtle motions, such as the pitch variations due to the small waves in the river, and in general, the fusion of IMU and GNSS is recommended. Although the standard deviation for the heading estimates were below the degree, the roll presented some errors of up to 10 • , likely due to a wrong fix, since the period of time analyzed (12:00-13:00 UTC) corresponded to the time in which the vessel maneuvered around the bridges.

Conclusions and Future Work
This article provided a discussion of GNSS carrier phase-based positioning and attitude determination. Generally, these two estimation problems are solved individually, with some estimation machinery behind each of them, or as part of the same navigation solution, but disregarding the cross-correlation existing between the observations. The joint positioning and attitude (JPA) problem was introduced, with special emphasis on the recursive solution to the estimation problem. Thus, applying the principles of Lie theory and the well-known error state Kalman filter (ESKF), a complete recursive solution for the JPA problem was presented. Besides the simplicity of the solution, the ESKF presented as assets: (i) minimal state parametrization, avoiding possible numerical instabilities on the estimated covariance matrix; (ii) straightforward use of the three-step decomposition of the mixed integer and real parameter estimation, where computationally intensive methods for the fix solution were avoided. The proposed ESKF and JPA were motivated by the expression of the GNSS attitude model in terms of the unit quaternion rotation, which was also presented here. Unlike traditional models based on the estimate of inter-antenna baselines or rotation matrices, the model presented respected the nonlinear constraints of the rotation throughout the three consecutive least-squares solutions. Moreover, the tight search for the integer candidates in conjunction with the weighted orthogonal Procrustes problem as in the acclaimed MC-LAMBDA was bypassed.
The experimentation of this work included real data collected for an ITS application, for which a vessel navigated for several hours around an inland waterway scenario. The results indicated a clear gain in precise navigation availability for the proposed recursive JPA, when compared to filtering solutions applied in parallel for the positioning and attitude cases. However, evidence showed that, since the GNSS attitude model was clearly stronger than the positioning one (at least when three or more GNSS antennas conformed the vehicle setup), the best strategy for a real application would probably be the parallel execution of an attitude-only and JPA filtering solutions. Thus, the availability of precise attitude and positioning estimates would be maximized.
Prospective lines of research on the JPA problem relate to the theoretical characterization of possible estimators. For any estimation problem, it is fundamental to obtain the minimum achievable performance, i.e., the derivation of the associated Cramér-Rao bound (CRB). This is especially challenging for JPA, since the model constitutes a mixture of real and integer unknown parameters. Moreover, such CRB should be invariant, since the performance shall not depend on the attitude parametrization employed. Hence, the ultimate performance (i.e., efficiency) of the proposed Lie algebra-based recursive JPA solution remains an open issue until such a CRB is derived.  composition [30]. Rotations are linear operations preserving the vector length and the relative vector orientation. Thus, the rotation group SO (3) can be defined as: The rotation operation can be represented in multiple forms, although the most well known correspond to the Euler angles, rotation matrix, and quaternion. The quaternion is a four-dimensional hyper-complex number, often used to represent the orientation of a rigid body in a 3D space under the unit-norm constraint. The quaternion is defined as: q q w q u = cos (θ/2) e sin (θ/2) , with q w and q u the real and complex parts, e a unit rotation axis, and θ a rotation angle. The unit quaternion can also be expressed with the exponential mapping of the rotation vector θ ∈ R 3 , such that: q{θ} exp (θ) = cos θ 2 /2 θ θ 2 sin θ 2 /2 .

(A3)
A quaternion only describes a proper rotation under the unit-norm constraint of its components: Quaternions obeying this constraint can be said to belong to the unit-quaternion group, described by the three-sphere S 3 . The transformation quaternion to rotation matrix R is given by the following homogeneous quadratic function: R{q} = q 2 w − q u q u I 3 + 2q u q u + 2q w [q u ×], where I 3 ∈ R 3 is the identity matrix and the skew operator [v×] defines the cross-product matrix: which is a skew-symmetric matrix, i.e., [v×] = −[v×], equivalent to the matrix form of the cross product [v×] w = v × w.
The quaternion product (quaternion product is generally denoted as ⊗, which this work reserves for the Kronecker product instead.) • is non-commutative, associative, and distributive over the sum: p • q = p w q w − p u q u p w q u + q w p u + p u × q u , The composition of quaternions is bilinear and can be expressed as two matrix products: with [q] L and [q] R the left and right quaternion product matrices, respectively. These product matrices are given by: [q] L = q 0 I 4 + 0 −q u q u [q u ×] , [q] R = q 0 I 4 + 0 −q u q u −[q u ×] .
The rotation operator over a vector r(v) using quaternion parametrization is expressed as: where q * is the inverse quaternion operation:

Jacobian with Respect to the Quaternion
This subsection derives the Jacobian matrix J q (v) for the rotation of the vector v with respect to the quaternion q, defined as: Extending the quaternion multiplication and using (A7), from here, we will pay attention solely to the imaginary part: Considering the vectorial product properties: we can group q w q u × v − q w v × q u = −2q w v × q u . Substituting the vectorial product for the skew operator [·×], to facilitate the matrix formulation, and applying Equation (A17), we can reformulate the quaternion operator as: From here, deriving the partial derivatives of the quaternion results in being uncomplicated: obtaining finally: J q (v) = 2 q w v − [v×]q u −q w [v×] + q u T vI 3 + q u v − vq u ∈ R 3,4 .