# Analytical and Iterative Solutions to GNSS Attitude Determination Problem in Measurement Domain

Authors: Baowei Chen; Guobin Chang; Shengquan Li; Kailiang Deng
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

Attitude determination using double-differenced GNSS carrier phase measurements is studied. A realistic stochastic model is employed to take the correlations among the double-differenced measurements into full consideration. Two important issues concerning iteratively solving the nonlinear least-squares attitude determination problem are treated, namely, the initial guess and the iteration scheme. An analytical and sub-optimal solution is employed to provide the initial guess. In this solution, the orthogonal and determinant constraints among the elements of the direction cosine matrix (DCM) of the attitude are firstly ignored, and hence a relaxed 3×3 matrix is estimated using the linear weighted least-squares method. Then a mathematically feasible DCM, i.e., orthogonal and with +1 determinant, is extracted from the relaxed matrix estimate, optimally in the sense of minimum Frobenius norm. This analytical initial guess estimation method can be used for all feasible cases, including some generated ones, e.g., the case with only 3 antennas and only 3 satellites, subject possibly to some necessary, yet minor modifications. In each iteration, an error attitude, whose DCM is parameterized using the Gibbs vector, is introduced to relate the previously estimated and the true DCM. By linearizing the measurement model at the zero Gibbs vector, the least-squares estimate of the Gibbs vector is obtained and then used to correct the previously estimated DCM. By repeating this process, the truly least-squares estimate of the attitude can be achieved progressively. These are in fact Gauss-Newton iterations. For the final estimate, the variance covariance matrix (VCM) of the attitude estimation error can be retained to evaluate or predict the estimation accuracy. The extraction of the widely used roll-pitch-yaw angles and the VCM of their additive estimation errors from the final solution is also presented. Numerical experiments are conducted to check the performance of the developed theory. For the case with 3 2-meter long and orthogonally mounted baselines, 5 visible satellites, and 5-millimeter standard deviations of the carrier phase measurements, the root mean squared errors (RMSE) of the roll-pitch-yaw angles in the analytical solution are well below 0.5 degrees, and the estimates converge after only one iteration, with all three RMSEs below 0.2 degrees.


---

## Body

## 1. Introduction

Determination/estimation of the attitude, i.e., the rotational information between a vehicle’s body frame and a reference frame, is of significant importance in navigation, guidance, and control [1]. In general, the methods can be roughly grouped into two categories, i.e., the dynamic or filtering ones and the static or point ones [2], in both of which global navigation satellite systems can be involved [3–5]. Since about 1990s, the GNSS attitude determination system or the GNSS attitude sensor emerged as an attractive alternative [6]. It uses carrier phase signals differenced across more than 2 baselines. It is attractive because that the accuracy is moderate and hence satisfactory in many situations, that it can employ directly off-the-shelf antennas and receivers and hence is low-cost and fast-to-build, and that it is drift-free and hence suitable for long-term use. In recent years, it is receiving more and more attention in space/air/sea/land applications [7–10]. It can be used in both the dynamic and the static methods; however, in this contribution, the latter is the focus.As only the fractional parts of the carrier phases are measured, accurate pseudo range type information can be extracted only when the integer cycles are resolved with sufficiently high success rate. In GNSS attitude determination, besides the attitude calculation itself, integer ambiguity resolution is an important issue [11–14]. Though it is claimed that the two sub-problems, i.e., the ambiguity resolution and the attitude calculation, should be combined, there are still many reports in the literature which are focused on the latter only. This is worthy in twofold. First it is worthy in attitude determination itself because the integers remain known once they are resolved without cycle slips. Second, it can help to resolve the yet ambiguous carrier phase measurements. As is shown in the sequel, the developed method can determine uniquely the attitude with un-ambiguous carrier phase signals from only 3 visible satellites. In fact, it can be shown that, for dedicated receiver, only 2 visible and un-ambiguous satellites are needed. The determined attitude, as extra information, can hopefully help to fix the integer cycles of the yet ambiguous signals from other satellites. In this contribution, only the attitude calculation is studied assuming that the integer ambiguities have been resolved.The attitude determination problem using the un-ambiguous carrier phase measurements is often solved using the least-squares method, or the maximum likelihood method under the Gaussian assumption [15]. The solution is often obtained iteratively due to the inherent nonlinearity of the problem; see, e.g., [16]. Analytical solutions were also reported in the literature; however, they cannot be strictly least-squares optimal in general. By ignoring the constraints among the 9 elements of the attitude DCM, a relaxed attitude matrix, not necessarily orthogonal, was estimated in [17]. The relaxation apparently introduces approximation or in other words costs the optimality. In fact, a mathematically feasible DCM can be extracted from an arbitrary 3×3 matrix by fulfilling certain optimality criterion; see, e.g., [18]. The problem can be transformed into Wahba’s problem [19] to which analytical solutions exist; see, e.g., [20, 21]. However, the transformation is only approximate, because the weight matrix therein fails to be the inverse of the variance-covariance matrix (VCM) of the measurement errors [6, 22]. The problem is solved analytically using the total least-squares method in [23]. The carrier phase measurements and the sightlines are equally weighted; however, the latter is much more accurate than the former. For orthogonally mounted 3 equal-length baselines, and under the assumption that 3 differenced carrier phase measurement errors are of same accuracies and un-correlated, it was found that the problem becomes Wahba’s problem [24]. However, the un-correlation assumption cannot hold in general, because different difference measurements are constructed using the same reference station which clearly introduces correlations [15]. Besides, as single differenced measurements are used, the method can only be applicable for dedicated systems [25–27]. Optimal analytical solution can exist only for certain configurations of the baselines [28]. Algorithms for GNSS attitude determination are reviewed recently in [29]. Analytical methods, though sub-optimal, have merits not only in themselves as providing fast solutions, but also in that they can provide efficient initial guesses for the iterative methods.GNSS attitude determination problem can be formulated in baseline domain, also called a vectorization approach [30, 31], or in measurement domain [8]. In this contribution, the analytical and the iterative methods in measurement domain are both studied, with the former providing an initial guess for the latter. In the analytical solution, similar methodology is adopted as in [17], but the orthogonality and +1 determinant properties of the estimated attitude matrix are strictly satisfied by the following [18]. Note here that only the attitude estimate is of interest as an intermediate quantity; its error analysis as presented in [18] is not necessary. The following three differences are noted. First, a more statistically more rigorous error model for the double differenced measurements is followed as in [15]. Second, a mathematically true DCM, i.e., orthogonal and of +1 determinant, is extracted from the above estimated relaxed matrix. Third, some degenerated cases with only 2 baselines or/and only 2 sightlines are treated in detail. In the iterative solution, an error attitude DCM is defined relating the previously estimated DCM to the true DCM. This error attitude DCM, essentially a multiplicative error, is formulated using the Gibbs vector or the Rodrigues parameters. The reason behand is in that the Gibbs vector avoids trigonometric function evaluations compared to the rotation vector or Euler angles. Then the simple Gauss-Newton iteration is carried out. As, in attitude determination applications, the intuitive roll-pitch-yaw angles, one of the 12 kinds of Euler angles, are often of interest, their estimates along with the VCM of their additive estimation errors are also extracted from the final solution.In the next section, the measurement model is introduced, wherein the functional model relates the attitude to the double differenced carrier phase measurements, and a more realistic stochastic model is employed to properly take the correlations among different measurements into account. Then the least-squares problem formulation naturally followed. In Section3, an analytical, sub-optimal solution is derived. The general cases with more than 2 baselines and more than 2 line-of-sight differences (LOSD) are first treated and then the degenerated cases with 2 baselines and/or 2 LOSDs follow. In Section 4, based on the initial guess provided in the previous section, the problem is transformed to the one with the Gibbs vector of the error attitude as the arguments and solved using the simple Gauss-Newton iteration. Once the final solution, i.e., the truly least-squares estimate, is obtained, the roll-pitch-yaw angle estimates, along with the VCM of their estimation errors, are extracted. Numerical experiments are conducted in Section 5 to check the estimation accuracy of the proposed method and also the convergence property of the iteration with the proposed initial guess. The paper is concluded in Section 6.


## 2. Measurement Model and Least-Squares Problem Formulation

In this section, the measurement model of the GNSS attitude determination, including the functional model and the stochastic or error model, is introduced. For the functional model, double differenced carrier phase measurements are used indicating that the theory can be applied to non-dedicated receivers. A more realistic stochastic model is employed wherein the correlations among different measurements are taken into full consideration. Then based on the measurement model, the loss function of the least-squares method is naturally constructed which is to be minimized in the subsequent sections.Assume that carrier phase signals from totallym+1 satellites, indexed asj=0, 1, …,m, to totallyn+1 antennas, indexed ask=0, 1, …,n, are measured at an epoch. Denote the one from thejth satellite to thekth antenna as ϕkj and the corresponding measurement error as with varekj=σk,j2. Note here ekj includes all error terms remaining after the double difference operation in (1). If it cannot be assumed unbiased statistically implying the presence of the unknown systematic errors, the variance should be replaced with the mean squared error which equals the sum of the variance and the squared bias [32]. In general, we can safely assume that these measurements are independent of each other. Without loss of generality, let the 0th satellite and the 0th antenna be the reference ones, then we have the following double differenced measurement model without considering the integer ambiguities.(1)yk,j=ϕkj-ϕ0j-ϕk0+ϕ00=bkTRsj+εk,jwith(2)εk,j=ekj-e0j-ek0+e00In (1), the indices vary asj=1, 2, …,m andk=1, 2, …,n. Baseline vector bk, expressed in the body frame, is defined from thekth antenna to the 0th one. Assume the antennas are rigidly mounted on the vehicle and hence bk is constant over time and can be precisely surveyed off-line. LOSD vector sj, expressed in the reference frame, denotes the difference between the lines-of-sight from the antenna to thejth satellite and that to the 0th satellite. The LOSD can be safely assumed without errors, because the positioning and orbit errors (say about several meters) are far less than the range to the satellites (say about 20 million meters). Matrix R (3×3) is the DCM of the attitude from the reference frame to the body frame. This matrix should be mathematically orthogonal and with +1 determinant. Assume the axes of the reference frame point to north, east, and downward, respectively, and that of the body frame to forward, rightward, and downward, respectively. Under this convention, the Euler angles with the 3-2-1 axial rotation sequence are exactly the widely used roll-pitch-yaw angles. From (2), we have the following variances or covariances [15]:(3)covεk,j,εl,i=δijδklσk,j2+δijσ0,j2+δklσk,02+σ0,02whereδij denotes the Kronecker delta function. Define the following variables:(4)Y=yk,jk=1,2,⋯,n;j=1,2,⋯,mE=εk,jk=1,2,⋯,n;j=1,2,⋯,m(5)B=b1b2⋯bnS=s1s2⋯smThen we have the following matrix-form measurement model: (6)Y=BTRS+EDefine the following vectors:(7)y=vecYε=vecEr=vecRwhere vec(X) denotes the vectorization of the matrix X formed by stacking the columns of X into a single column vector. Then, according to the Kronecker product rule, we have the following vector-form measurement model:(8)y=ST⊗BTr+εwhere ⊗ denotes the Kronecker product. From (3), the VCM of the vectorized measurement error, i.e., Q=covε, can be easily constructed.Based on the measurement model (8), the following least-squares loss or cost function can be constructed:(9)μ=y-ST⊗BTrQ-1with(10)xX=xTXxThe least-squares estimate should be the one that minimizes (9) of course subject to the orthogonal and determinant constraints among the elements of the DCM. Under the Gaussian assumption about the distribution of the measurement errors, (9) becomes the minus log likelihood function (without considering some parameter-independent constants), and hence the least-squares estimate becomes the maximum likelihood one. The subsequent 2 sections are devoted to solve this problem, approximately but analytically, and rigorously but iteratively, respectively.


## 3. An Analytical Sub-Optimal Solution

In this section, the least-squares problem with the loss function as (9) is solved approximately and analytically. The resulting estimate, which cannot be rigorously least-squares in general, will be used as initial guess for the iterative method in the subsequent section. This solution is achieved by ignoring the inherent constraints among the elements of the DCM. The derivation begins with the general cases and then two degenerated cases are treated.


### 3.1. General Cases

Without considering the inherent constraints imposed onr, (8) actually represents a simple linear measurement model, or a Gauss-Markov model under the Gaussian assumption. For the general cases with no less than 9 measurements, the least-squares solution, whose loss function is as (9), can be easily obtained as follows:(11)Prr=S⊗BQ-1ST⊗BT-1(12)r-=PrrS⊗BQ-1yA 3×3 matrix R-, as an estimate of R, can be restored from r-. This estimate is not mathematically feasible, because it cannot always be orthogonal. In the sequel, a mathematically feasible DCM estimate, denoted as R^, is constructed as the nearest one to R- in the sense of the Frobenius norm. The following loss function is to be minimized:(13)φ=R-R-F=trR-R-R-R-Twhere the subscript F denotes the Frobenius norm. With the orthogonality of R in mind, (13) can be rearranged as follows.(14)φ=3-2trRR-T+trR-R-TAs only the second term in the right-hand side depends on the parameter of interest, it is equivalent to maximizing trRR-T. Fortunately analytical solution exists to this problem [33]. Perform the following singular value decomposition:(15)R-=UDVTThen the DCM estimate that minimizes (14) can be calculated as follows:(16)R^=U10001000detUdetVVTThe first-order error analysis of the estimate in (16) can be performed. However, we would not do it, because this estimate is primarily used as an initial guess for the iteration in the sequel to get more accurate estimate, and the VCM of the error in this estimate, actually also in every intermediate estimate, is not involved at all in the iteration. Only the error analysis of the final estimate, i.e., that in the last iteration, is necessary.In order to avoid being underdetermined, the method in this subsection necessitates thatmn be equal to or larger than 9. Also note that bothm andn should be equal to or larger than 2 in order to uniquely define the reference and the body frames, respectively. Accordingly, some feasible number sets of the baselines/antennas and LOSDs/satellites for the method in this subsection can be summarized as in Table 1. However, the following two are noted concerning Table 1. First, for all the sets in Table 1, the baselines should further be linearly independent, and the same goes for the LOSDs. Second, it is only for the analytical in this subsection. Actually at least 3 non-collinear antennas, i.e., 2 non-parallel baselines, and 3 non-collinear satellites, i.e., 2 non-parallel LOSDs, are sufficient to uniquely determine the attitude. Some degenerated cases, i.e., with less numbers of antennas and/or satellites as in Table 1, are treated in the following two subsections.Table 1
Feasible numbers of satellites and antennas for the analytical solution in Section3.1.
Number of antennas (n+1) Number of satellites (m+1) = 3 ⩾6 = 4 ⩾4 ⩾5 ⩾3


### 3.2. Degenerated Cases with 2 Baselines or 2 LOSDs

In this subsection, the following two degenerated cases are treated, i.e., 2 baselines with 3 or 4 LOSDs and 2 LOSDs with 3 or 4 baselines. The treatments of the two cases are the same, so without loss of generality, we take the former as the example. From the LOSDs, select arbitrarily 3 non-coplanar ones; says1, s2, and s3. Define the following variables: (17)S=s1s2s3(18)yk=yk,1yk,2yk,3T(19)εk=ek,1ek,2ek,3TThen we have the following measurement equation:(20)yk=STRbk+εkNote that, in (18), (19), and (20),k=1 and 2. For an arbitrary non-singular 3×3 matrix X and an arbitrary 3×1 vectorx, the following can be validated: (21)Xx×=detXX-Tx×X-1where [x×] denotes the cross-product matrix or anti-symmetric/skew-symmetric matrix of x. Using (21), we have the following:(22)STRb1×STRb2=STRb1×STRb2=detSTRS-1Rb1×RTS-TSTRb2=detSS-1Rb1×b2=detSS-1Rb1×b2=detSS-1Rb3So we have the following pseudo measurement model: (23)detSS-1Rb3=y1-ε1×y2-ε2≈y1×y2+y2×ε1-y1×ε2=y3-ε3or(24)y3=detSS-1Rb3+ε3where the approximation is introduced by omitting the cross product of the two error vectors. Note that the simultaneous equations (20) and (24), with vec[R] as the arguments, are of full column rank, so they are well-determined and hence can be solved. Also note that as the (pseudo) measurement error vector in (24) is the linear combination of the two in (20), the VCM of the overall measurement error vector is singular. Singular VCM implies equality constraints. After the estimate of vec[R] is obtained, a DCM estimate can be extracted from it as in the previous subsection.


### 3.3. Degenerated Case with 2 Baselines and 2 LOSDs

The case with only 2 baselines and only 2 LOSDs is of both theoretical and practice importance. It is the minimum requirement for a unique 3-dimensional attitude estimate. It also represents the most challenging situation in which the GNSS attitude sensor can work. For this case, there are only 4 raw measurements, much less than the required 9 in the algorithm for the general case. It is necessary and also possible as shown in the sequel, to construct 5 more (pseudo) measurements in order to get an analytical solution.First, let us orthonormalize both the baselines and the LOSDs through the Gram-Schmidt process, as follows [34]: (25)b-1=1b1b1=k1b1b-2=1b2-b2Tb-1·b-1b2-b2Tb-1·b-1=k21b1+k22b2(26)s-1=1s1s1=κ1s1s-2=1s2-s2Ts-1·s-1s2-s2Ts-1·s-1=κ21s1+κ22s2where x denotes the length of x. Then we have the following pseudo measurement equations:(27)cos⁡θ11=b-1TRs-1=k1κ1b1TRs1=k1κ1y1,1-k1κ1e1,1=y-1,1-e-1,1(28)cos⁡θ21=b-2TRs-1=k21κ1b1TRs1+k22κ1b2TRs1=k21κ1y1,1+k22κ1y2,1-k21κ1e1,1-k22κ1e2,1=y-2,1-e-2,1(29)cos⁡θ12=b-1TRs-2=κ21k1b1TRs1+κ22k1b1TRs2=κ21k1y1,1+κ221ky1,2-κ21k1e1,1-κ22k1e1,2=y-1,2-e-1,2(30)cos⁡θ22=b-2TRs-2=k21κ21b1TRs1+k22κ21b2TRs1+k21κ22b1TRs2+k22κ22b2TRs2=k21κ21y1,1+k22κ21y2,1+k21κ22y1,2+k22κ22y2,2-k21κ21e1,1-k22κ21e2,1-k21κ22e1,2-k22κ22e2,2=y-2,2-e-2,2Then the following 5 pseudo measurement equations can be constructed:(31)cos⁡θ3j=b-3TRs-j=b-1×b-2TRs-j=1-cos2θ1j-cos2θ2j≈1-y-1,j2-y-2,j2+e-1,j2+e-2,j21-y-1,j2-y-2,j2=y-3,j-e-3,jforj= 1 and 2. (32)cos⁡θk3=b-kTRs-3=b-kTRs-1×s-2=1-cos2θk1-cos2θk2≈1-y-k,12-y-k,22+e-k,12+e-k,221-y-k,12-y-k,22=y-k,3-e-k,3fork= 1, 2, and 3. Now we have 9 linear independent equations, and hence a least-squares estimate of vec[R] can be obtained. After the estimate of vec[R] is obtained, a DCM estimate can be extracted from it as in Section 3.1.


## 3.1. General Cases

Without considering the inherent constraints imposed onr, (8) actually represents a simple linear measurement model, or a Gauss-Markov model under the Gaussian assumption. For the general cases with no less than 9 measurements, the least-squares solution, whose loss function is as (9), can be easily obtained as follows:(11)Prr=S⊗BQ-1ST⊗BT-1(12)r-=PrrS⊗BQ-1yA 3×3 matrix R-, as an estimate of R, can be restored from r-. This estimate is not mathematically feasible, because it cannot always be orthogonal. In the sequel, a mathematically feasible DCM estimate, denoted as R^, is constructed as the nearest one to R- in the sense of the Frobenius norm. The following loss function is to be minimized:(13)φ=R-R-F=trR-R-R-R-Twhere the subscript F denotes the Frobenius norm. With the orthogonality of R in mind, (13) can be rearranged as follows.(14)φ=3-2trRR-T+trR-R-TAs only the second term in the right-hand side depends on the parameter of interest, it is equivalent to maximizing trRR-T. Fortunately analytical solution exists to this problem [33]. Perform the following singular value decomposition:(15)R-=UDVTThen the DCM estimate that minimizes (14) can be calculated as follows:(16)R^=U10001000detUdetVVTThe first-order error analysis of the estimate in (16) can be performed. However, we would not do it, because this estimate is primarily used as an initial guess for the iteration in the sequel to get more accurate estimate, and the VCM of the error in this estimate, actually also in every intermediate estimate, is not involved at all in the iteration. Only the error analysis of the final estimate, i.e., that in the last iteration, is necessary.In order to avoid being underdetermined, the method in this subsection necessitates thatmn be equal to or larger than 9. Also note that bothm andn should be equal to or larger than 2 in order to uniquely define the reference and the body frames, respectively. Accordingly, some feasible number sets of the baselines/antennas and LOSDs/satellites for the method in this subsection can be summarized as in Table 1. However, the following two are noted concerning Table 1. First, for all the sets in Table 1, the baselines should further be linearly independent, and the same goes for the LOSDs. Second, it is only for the analytical in this subsection. Actually at least 3 non-collinear antennas, i.e., 2 non-parallel baselines, and 3 non-collinear satellites, i.e., 2 non-parallel LOSDs, are sufficient to uniquely determine the attitude. Some degenerated cases, i.e., with less numbers of antennas and/or satellites as in Table 1, are treated in the following two subsections.Table 1
Feasible numbers of satellites and antennas for the analytical solution in Section3.1.
Number of antennas (n+1) Number of satellites (m+1) = 3 ⩾6 = 4 ⩾4 ⩾5 ⩾3


## 3.2. Degenerated Cases with 2 Baselines or 2 LOSDs

In this subsection, the following two degenerated cases are treated, i.e., 2 baselines with 3 or 4 LOSDs and 2 LOSDs with 3 or 4 baselines. The treatments of the two cases are the same, so without loss of generality, we take the former as the example. From the LOSDs, select arbitrarily 3 non-coplanar ones; says1, s2, and s3. Define the following variables: (17)S=s1s2s3(18)yk=yk,1yk,2yk,3T(19)εk=ek,1ek,2ek,3TThen we have the following measurement equation:(20)yk=STRbk+εkNote that, in (18), (19), and (20),k=1 and 2. For an arbitrary non-singular 3×3 matrix X and an arbitrary 3×1 vectorx, the following can be validated: (21)Xx×=detXX-Tx×X-1where [x×] denotes the cross-product matrix or anti-symmetric/skew-symmetric matrix of x. Using (21), we have the following:(22)STRb1×STRb2=STRb1×STRb2=detSTRS-1Rb1×RTS-TSTRb2=detSS-1Rb1×b2=detSS-1Rb1×b2=detSS-1Rb3So we have the following pseudo measurement model: (23)detSS-1Rb3=y1-ε1×y2-ε2≈y1×y2+y2×ε1-y1×ε2=y3-ε3or(24)y3=detSS-1Rb3+ε3where the approximation is introduced by omitting the cross product of the two error vectors. Note that the simultaneous equations (20) and (24), with vec[R] as the arguments, are of full column rank, so they are well-determined and hence can be solved. Also note that as the (pseudo) measurement error vector in (24) is the linear combination of the two in (20), the VCM of the overall measurement error vector is singular. Singular VCM implies equality constraints. After the estimate of vec[R] is obtained, a DCM estimate can be extracted from it as in the previous subsection.


## 3.3. Degenerated Case with 2 Baselines and 2 LOSDs

The case with only 2 baselines and only 2 LOSDs is of both theoretical and practice importance. It is the minimum requirement for a unique 3-dimensional attitude estimate. It also represents the most challenging situation in which the GNSS attitude sensor can work. For this case, there are only 4 raw measurements, much less than the required 9 in the algorithm for the general case. It is necessary and also possible as shown in the sequel, to construct 5 more (pseudo) measurements in order to get an analytical solution.First, let us orthonormalize both the baselines and the LOSDs through the Gram-Schmidt process, as follows [34]: (25)b-1=1b1b1=k1b1b-2=1b2-b2Tb-1·b-1b2-b2Tb-1·b-1=k21b1+k22b2(26)s-1=1s1s1=κ1s1s-2=1s2-s2Ts-1·s-1s2-s2Ts-1·s-1=κ21s1+κ22s2where x denotes the length of x. Then we have the following pseudo measurement equations:(27)cos⁡θ11=b-1TRs-1=k1κ1b1TRs1=k1κ1y1,1-k1κ1e1,1=y-1,1-e-1,1(28)cos⁡θ21=b-2TRs-1=k21κ1b1TRs1+k22κ1b2TRs1=k21κ1y1,1+k22κ1y2,1-k21κ1e1,1-k22κ1e2,1=y-2,1-e-2,1(29)cos⁡θ12=b-1TRs-2=κ21k1b1TRs1+κ22k1b1TRs2=κ21k1y1,1+κ221ky1,2-κ21k1e1,1-κ22k1e1,2=y-1,2-e-1,2(30)cos⁡θ22=b-2TRs-2=k21κ21b1TRs1+k22κ21b2TRs1+k21κ22b1TRs2+k22κ22b2TRs2=k21κ21y1,1+k22κ21y2,1+k21κ22y1,2+k22κ22y2,2-k21κ21e1,1-k22κ21e2,1-k21κ22e1,2-k22κ22e2,2=y-2,2-e-2,2Then the following 5 pseudo measurement equations can be constructed:(31)cos⁡θ3j=b-3TRs-j=b-1×b-2TRs-j=1-cos2θ1j-cos2θ2j≈1-y-1,j2-y-2,j2+e-1,j2+e-2,j21-y-1,j2-y-2,j2=y-3,j-e-3,jforj= 1 and 2. (32)cos⁡θk3=b-kTRs-3=b-kTRs-1×s-2=1-cos2θk1-cos2θk2≈1-y-k,12-y-k,22+e-k,12+e-k,221-y-k,12-y-k,22=y-k,3-e-k,3fork= 1, 2, and 3. Now we have 9 linear independent equations, and hence a least-squares estimate of vec[R] can be obtained. After the estimate of vec[R] is obtained, a DCM estimate can be extracted from it as in Section 3.1.


## 4. Gauss-Newton Iteration for Least-Squares Optimal Solution Using Gibbs Vector

The analytical solution in the previous section is approximate and cannot be optimal in the least-squares sense in general. So it is necessary to modify or correct it to approach the truly least-squares solution. Due to the inherent nonlinearity of the problem, correction more than once, i.e., iteration, may be necessary. In this section, the iteration scheme is developed. First, an error attitude DCM, parameterized using the 3×1 Gibbs vector, is employed to relate the previous DCM estimate to the true DCM. The problem becomes a free one without constraints inherent among the parameters. Then a simple Gauss-Newton iteration is employed to solve the problem.An error attitude or a multiplicative error, with DCM being denoted asδR, is introduced to relate the true and the estimated DCM, as follows:(33)R=δRR^This error attitude, as any attitude, can be parameterized in several forms, e.g., the DCM itself as in (33), the quaternions, the rotation vectors, the Gibbs vectors/Rodrigues parameters, the modified Rodrigues parameters, and the Euler angles (12 different forms with different axial rotation sequences) [35]. In this section, the Gibbs vector is employed because first, compared to the DCM or quaternion parameterization; it avoids constraints within its elements. Second, compared to the rotation vector and the Euler angles, it avoids transcendental (triangular) function evaluations. This advantage is shared with the modified Rodrigues parameters; however, the Gibbs vector is of simpler form and involves weaker nonlinearities. Note that the only disadvantage of using the Gibbs vector, i.e., being singular with a 180-degree rotation angle, cannot be a problem here. It is because the error attitude, rather than the attitude itself, should be with small angles, much less than 180 degrees, provided that the initial guess is already with relatively good accuracy. As shown in the simulation study in next section, the errors of the initial guess are less than 0.5 degrees, which is far away from the Gibbs vector’s singularity point. With the Gibbs vector parameterization, (33) is rearranged as follows:(34)R=I3-θ×I3+θ×-1R^With this parameterization of the attitude’s DCM, the original measurement model, i.e., (6), becomes a nonlinear function with arguments being the Gibbs vector θ. The corresponding least-squares problem is solved using the Gauss-Newton iteration here, with initial guess of the Gibbs vector being the zero vector. The Gauss-Newton iteration of the least-squares problem is equivalent to first linearizing the measurement model followed by a standard (linear) weighted least-squares method. The linearization of the measurement model is essentially the linearization of the attitude parameterization (34), around the zero Gibbs vector, as follows:(35)R≈I3-2θ×R^With this approximation, the measurement model (6) becomes the following:(36)Y=BTI3-2θ×R^S+Ei.e., (37)ΔY=-2BTθ×R^S+Ewith ΔY=Y-BTR^S. The vector form of (37) is follows:(38)Δy=Jθ+εwith ∆y = vec[∆Y], and (39)J=2BTR^s1×2BTR^s2×⋮2BTR^sm×Then the least-squares solution follows naturally as follows:(40)θ^=JTQ-1J-1JTQ-1ΔyUsing this estimate, the previous DCM estimate can be updated or corrected. This linearization-estimation-correction process is repeated to progressively approach the truly least-squares solution. Note that the correction should be performed using the rigorous one, i.e., (34), rather than the approximate one, i.e., (35). This is important, because the Gibbs vector estimate cannot be seen as sufficiently small, say several tenth of degrees. After the iteration terminates, the VCM of the Gibbs vector’s estimation error can be retained in the final iteration, as follows:(41)Pθθ=JTQ-1J-1Note that this is only a first-order approximation, whose bias, along with the bias of the least-squares estimate itself, can be further corrected by considering the higher-order terms, mainly the second-order one [36]. However, this will not be treated further in this contribution, only the first-order approximation is considered. Also note that, for small errors of the final solution, and within the first-order approximation, the error Gibbs vector in the VCM (41) can also be any one of the following error terms, half the Euler angles (any one of the 12 different forms), half the rotation vector, the vector part of the quaternion, 2 times the modified Rodrigues parameters. All these parameterizations of the error attitude are of multiplicative type rather than the plain additive type, i.e., simply the additive errors in the estimate of the attitude itself. While the DCM estimate and the VCM of the error attitude in (41) are often sufficient for use in frame rotation, attitude correction, etc., in some situations, the estimate and the VCM of its additive estimation errors, of the intuitively and physical meaningful Euler angles, especially the roll-pitch-yaw angles, are of interest. Fortunately, they can be readily extracted from the solution in this subsection; see, e.g., [1, 35].


## 5. Numerical Experiments

The main objective of the simulation is to check the quality of the analytical solution and the convergence property of the iterative method with the analytical solution as its initial guess and also to check the attitude determination accuracies for varying attitudes. Without loss of generality, only the general case, which is also most common in practice, is considered.


### 5.1. Experiments Setup

Assume that 3 baselines are deployed whose coordinates in the body frame are200T, 020T, and 002T, respectively. The true pitch-roll-yaw angles vary across time as follows:(42)αt=12∘sin⁡π9tβt=10∘sin⁡π6tγt=45∘sin⁡π300tThe time span of the simulation is set as 30 minutes, in which the receiver-to-satellite sightlines are assumed constant, of course without loss of generality. The GNSS output rate is set as 1 Hz. Totally 5 receiver-to-satellite sightlines are generated randomly as follows. Generate the azimuth angleφ from a uniform distribution in interval [0 2π]. Generate the elevation angleθ from a uniform distribution [0.1π 0.5π], considering the mask angle [37]. We have the following receiver-to-satellite sightlines:(43)lj=sin⁡θjcos⁡φjsin⁡θjsin⁡φjcos⁡θjwithj=0, 1, 2, 3, 4. Then totally four LOSDs can be obtained as sj=lj-l0, withj=1, 2, 3, 4. Assume the standard deviation of the random errors (in length) in raw carrier phases is 5 mm. The following root mean squared errors (RMSE) are employed to illustrate the attitude determination accuracy, without loss of generality, with the pitch as the example:(44)RMSEi=130×60∑t=130×60dαit2wheret denotes thetth epoch, andi denotes theith iteration.


### 5.2. Results and Analyses

The RMSEs of the roll-pitch-yaw angles are depicted in Figure1. From the simulation results as shown in Figure 1, the following four can be observed. First, the analytical solution is with rather good quality, with the RMSEs below 0.5 degrees for all three channels. This is meaningful in twofold: (1) the analytical solution itself can be used as the final solution in situations when computation resource is limited; (2) good-quality analytical solution, used as initial guess for iterative method, can assure and speed the convergence. Second, the iteration converges rather fast. Actually, for all three channels, only one iteration is sufficient for the estimate to converge. Third, the errors in the iterative solution can be reduced by about 50% compared to the analytical solution in terms of RMSE, specifically from about 0.37-0.33-0.25 degrees to about 0.16-0.17-0.14 degrees for the roll-pitch-yaw angles. Fourth, the final errors in terms of RMSE for all three channels are well below 0.2 degrees. This is consistent with that reported in the literature [7].Figure 1
Mean squared errors of the estimated roll-pitch-yaw angles versus the numbers of iterations.To check the real-time performance of the proposed method, the attitude determination errors across time is depicted in Figures2, 3, and 4 for roll, pitch, and yaw, respectively. From the results about the RMSEs in the above, at every epoch, only one iteration is performed. It can be seen clearly from these figures that the errors for all three channels and at all epochs are well below 1 degree in magnitude. Actually, in most of the time, they are within 0.5 degrees in magnitude. Hopefully, the accuracy shown in this simulation, along with the theory developed above, can be used as reference in the mission design phase in practical applications.Figure 2
Roll angle estimation errors versus epochs.Figure 3
Pitch angle estimation errors versus epochs.Figure 4
Yaw angle estimation errors versus epochs.


## 5.1. Experiments Setup

Assume that 3 baselines are deployed whose coordinates in the body frame are200T, 020T, and 002T, respectively. The true pitch-roll-yaw angles vary across time as follows:(42)αt=12∘sin⁡π9tβt=10∘sin⁡π6tγt=45∘sin⁡π300tThe time span of the simulation is set as 30 minutes, in which the receiver-to-satellite sightlines are assumed constant, of course without loss of generality. The GNSS output rate is set as 1 Hz. Totally 5 receiver-to-satellite sightlines are generated randomly as follows. Generate the azimuth angleφ from a uniform distribution in interval [0 2π]. Generate the elevation angleθ from a uniform distribution [0.1π 0.5π], considering the mask angle [37]. We have the following receiver-to-satellite sightlines:(43)lj=sin⁡θjcos⁡φjsin⁡θjsin⁡φjcos⁡θjwithj=0, 1, 2, 3, 4. Then totally four LOSDs can be obtained as sj=lj-l0, withj=1, 2, 3, 4. Assume the standard deviation of the random errors (in length) in raw carrier phases is 5 mm. The following root mean squared errors (RMSE) are employed to illustrate the attitude determination accuracy, without loss of generality, with the pitch as the example:(44)RMSEi=130×60∑t=130×60dαit2wheret denotes thetth epoch, andi denotes theith iteration.


## 5.2. Results and Analyses

The RMSEs of the roll-pitch-yaw angles are depicted in Figure1. From the simulation results as shown in Figure 1, the following four can be observed. First, the analytical solution is with rather good quality, with the RMSEs below 0.5 degrees for all three channels. This is meaningful in twofold: (1) the analytical solution itself can be used as the final solution in situations when computation resource is limited; (2) good-quality analytical solution, used as initial guess for iterative method, can assure and speed the convergence. Second, the iteration converges rather fast. Actually, for all three channels, only one iteration is sufficient for the estimate to converge. Third, the errors in the iterative solution can be reduced by about 50% compared to the analytical solution in terms of RMSE, specifically from about 0.37-0.33-0.25 degrees to about 0.16-0.17-0.14 degrees for the roll-pitch-yaw angles. Fourth, the final errors in terms of RMSE for all three channels are well below 0.2 degrees. This is consistent with that reported in the literature [7].Figure 1
Mean squared errors of the estimated roll-pitch-yaw angles versus the numbers of iterations.To check the real-time performance of the proposed method, the attitude determination errors across time is depicted in Figures2, 3, and 4 for roll, pitch, and yaw, respectively. From the results about the RMSEs in the above, at every epoch, only one iteration is performed. It can be seen clearly from these figures that the errors for all three channels and at all epochs are well below 1 degree in magnitude. Actually, in most of the time, they are within 0.5 degrees in magnitude. Hopefully, the accuracy shown in this simulation, along with the theory developed above, can be used as reference in the mission design phase in practical applications.Figure 2
Roll angle estimation errors versus epochs.Figure 3
Pitch angle estimation errors versus epochs.Figure 4
Yaw angle estimation errors versus epochs.


## 6. Concluding Remarks

Attitude determination using GNSS carrier phase difference or interference technology is receiving more and more applications in recent years. The processing of the un-ambiguous signals is studied in this contribution, assuming that the integer ambiguities have been resolved with sufficiently high success rate maybe at previous epochs. First of all, the problem should be properly modeled. The functional model, relating the attitude with the double differenced carrier phase measurements, is well known. The stochastic model, which should effectively tell the statistical property of the measurement errors, is emphasized here. A realistic stochastic model is employed in order to fully take the correlations among the double differenced measurements into account. To get the least-squares estimate of the attitude, one is in fact solving a nonlinear optimization problem, maybe also a constrained one if the attitude is formulated using redundant parameters, e.g., the direction cosine matrix (DCM) used in this work. In general, no analytical solution exists, but iterations are necessary. Then two important issues arise naturally needing to be treated properly, i.e., the initial guess to start the iteration and the iteration scheme. To get the initial guess is essential to develop an analytical solution which may not be optimal. However, it should be approximately optimal or sub-optimal; otherwise the converging speed, even the convergence itself, may not be assured. In this work, an analytical, sub-optimal, two-step solution is employed to provide the initial guess. In the first step, the orthogonal and determinant constraints on the DCM are ignored, and hence the least-squares estimate of it can be readily obtained. In the second step, a mathematically feasible DCM, i.e., orthogonal and of +1 determinant, is extracted from the above relaxed estimate. The extraction is optimal in that the Frobenius norm of the difference between the two matrices is a minimum. For the iteration, the following two are noted, namely, the parameterization of the attitude and the iterating strategy. With the assumption that the initial guess or the previous estimate is already of relatively good quality, an error, or delta attitude is introduced to relate the available estimate to the unknown true attitude. This error attitude is parameterized as the Gibbs vector. The advantages are mainly twofold. First, the problem becomes an unconstrained one. Second, it avoids transcendental (trigonometric) function evaluations. The simple Gauss-Newton iteration is employed here. For this iteration, we simply first linearize the DCM with respect to the Gibbs vector at the zero vector and then perform the (linear) weighted least-squares estimation.A numerical experiment is conducted to check the performance of the developed method. In this experiment, 3 2-meter long baselines are mounted orthogonally on the vehicle. Totally 5 satellites are assumed visible, and hence there are 4 antenna-to-satellite line-of-sight differences. The standard deviation of the un-differenced carrier phase measurement errors is assumed as 5 millimeter. This error is mainly due to multi-path effect, while other errors are assumed to have been effectively canceled by double differencing. The root mean squared errors (RMSEs) are chosen as the accuracy index. From the simulation results, it is found that the RMSEs of the analytical solutions, i.e., the initial guess, for all roll, pitch, and yaw angles are well below 0.5 degrees. This shows the relatively good quality of the initial guess. A rather fast convergence is also found; more specifically, the method converges after only one iteration, bringing all three RMSEs down below 0.2 degrees. The estimation errors versus time are also depicted. It is found that the errors, for all three channels and all epochs, are uniformly below 1 degree in magnitude; in fact, in most of the time they are within ±0.5 degrees. It is hoped that the developed theory, along with simulation results, can help to analyze the feasibility of the GNSS attitude determination in some potential missions.We want to make a final remark about the significant integer ambiguity resolution issue. This issue is not considered here for which the reason has been discussed before. However, the dependencies between the attitude calculation and the integer resolution can be two-way in some sense. While the dependence of the former on the latter is well known, the reverse one is briefly discussed as follows. As is shown before, ambiguous carrier phase signals from only 3 visible satellite can ensure the unique determination of the attitude. For dedicated receivers, even 2 satellites can suffice. The determined attitude can provide useful information for the resolution of the integer cycles of the yet ambiguous signals from other satellites. Hopefully, this topic can be a potential subject of future study.



---

*Source: 7908675-2019-04-03.xml*