Integrated GNSS Attitude Determination and Positioning for Direct Geo-Referencing

Direct geo-referencing is an efficient methodology for the fast acquisition of 3D spatial data. It requires the fusion of spatial data acquisition sensors with navigation sensors, such as Global Navigation Satellite System (GNSS) receivers. In this contribution, we consider an integrated GNSS navigation system to provide estimates of the position and attitude (orientation) of a 3D laser scanner. The proposed multi-sensor system (MSS) consists of multiple GNSS antennas rigidly mounted on the frame of a rotating laser scanner and a reference GNSS station with known coordinates. Precise GNSS navigation requires the resolution of the carrier phase ambiguities. The proposed method uses the multivariate constrained integer least-squares (MC-LAMBDA) method for the estimation of rotating frame ambiguities and attitude angles. MC-LAMBDA makes use of the known antenna geometry to strengthen the underlying attitude model and, hence, to enhance the reliability of rotating frame ambiguity resolution and attitude determination. The reliable estimation of rotating frame ambiguities is consequently utilized to enhance the relative positioning of the rotating frame with respect to the reference station. This integrated (array-aided) method improves ambiguity resolution, as well as positioning accuracy between the rotating frame and the reference station. Numerical analyses of GNSS data from a real-data campaign confirm the improved performance of the proposed method over the existing method. In particular, the integrated method yields reliable ambiguity resolution and reduces position standard deviation by a factor of about 0.8, matching the theoretical gain of 3/4 for two antennas on the rotating frame and a single antenna at the reference station.


Introduction
The acquisition and interpretation of three-dimensional (3D) spatial data are important assets for scientific and industrial applications, such as 3D city modeling, facility management, construction engineering, navigation and forensic investigations. Direct geo-referencing, which does not require dedicated control points, is an efficient methodology for the fast acquisition of 3D spatial data by means of a 3D laser scanner. It can be performed either using additional backsight targets [1][2][3][4] or using external sensors [5][6][7][8]. The latter requires the fusion of spatial data acquisition sensors and navigation sensors, such as Global Navigation Satellite System (GNSS) sensors. In this contribution, we consider an integrated GNSS navigation system to provide estimates of the position and attitude (orientation) of a 3D laser scanner.
The use of GNSS for geo-referencing has been explored in various studies. Direct geo-referencing is demonstrated using GNSS integrated with inertial sensors [9,10] and a digital compass [5]. In this work, we explore pure a GNSS-based navigation solution as in [8,11]. In [11], a single rotating antenna is used to provide a post-processing navigation solution. As in [7,8], the proposed multi-sensor system (MSS) consists of multiple GNSS antennas rigidly and symmetrically mounted on the frame of a rotating laser scanner and a reference GNSS station with known coordinates.
The proposed method uses the multivariate constrained integer least-squares (MC-LAMBDA) method [12][13][14][15][16][17] for the estimation of rotating frame ambiguities and attitude angles. MC-LAMBDA makes use of known antenna geometry to strengthen the underlying attitude model, enabling reliable instantaneous ambiguity resolution and attitude determination of the rotating frame. The reliable estimation of rotating frame ambiguities is consequently utilized to enhance the positioning of the rotating frame. This array-aided positioning method [15,[18][19][20] naturally yields the estimates of the rotating frame center (centroid of antennas' reference points) and improves ambiguity resolution, as well as the positioning accuracy of the relative position between the rotating frame and the reference station.
The numerical studies considered in this contribution include performance analyses of the proposed method with GNSS data from two real data campaigns. Comparison studies using epoch-by-epoch processing and filtering confirm the improved performance of the proposed method over the existing method from [8]. This contribution is organized as follows: Section 2 describes the multi-sensor system considered and defines the problem. Section 3 describes our attitude determination and filtering approaches for the rotating frame. Section 4 describes the array-aided positioning and filtering methods for the positioning of the rotating frame. Section 5 presents real data analyses demonstrating the improved performance of the proposed method. Finally, Section 6 contains the summary and conclusions of this contribution.

Background
The multi-sensor system (MSS) considered for geo-referencing in this contribution consists of a laser scanner and two GNSS antennas/receivers. As shown in Figure 1a, the laser scanner is the core sensor of the MSS, which rotates about its vertical axis with a constant angular velocity. GNSS receivers are connected to two eccentric GNSS antennas, which are mounted such that the centroid of the antenna reference points (ARPs) coincides with the scanner rotating axis. In addition to these GNSS receivers, it is assumed to have a nearby reference GNSS station with a known position (Figure 1b). During the data acquisition, the MSS makes a complete 360 degree rotation about its vertical axis, collecting both laser and GNSS measurements, which are synchronized through a GNSS receiver event marker. The objective of the navigation system is to provide the position (the centroid of ARPs) and the pointing direction (heading) of the laser scanner. In [8], standard real-time kinematic (RTK) positioning [21] is used to estimate individual rotating antenna positions, and then, a constrained nonlinear filtering method, in particular an extended Kalman filter, is used to obtain the above parameters. In this contribution, we use constrained integer least-squares (Section 3) and array-aided positioning (Section 4), enabling improved ambiguity resolution and improved positioning accuracy. In the following sections, we formulate a more general problem, estimating attitude angles and relative position between two platforms with multiple GNSS antennas/receivers, which enables us to demonstrate the potential of array-aided positioning. As shown in the following sections, array-aided positioning utilizes the reliable estimation of rotating frame ambiguities, which are obtained from array processing (MC-LAMBDA), improving the estimation of the relative position of the rotating frame with respect to the given reference station.

Attitude Determination
This section describes the platform processing involving attitude determination for a small-sized array of GNSS receivers/antennas with a known local body frame antenna geometry. First, the multi-baseline attitude model is introduced using the multivariate formulation of [12]. This formulation makes frequent use of the Kronecker product and the vec-operator [22]. Then, we include the local body frame antenna-geometry and show how the constrained attitude model can be solved in a step-wise manner.

The Multivariate Model
Let us consider the k-th platform equipped with a set of n k + 1 antennas simultaneously tracking m + 1 satellites on f frequencies. The set of linearized double difference (DD) GNSS phase and code observations obtained on the n k baselines formed by these antennas at an observation epoch forms a multivariate Gauss-Markov model [12,19]: where E(·) and D(·) denote the expectation and dispersion operator, ⊗ denotes the Kronecker product, . . , b k n k ] the 3 × n k matrix of n k unknown baseline vectors b j , G the 2f m × 3 geometry matrix that contains the unit line-of-sight vectors, A the 2f m × f m matrix that links the DD data to the integer ambiguities and P n k and Q yy the known matrices of order n k × n k and 2f m × 2f m, respectively. Here, vec(·) denotes the vec-operator, which transforms a matrix into a vector by stacking the columns of the matrix, one underneath the other. Note that, for the simplicity of the formulation, we assumed that all receivers/antennas track the same set of satellites. However, this restriction is relaxed in the software implemented using MATLAB. Since the unit line-of-sight vectors of two antennas to the same satellite on a short baseline considered in this work (≤10 km) are the same for all practical purposes, the geometry matrix G is considered the same for different platforms, as well as for the between-platform baseline at a given time instant.
For the stochastic model, we assumed that all receivers/antennas have similar (noise) characteristics. However, the results in the following are also applicable for dissimilar receivers/antennas [19]. The correlation matrix P n k takes care of the correlation that follows from the fact that the n k baselines share the observations from the reference receiver. For similar receivers/antennas, it is given as: with I n k the identity matrix of size n k and e n k the n k -vector of ones. Matrix Q yy takes care of the precision of the phase and code data and is given as: where: the single difference operator, "blockdiag" referring to the block diagonal matrix formed by given arguments and "diag" referring to the diagonal matrix formed by given arguments. The factor two in Equation (5) is due to the between-receiver difference of similar receivers. We assume elevation-dependent noise characteristics [23] for the undifferenced observables. That is, the standard deviation of the undifferenced observable can be written as: where θ s is the elevation angle of satellite s and σ f :t 0 , a f :t 0 and θ f :t 0 are the elevation-dependent model parameters.

The Body-Frame Antenna-Geometry as Multivariate Constraints
The strength of the above model can be improved by including information about the geometry of the antenna configuration. The known body-frame antenna-geometry can be included into the above model through the parametrization: describing the known geometry of the antenna configuration in the body frame. Here, q k is the degree of geometrical independence of the GNSS baselines, for example, q k = 1 for co-linearly installed antennas, q k = 2 for co-planarly installed antennas and q k = 3 for antennas not installed in a single plane. For q k = 3, R k is related to the Euler attitude angles ϑ = [φ θ ψ] T as follows: with φ the heading, θ the elevation, ψ the bank and where s α = sin(α) and c α = cos(α). Note that for q < 3, only the first q columns of R are defined. For example, for a linear antenna array (q = 1), only the first column is defined, and hence, only heading and elevation are estimable. For q > 1 (an array with more than two antennas that are not in a straight line), all three angles are estimable. Substitution of Equation (7) into Equation (1) leads to the constrained GNSS attitude model [19,24]: Our objective is to solve for the attitude matrix R k in a least-squares sense, thereby taking the integer constraint on matrix Z k ∈ Z f m×n k and the orthonormality constraint on matrix R k ∈ O 3×q k into account. Hence, the least-squares minimization problem that will be solved reads: with || · || 2 Q = (·) T Q −1 (·). This is a mixed integer nonlinear least-squares problem that does not permit a closed-form solution. We now describe how Equation (11) can be solved.

The Real-Valued Float Solution
The float solution is defined as the solution of Equation (11) without the constraints. When we ignore the integer constraint on Z k and the orthonormality constraint on R k , the float solutionsẐ k andR k and their variance-covariance matrices are obtained from solving the system of normal equations: with: The Z k -constrained solution of R k and its variance-covariance matrix can be obtained from the float solution as follows: Using the above estimators, the original problem in Equation (11) can be decomposed as: being the matrix of least-squares residuals. Note that the first term on the right-hand side is constant, as it does not depend on the unknown matrices Z k and R k .

The Integer Ambiguity Solution
Based on the orthogonal decomposition (15), the multivariate constrained integer minimization can be formulated as:Ž where: with:Ř The ambiguity objective function C k (Z k ) is the sum of two coupled terms: the first weighs the distance from the float ambiguity matrixẐ k to the nearest integer matrix Z k in the metric of QẐkẐk, while the second weighs the distance from the conditional float solutionR k (Z k ) to the nearest orthonormal matrix R k in the metric of QRk (Z k )R k (Z k ) . Unlike with the standard LAMBDA method [25], the search space of the above integer minimization problem is non-ellipsoidal, due to the presence of the second term in C k (Z k ). This second term is a consequence of having the orthonormality constraints rigorously included. The evaluation of C k (Z k ) requires the computation of a nonlinear constrained least-squares problem (18) for every integer matrix in the search space. In the MC-LAMBDA method, this problem is mitigated through the use of easy-to-evaluate bounding functions [17].

The Ambiguity Resolved Attitude Solution
Finally, we obtain the integer ambiguity-resolved attitude solution by substitutingŽ k into Equation (13), thus givingR k (Ž k ). The sought-for attitude angles ϑ k Ž k are then given by the reparametrized solution of Equation (18). Using a first order approximation, the formal variance-covariance matrix of the attitude angles is given by: where J R k ,ϑ k is the Jacobian of ϑ k (R k ).

Attitude Filtering
The epoch-by-epoch MC-LAMBDA attitude solution is further processed using an unscented Kalman filter (UKF) [26]. For a leveled platform (i.e., for small θ and ψ), the kinematic equations of the attitude angles are given as [27]: where the state vector α i = φ iφi θ iθi ψ iψi T consists of attitude angles and angular rates, and the state transition matrix F is given as: where T is the sampling interval. The process noise v α i−1 has a zero mean normal distribution with variance-covariance matrix Q v α v α ,i−1 , which is given as: with σ φ , σ θ and σ ψ the process noise standard deviations. The observation model reads: with ζ i given by R k (Ž k ) at epoch i. The nonlinear observational function h(α i ) is defined by Equation (8), and the observation noise w α i is assumed to have a zero mean normal distribution with covariance matrix Q w α w α ,i , which is given by The use of the above constant-velocity model [28] reflects the fact that the frame is rotating at a constant rate. For the two-antenna scenario considered in real-data analyses (Figure 2), only heading and elevation angles are estimable. Hence, a reduced state space model consisting of only heading, elevation and their rates is used in Section 5. The recursive filter is initialized with two-point initialization [28] and propagated with process noise standard deviations of σ φ = 0.01 • s − 3 2 and σ θ = 0 • s − 3 2 (i.e., dead reckoning filtering for elevation constraining to horizontal 1D rotation).

Integrated Positioning
This section describes the between-platform processing involving relative positioning between two platforms equipped with arrays of GNSS receivers/antennas. The array-aided positioning described in the following is a novel positioning concept improving between-platform positioning using an array of antennas on the platforms [15,[18][19][20]. Unlike the formulations in [18][19][20], the formulation in this contribution explicitly considers different numbers of antennas on the reference and user platforms. Unlike the parameter space formulation in [15], the current contribution considers a simplified, double-difference observation space formulation elegantly demonstrating the advantages of array-aided positioning. First, the combined observation model for all independent baselines among all receivers on both platforms is described. Then, we describe attitude-bootstrapping, showing how platform arrays improve the between-platform baseline estimate.
Let us consider two platforms carrying n 1 + 1 and n 2 + 1 receivers/antennas. The functional and stochastic models for the between-platform baseline formed by the first antennas (pivot antennas) read: E(y 12 ) = Az 12 + Gb 12 z 12 ∈ Z f m (24) D(y 12 ) = Q yy (25) where y 12 is the between-platform double-difference observables, z 12 is the unknown between-platform double-difference ambiguities and b 12 is the unknown between-platform baseline. Note that atmosphere delays are not considered in this formulation, as troposphere delays and ionosphere delays can be ignored for the short baseline (<10 km) considered in this work. However, these atmosphere delays must be taken into account for general long baseline scenarios [19]. In standard positioning, the LAMBDA method yields the optimal estimates for the ambiguities and, hence, for the baseline.

Integrated Between-Platform Model
By combining between-platform observables in Equation (24) and platform array observables in Equation (9), the functional and stochastic models of the integrated system read: is the combined rotation matrices and between-platform baseline, Z = [Z 1 Z 2 z 12 ] ∈ Z f m×nt is the combined ambiguity matrix, B 0 = blockdiag (B 1 0 , B 2 0 , 1) is the combined local geometry matrix and: is the combined correlation matrix with n t = n 1 +n 2 +1. The above system consists of attitude models of both platforms with unknowns Z k and R k and the between-platform baseline model with unknowns z 12 and b 12 . Even though these three subsystems do not have any parameter in common, they are correlated as in Equation (28), due to the use of common observations from pivot antennas.

Attitude Bootstrapping
The attitude bootstrapping method [18,19] uses a decorrelation technique to decouple the combined system in Equation (26), such that the subsystems still yield the optimal solution. Using decorrelation matrix: The decorrelated system reads: is the combined rotation matrices and between-platform baseline after decorrelation, Z = Z 1 Z 2 z 12 is the combined ambiguity matrix after decorrelation and: with: This decorrelation keeps the platform processing intact as in Equation (16)  where: Due to the reduction of variance-covariance by a factor of η, this model yields improved ambiguity resolution and baseline estimation compared to the standard positioning model in Equation (24). That is, the use of array-aided positioning reduces the variance-covariance matrices of the float ambiguities and ambiguity-fixed baseline estimators by a factor of η. For the rotating frame with two antennas and a single antenna at the reference station in Figure 1, the variance reduction factor is η = 3 4 . Note that the between-platform baseline estimate in Equation (35) corresponds to between-centroids of antenna arrays, naturally yielding the parameter of interest for the geo-referencing system in Figure 1. The unconstrained mixed-integer least-squares problem defined in Equations (37) and (38) can be solved efficiently using the LAMBDA method [25] providing ambiguity-fixed baseline estimateb 12 (Ž ) and associated variance-covariance matrix Qb 12 (Ž )b 12 (Ž ) .

Baseline Filtering
The epoch-by-epoch baseline solution in Section 4.2 is further processed to obtain the filtered estimates for the center of the MSS (Figure 1a). Unlike the previous method in [8], which uses constrained nonlinear filtering for antenna positions, the integrated method in Section 4.2 yields estimates of the center position, which is assumed to be stationary for a rotating, leveled frame. Hence, dead reckoning (linear) Kalman filtering yields the filtered estimates for the stationary center position. The kinematic equation reads: where the state vector β i = b 12 i consists of position components. The process noise v β i−1 has a zero mean normal distribution with variance-covariance matrix Q v β v β ,i−1 , which is given as: with σ bx , σ by and σ bz the process noise standard deviations. Since the center position of the rotating frame is stationary, a dead reckoning filter is used (i.e., σ bx = σ by = σ bz = 0, which is equivalent to the recursive least-squares estimation of constant parameter vector b 12 ). The observation model reads: with ξ i given byb 12 (Ž ) at epoch i. Observation noise w β i is assumed to have a zero mean normal distribution with covariance matrix Q w β w β ,i , which is given by Qb 12 (Ž )b 12 (Ž ) at epoch i.

Analyses
For numerical analyses, we used the data from a static and a kinematic experiment on the roof of the Geodetic Institute (Messdach) building, Leibniz Universität Hannover, Germany. The MSS is mounted on Pillar 5 (cf. Figure 2) and equipped with a terrestrial laser scanner (TLS) Z+F Imager 5006 and two LEIAX1202GG GNSS antennas about 0.6 m apart. These antennas are connected to two JAVAD TRE G3TH DELTA GNSS receivers. The reference station is located on Pillar 8 (about 20 m from the MSS and equipped with a JAVAD TRE G3TH DELTA GNSS receiver and a LEIAR25.R3 LEIT antenna. For the kinematic experiment, we also considered another reference station equipped with a LEICA GRX1200 GNSS receiver and a LEIAR25.R4 LEIT antenna and located about 6 km away from the MSS. The static experiment was conducted on 4 October 2011, for about six hours with the collection of GPS data at a rate of 1 Hz. In the kinematic experiment on 7 October 2011, we collected GPS data for five subsequent rotations (about an hour) at a rate of 10 Hz. For all of our analyses, we used the elevation-dependent stochastic model with the parameters given in Table 1. Table 1. Elevation-dependent stochastic model parameters defined in Equation (6).

Static Data
This section presents the analyses of the static data demonstrating the benefits of using the knowledge of antenna geometry for attitude determination and array-aided positioning. Figure 3 shows the satellite visibility: the skyplot, the number of satellites and the position dilution of precision (PDOP) values. With a 10 • elevation cut-off, we have a moderate GPS satellite geometry (PDOP < 3) during the experiment.   Table 2 summarizes the empirical instantaneous integer ambiguity resolution success rate for attitude determination, indicating the advantage of using MC-LAMBDA for the case of single-frequency processing. Similarly, Table 3 demonstrates the improved success rate performance of the proposed array-aided positioning with two antennas/receivers on the frame. Note that further improvement is possible by having more antennas/receivers. Table 2. Empirical instantaneous ambiguity resolution success rate (%) for attitude determination using static data (0.6-m baseline).

Processing Standard Positioning Array-Aided Positioning
Single-frequency 85.8 90.0 Dual-frequency 100 100 Figure 4 depicts the scatter plot of the ambiguity-fixed attitude angles, while Figure 5 shows the plots (scatter plot of the horizontal components and time series of the down component) for the ambiguity-fixed baseline estimates. Tables 4 and 5 summarise the corresponding empirical standard deviations. Note that, once the ambiguities have been resolved, the precision of the attitude solution is driven by the high precision of the carrier phase observations [16]. That is, the accuracy of the unconstrained attitude solution (using the LAMBDA method) is comparable to that of the constrained solution (using the MC-LAMBDA method), provided that ambiguities are correctly fixed. However, the knowledge of the antenna geometry plays an important role by strengthening the underlying model and, hence, improving the ambiguity resolution (see Table 2). In the case of baseline estimation (Table 5), the proposed method yields slightly improved estimates, due to the integrated processing with an array of antennas.

Kinematic Data
This section presents the analyses of the kinematic data comparing the proposed method with the existing method. These dual-frequency GPS data analyses aim to compare the estimation of the parameters of interest for geo-referencing, namely the heading and center point of the rotating frame. Figure 6 shows the satellite visibility: the skyplot, the number of satellites and the PDOP values. With a 10 • elevation cut-off, we have a good GPS satellite geometry (PDOP ≈ 2) during the experiment. Table 6 summarizes the root mean square (RMSE) values of the heading angle for all five rotations, where the ground truth is determined using the fact that the frame was rotating at a constant rate and synchronized though a GNSS receiver event marker. Since the precision of ambiguity-fixed attitude angles is driven by the high precision of the carrier phase observations, both the proposed method and the previous method [8] have a similar angular accuracy. Filtering further improves the estimates. Based on these analyses, the achievable heading angular accuracy using a 0.6-m baseline on the rotating frame is about 0.2 • RMSE.   Figure 7, while the 3D position RMSE values for all five rotations are reported in Table 7. Here, smoothing estimates based on all five rotations are considered as the ground truth. The apparent improved performance of the proposed method is due to the use of novel integrated processing. Since the proposed integrated method naturally yields the estimates of the center point (see Section 4.2), the proposed filtering has a simple and strong dynamic model compared to the nonlinear, constrained filtering of the existing method [8] and, hence, yields improved estimates. Based on these analyses, the achievable position accuracy using a 20-m baseline with two antennas on the rotating frame is about 3 mm RMSE. Finally, we considered the determination of the center point using a reference station at about 6 km away. The baseline estimation results for this long baseline using the proposed method are provided in Table 8. The significant increase of the position RMSE of this baseline compared to the short baseline case in Table 7 is due to the presence of residual atmosphere delays for this long baseline. Hence, the achievable position RMSE for this baseline is about 20 mm.

Summary and Conclusions
In this contribution, we explored the use of an array of GNSS antennas for attitude determination and relative positioning for direct geo-referencing. The MC-LAMBDA method exploits the known antenna geometry to improve the reliability of resolving rotating frame ambiguities and, hence, to improve the reliability of the rotating frame attitude determination. Furthermore, the reliable estimation of rotating frame ambiguities enables the strengthening of the estimation of the baseline between the rotating frame and a reference station. Our analysis includes epoch-by-epoch processing, as well as recursive filtering. We demonstrated the improved performance of the proposed method using data from two experiments with a prototype MSS representing a rotating frame. The use of constrained attitude determination and array-aided positioning increases the reliability (in terms of ambiguity resolution) and improves the achievable position accuracy. It enables instantaneous ambiguity resolution, which is immune to cycle slips, and, hence, enables instantaneous mapping. Furthermore, the reliability and accuracy can further be improved by employing more antennas on the rotating frame and at the reference station. With a sufficient number of low-cost GNSS receivers, the potential of instantaneous mobile mapping for low-cost applications will be explored in future studies.