Rotation Matrix Method Based on Ambiguity Function for GNSS Attitude Determination

Global navigation satellite systems (GNSS) are well suited for attitude determination. In this study, we use the rotation matrix method to resolve the attitude angle. This method achieves better performance in reducing computational complexity and selecting satellites. The condition of the baseline length is combined with the ambiguity function method (AFM) to search for integer ambiguity, and it is validated in reducing the span of candidates. The noise error is always the key factor to the success rate. It is closely related to the satellite geometry model. In contrast to the AFM, the LAMBDA (Least-squares AMBiguity Decorrelation Adjustment) method gets better results in solving the relationship of the geometric model and the noise error. Although the AFM is more flexible, it is lack of analysis on this aspect. In this study, the influence of the satellite geometry model on the success rate is analyzed in detail. The computation error and the noise error are effectively treated. Not only is the flexibility of the AFM inherited, but the success rate is also increased. An experiment is conducted in a selected campus, and the performance is proved to be effective. Our results are based on simulated and real-time GNSS data and are applied on single-frequency processing, which is known as one of the challenging case of GNSS attitude determination.


Introduction
Given that global navigation satellite system (GNSS) attitudes are unaffected by drifts and do not require any alignment, GNSS are well suited for attitude determination. The attitude could be noted as yaw (ψ), pitch (θ) and roll (φ). One baseline vector composed by two antennas comprises two attitude angles like (ψ, θ). Thus, we could get three attitude angles from two baselines which cannot be arranged in a parallel frame. The basic measurement used for GNSS attitude determination is the phase difference (∆ϕ) between the signals received by two antennas. However, GNSS attitude determination includes unknown integer ambiguities, which should be solved at first. The baseline length is usually known and it will help in integer ambiguities search. Many approaches have been studied for resolving the GNSS attitude ambiguity resolution problem [1][2][3]. In this study, we use the rotation matrix method (RMM) combined with the AFM method to solve the problem. This method is efficiency used for GNSS ambiguity resolution. When the integer ambiguities are fixed, the attitude angles can be calculated on the basis of the AFM method. Although the AFM is more flexible, it is lack of analysis on the relationship of the geometric model and the noise error. In this study, the influence of the satellite geometry model on the success rate is analyzed. The computation error and the noise error are effectively treated. It is introduced in detail and the performance of this strategy is presented in the sections below. Not only the flexibility of the AFM is inherited, but also the success rate is increased.
ϕ " 1 λ pρ`δρ`cδt r´c δt s`δ ρ t´δ ρ i q´N`ε (1) where ϕ is the GNSS receiver carrier-phase observation; λ is the GNSS carrier wavelength; ρ is the range between the receiver antenna and GNSS satellite; δρ is the orbital error along the line of sight from the satellite to station; c is the speed of light; δt r is the receiver clock offset from GNSS time; δt s is the satellite clock offset from GNSS time; δρ t is the troposphere delay; δρ i is the ionosphere delay; N is the carrier-phase integer ambiguity; and ε is an error term that includes the measurement noise, multi-path errors, others. An attitude determination system based on GNSS often consists of two receivers to receive the GNSS signals from independent antennae. Given that the common clock is used in the system, the satellite clock error can be removed by a single difference (SD) [4]. For a baseline length of a few meters, the orbital and atmosphere errors in Equation (1) are actually the same, so that these errors can be removed by a single difference. However, the receiver clock error δt r still exists. The measurement of the single difference is expressed as: ∆ϕ " 1 λ p∆ρ`c∆δt r q´∆N`∆ε (2) As shown in Figure 1, the SD model is built in local level frame (LLF).
Sensors 2016, 16, 841 2 of 13 the AFM is inherited, but also the success rate is increased. The results are based on simulated and real-time GNSS data and are applied on single-frequency processing, which is known as one of the challenging case of GNSS attitude determination.

Model of GNSS Attitude Determination
The carrier-phase equation for GNSS can be written as follows: where  is the GNSS receiver carrier-phase observation;  is the GNSS carrier wavelength;  is the range between the receiver antenna and GNSS satellite;  is the orbital error along the line of sight from the satellite to station; c is the speed of light; r t  is the receiver clock offset from GNSS time; s t  is the satellite clock offset from GNSS time; t  is the troposphere delay; i  is the ionosphere delay; N is the carrier-phase integer ambiguity; and  is an error term that includes the measurement noise, multi-path errors, others. An attitude determination system based on GNSS often consists of two receivers to receive the GNSS signals from independent antennae. Given that the common clock is used in the system, the satellite clock error can be removed by a single difference (SD) [4]. For a baseline length of a few meters, the orbital and atmosphere errors in Equation (1) are actually the same, so that these errors can be removed by a single difference. However, the receiver clock error r t  still exists. The measurement of the single difference is expressed as: As shown in Figure 1, the SD model is built in local level frame (LLF). In Figure 1, b  is the baseline vector formed by the antennae A and B, which contains the attitude parameters. The SD carrier phase measurement equation is expressed as: where (  In Figure 1, Ñ b is the baseline vector formed by the antennae A and B, which contains the attitude parameters. The SD carrier phase measurement equation is expressed as: where Ñ b " |b|pcosθsinψ, cosθcosψ, sinθq is the baseline vector, ψ, θ are the yaw and pitch, respectively; Ñ s i " pcosβ i sinα i , cosβ i cosα i , sinβ i q is the satellite I vector; α i , β i are the yaw and pitch, respectively; η i is the included angle between the baseline vector and the satellite I vector; and ∆N i , ∆ϕ i are the integral and decimal part of the SD carrier measurement, respectively. After the measurement of the double difference (DD) between satellites I and J, the receiver clock error can be removed. The DD measurement is expressed as: where Ñ b " |b|pcosθsinψ, cosθcosψ, sinθq is the baseline vector, ψ, θ are the yaw and pitch, respectively; Ñ s j´Ñ s i " | Ñ s j´Ñ s i |pcosβ ij sinα ij , cosβ ij cosα ij , sinβ ij q is the difference between the satellite I and J vectors; α ij , β ij are the yaw and pitch, respectively; η ij is the included angle between the baseline vector and the satellite ( Ñ s j´Ñ s i ) vector; and ∇∆N ij , ∇∆ϕ ij are the integral and decimal part of the DD carrier measurement, respectively. As shown in Figure 2, the DD model is built in LLF.
 are the integral and decimal part of the DD carrier measurement, respectively. As shown in Figure 2, the DD model is built in LLF. The parameters that are unknown in this function are N ij  ,  , and  . Assuming that the value of N ij  is known in the initial course, we can determine the attitude information using the RMM method. Thus, the first step we should discuss is how to find out the integer ambiguities.

Rotation Matrix Method in Resolving Equations
The important feature in this method is the use of the RMM to resolve the equations problem. This method is aimed at two equations. Equation (3) is a nonlinear equation, including sine and cosine functions. The solution of the equations is given by the analytical resolution method [5,6], but the solution process is very complex and the error angles of different scales in trigonometric functions have not been analyzed. The rotation matrix method (RMM) is very convenient to obtain. It will be discussed from two aspects.

Basic Model of Space
The basic model of spatial geometry is constructed according to the spatial relationship between the baseline vector and the satellite vector, as shown in Figure 3.

In this figure, AB
 is the baseline vector; AC  is the satellite vector;  is the pitch of the satellite vector;  is the included angle between the yaw of the baseline vector and the yaw of the satellite vector; and  is the included angle between the baseline vector and the satellite vector.
The geometric relations ( cos cos cos     ) can be proved in the basic model. This model has a regular structure and a clear geometric relationship. The following operations need to be conducted is to incorporate the actual situation of the baseline vector and the satellite vector into the model. The parameters that are unknown in this function are ∇∆N ij , ψ, and θ. Assuming that the value of ∇∆N ij is known in the initial course, we can determine the attitude information using the RMM method. Thus, the first step we should discuss is how to find out the integer ambiguities.

Rotation Matrix Method in Resolving Equations
The important feature in this method is the use of the RMM to resolve the equations problem. This method is aimed at two equations. Equation (3) is a nonlinear equation, including sine and cosine functions. The solution of the equations is given by the analytical resolution method [5,6], but the solution process is very complex and the error angles of different scales in trigonometric functions have not been analyzed. The rotation matrix method (RMM) is very convenient to obtain. It will be discussed from two aspects.

Basic Model of Space
The basic model of spatial geometry is constructed according to the spatial relationship between the baseline vector and the satellite vector, as shown in Figure 3.
In this figure, Ñ AB is the baseline vector; Ñ AC is the satellite vector; β is the pitch of the satellite vector; α is the included angle between the yaw of the baseline vector and the yaw of the satellite vector; and η is the included angle between the baseline vector and the satellite vector. The geometric relations (cosη " cosαcosβ) can be proved in the basic model. This model has a regular structure and a clear geometric relationship. The following operations need to be conducted is to incorporate the actual situation of the baseline vector and the satellite vector into the model.

Generation of Rotation Matrix
Two DD equations are given below: The spatial relation of two DD equations is shown in Figure 4. Figure 4. The spatial relation of two DD equations.
, , , , First, the ACD plane determined by the AC  vector and the AD  vector is defined as the level plane of the new coordinate system. Attitude rotation is then conducted. The rotating method is applied from the local level frame (LLF) to the new coordinate system (b): vector;  is the included angle between the ACD plane and the 2 2 X Y plane. In the new coordinate system  , as shown in Figure 5.

Generation of Rotation Matrix
Two DD equations are given below: The spatial relation of two DD equations is shown in Figure 4.

Generation of Rotation Matrix
Two DD equations are given below: The spatial relation of two DD equations is shown in Figure 4. Figure 4. The spatial relation of two DD equations.
, , , , First, the ACD plane determined by the AC  vector and the AD  vector is defined as the level plane of the new coordinate system. Attitude rotation is then conducted. The rotating method is applied from the local level frame (LLF) to the new coordinate system (b): vector;  is the included angle between the ACD plane and the 2 2 X Y plane. In the new coordinate system  , as shown in Figure 5.
First, the ACD plane determined by the Ñ AC vector and the Ñ AD vector is defined as the level plane of the new coordinate system. Attitude rotation is then conducted. The rotating method is applied from the local level frame (LLF) to the new coordinate system (b): where α is the yaw of the  ). At this point, the DD equations are given by:   are the yaw and pitch of the baseline vector in ( , respectively. According to Equation (6), the results can be obtained by: Finally, the baseline vector ( ; ; x y z is converted through the rotation matrix from the new coordinate system to the local level frame (LLF). This algorithm reduces the computational complexity compared with the analytical resolution method [5,6]:

Method of Satellite Selection in Resolving Equations
According to Equation (7), the noise effect is larger when the | | value are smaller. Thus, the following settings should be applied for satellite selection:

Integer Ambiguities Determination Method Based on Constraint Conditions
Many solutions have been studied to determine integer ambiguities, and all kinds of constraint conditions exist, such as the fixed baseline or micro-electromechanical system (MEMS) that are the yaw and pitch of the baseline vector in (O´X b Y b Z b ), respectively. According to Equation (6), the results can be obtained by: Finally, the baseline vector px b ; y b ; z b q is converted through the rotation matrix from the new coordinate system to the local level frame (LLF). This algorithm reduces the computational complexity compared with the analytical resolution method [5,6]:

Method of Satellite Selection in Resolving Equations
According to Equation (7), the noise effect is larger when the | Ñ s x´Ñ s i | value and the α ik b value are smaller. Thus, the following settings should be applied for satellite selection:

Integer Ambiguities Determination Method Based on Constraint Conditions
Many solutions have been studied to determine integer ambiguities, and all kinds of constraint conditions exist, such as the fixed baseline or micro-electromechanical system (MEMS) that provides a small-angle search region [7,8]. For the DD equation, an integer ambiguity value that we should first determine, as well as its span, is described as follows: Finally, the integer candidates that pass the constraint conditions are incorporated into the equations, and the candidate solution ( Ñ b " px; y; zq) can be calculated. These candidate solutions are evaluated and distinguished by the AFM. The AFM was originally proposed by Counselman and Gourevitch and later implemented by Remondi and Mader [9]. The candidate solutions can more easily be determined using RMM directly than by searching for the maximum in the full 2D space. From n pairs of integer ambiguity candidates, m pairs of preliminary solutions exist: px 1 ; y 1 ; z 1 q, px 2 ; y 2 ; z 2 q,¨¨¨, px m ; y m ; z m q. Only one of these solutions is correct, that is, the one that passes via AFM: where ( Ñ b " px; y; zq), N is the number of satellites and (i, j) are the master satellite and the concomitant satellite, respectively. Assuming that m pairs of preliminary solutions px 1 ; y 1 ; z 1 q, px 2 ; y 2 ; z 2 q,¨¨¨, px m ; y m ; z m q exist in k epoch, the float ambiguity of one group ( Ñ b p " px p ; y p ; z p q) can be described as: Thus, the ambiguity function is: Considering the measurement noise, a threshold T near 1 is necessary to filter out the incorrect solution [5]. According to the first-order difference of the DD carrier measurement, its fluctuation range is˘0.9 cm (« 0.05 λ) as shown in Figure 6. There are (2˚0.05 λ) being introduced in Equation (6). If the noise threshold is set to 0.2 λ, the T value is 0.8: For N satellites tracked, the integer ambiguity vector is described as: where |¨| denotes a rounding calculation [5].
In the initial course, the influence of the former epochs should be reduced. The solution is given by: The M value is the memory decline factor. After a few epochs, two cases always occur: (a) only one candidate satisfies Equation (14) at epoch k, or (b) a number of solutions satisfy the threshold. The real value is then selected with the following method. Ideally, for example, the number of satellites is Sensors 2016, 16, 841 7 of 13 greater than 8, and the geometric relationship is relatively good. The optimal value is obtained when one of the solutions is twice the suboptimal value [6]. The inequality is given by: In practical application, this ratio is slightly less than 2, while the ratio remains stable. If the optimal value is greater than 0.9, and the difference between the optimal value and the suboptimal value is greater than 0.3, the optimal solution is determined after 50 epochs, as shown in Figure 7: If the computation error and the noise error are not effectively treated, the success rate of the solution will be decreased. Thus, the influence of the satellite geometry model on the error should be analyzed in detail. In contrast to the AFM, the geometric correlation methods, such as the LAMBDA method [10], the C-LAMBDA method [11][12][13][14] and the M-LAMBDA method [15], get better performance in solving the relationship of the geometric model and the noise error. This relationship about the AFM is analyzed in the sections below. The computation error and the noise error are effectively treated. Not only the flexibility of the AFM is inherited, but also the success rate is increased.
In practical application, this ratio is slightly less than 2, while the ratio remains stable. If the optimal value is greater than 0.9, and the difference between the optimal value and the suboptimal value is greater than 0.3, the optimal solution is determined after 50 epochs, as shown in Figure 7: If the computation error and the noise error are not effectively treated, the success rate of the solution will be decreased. Thus, the influence of the satellite geometry model on the error should be analyzed in detail. In contrast to the AFM, the geometric correlation methods, such as the LAMBDA method [10], the C-LAMBDA method [11][12][13][14] and the M-LAMBDA method [15], get better performance in solving the relationship of the geometric model and the noise error. This relationship about the AFM is analyzed in the sections below. The computation error and the noise error are effectively treated. Not only the flexibility of the AFM is inherited, but also the success rate is increased.
In practical application, this ratio is slightly less than 2, while the ratio remains stable. If the optimal value is greater than 0.9, and the difference between the optimal value and the suboptimal value is greater than 0.3, the optimal solution is determined after 50 epochs, as shown in Figure 7: If the computation error and the noise error are not effectively treated, the success rate of the solution will be decreased. Thus, the influence of the satellite geometry model on the error should be analyzed in detail. In contrast to the AFM, the geometric correlation methods, such as the LAMBDA method [10], the C-LAMBDA method [11][12][13][14] and the M-LAMBDA method [15], get better performance in solving the relationship of the geometric model and the noise error. This relationship about the AFM is analyzed in the sections below. The computation error and the noise error are effectively treated. Not only the flexibility of the AFM is inherited, but also the success rate is increased.

The Influence of the Noise Error on the DD Equations
For the convenience of analysis, according to Equations (7) and (9), assuming that the satellite parameters of the DD equations are set to | Ñ s j´Ñ s i | " 0.7, | Ñ s k´Ñ s i | " 0.7, α ik b " 70˝. Now noise is added to Equation (7) and the parameters are substituted into Equation (7): Assuming that the vector ( (19) can be expressed as: According to Equation (4), although the DD model eliminates the receiver clock error and the satellite clock error, the cost of the DD measurement noise root mean square error is ? 2 times than the SD measurement, which is generally about 1 cm (i.e., roughly 0.05 GPS L1wavelength). Thus Then, Equation (23) is substituted into Equation (12) and the result is expressed as: . It represents the geometric relationship between the candidate vector and the satellite vector. If the correlation of geometric relationship is smaller, the value (rmsep∇∆N im q) of the float ambiguity is smaller. Thus, we need to find the suitable satellite vector ( Ñ s m´Ñ s x ), so that the value (rmsep∇∆N im q) of the float ambiguity is the smallest. This process is equivalent to ambiguity decorrelation adjustment of the LAMBDA method [10][11][12]. If this value (rmsep∇∆N im q) is smaller, the correlation interference of the noise error is smaller and the robustness of the ambiguity function (Fpx, y, zq) is better. In contrast to the LAMBDA method, this method gets better performance in reducing computational complexity.

Comparison of the Simulation Results
The simulation is performed from two aspects. In the first aspect, the algorithm is processed with ambiguity decorrelation adjustment of the noise error, and the other one is not do it. The simulation results are shown in Figures 8 and 9. where im  is the noise error of the DD equation (the satellite vector: m i s s  ); ˆi m N  is the float ambiguity. The Equation (21) is substituted into Equation (24), the noise root mean square error ) of the float ambiguity is expressed as:

Ambiguity Decorrelation Adjustment of the Geometric Relationship
According to the previous analysis, if the basic equations (Equation (6)

Comparison of the Simulation Results
The simulation is performed from two aspects. In the first aspect, the algorithm is processed with ambiguity decorrelation adjustment of the noise error, and the other one is not do it. The simulation results are shown in Figures 8 and 9.   The performance of the algorithm in Figure 9 is better than the algorithm in Figure 8. The algorithm is effective, and the robustness of the ambiguity function ( ( , , ) F x y z ) is improved. In the second aspect, the success rate of processing and not-processing is compared, as shown in Table 1.

Experimental Attitude Determination Results
In this section the RMM method is demonstrated using data collected in an experiment. The experiment is conducted on the top of a building. In this experiment, two receivers (COMNAV K500, Shanghai Siyue Technology Co. Ltd, Shanghai, China) are used, and both are connected to the GNSS antennae. The receiver electro-circuit is shown in Figure 10.  The performance of the algorithm in Figure 9 is better than the algorithm in Figure 8. The algorithm is effective, and the robustness of the ambiguity function (Fpx, y, zq) is improved. In the second aspect, the success rate of processing and not-processing is compared, as shown in Table 1.

Experimental Attitude Determination Results
In this section the RMM method is demonstrated using data collected in an experiment. The experiment is conducted on the top of a building. In this experiment, two receivers (COMNAV K500, Shanghai Siyue Technology Co. Ltd, Shanghai, China) are used, and both are connected to the GNSS antennae. The receiver electro-circuit is shown in Figure 10. The performance of the algorithm in Figure 9 is better than the algorithm in Figure 8. The algorithm is effective, and the robustness of the ambiguity function ( ( , , ) F x y z ) is improved. In the second aspect, the success rate of processing and not-processing is compared, as shown in Table 1.

Experimental Attitude Determination Results
In this section the RMM method is demonstrated using data collected in an experiment. The experiment is conducted on the top of a building. In this experiment, two receivers (COMNAV K500, Shanghai Siyue Technology Co. Ltd, Shanghai, China) are used, and both are connected to the GNSS antennae. The receiver electro-circuit is shown in Figure 10.  The baseline is constrained with 2.0 m and 0.5 m and the pitch angle search region of 20 degrees is provided by MEMS. During the initial step, the integer ambiguity combination is resolved by RMM method and the real value is work out after a few epochs, as is shown in Figure 11, and the real value is a point on the peak obviously. The baseline is constrained with 2.0 m and 0.5 m and the pitch angle search region of 20 degrees is provided by MEMS. During the initial step, the integer ambiguity combination is resolved by RMM method and the real value is work out after a few epochs, as is shown in Figure 11, and the real value is a point on the peak obviously. In this test, the number of locked GPS satellites is about eight, and the geometry of observations for this test is reasonably good. The experiment is performed from two groups in different arrangement. In Figure 12, with 2.0 m baseline, the standard deviation of the yaw and pitch angles are about 0.14° (1σ)and 0.18° (1σ), while with 0.5 m baseline the standard deviation reaches about 0.2° (1σ) and 0.25° (1σ) as shown in Figure 13. However, the calculation is reduced greatly.   In this test, the number of locked GPS satellites is about eight, and the geometry of observations for this test is reasonably good. The experiment is performed from two groups in different arrangement. In Figure 12, with 2.0 m baseline, the standard deviation of the yaw and pitch angles are about 0.14( 1σ) and 0.18˝(1σ), while with 0.5 m baseline the standard deviation reaches about 0.2˝(1σ) and 0.25( 1σ) as shown in Figure 13. However, the calculation is reduced greatly. The baseline is constrained with 2.0 m and 0.5 m and the pitch angle search region of 20 degrees is provided by MEMS. During the initial step, the integer ambiguity combination is resolved by RMM method and the real value is work out after a few epochs, as is shown in Figure 11, and the real value is a point on the peak obviously. In this test, the number of locked GPS satellites is about eight, and the geometry of observations for this test is reasonably good. The experiment is performed from two groups in different arrangement. In Figure 12, with 2.0 m baseline, the standard deviation of the yaw and pitch angles are about 0.14° (1σ)and 0.18° (1σ), while with 0.5 m baseline the standard deviation reaches about 0.2° (1σ) and 0.25° (1σ) as shown in Figure 13. However, the calculation is reduced greatly. The performance is good with high accuracy shown in Table 2.

Conclusions
This study describes the rotation matrix method and the relationship of the geometric model and the noise error for single-frequency and single-epoch GNSS attitude determination. The rotation matrix method reduces the computational complexity compared with the analytical resolution method [5,6]. In RMM, the calculation is reduced greatly and the error angles of different scales in trigonometric functions are effectively avoided. In contrast to the AFM, the geometric correlation methods [10][11][12][13][14][15] get better performance in solving the relationship of the geometric model and the noise error. Although the AFM is more flexible, there is a lack of analysis on this aspect. In the study, this relationship about the AFM is analyzed in detail. The computation error and the noise error are effectively treated. Not only is the flexibility of the AFM inherited, but the success rate is also increased. According to our simulations and real-time experiments, this method is verified as very reliable and effective. The computational complexity is greatly reduced and the success rate is effectively increased. In future studies, we plan to combine the method with an inertial navigation system (INS) for tight combination. The performance is good with high accuracy shown in Table 2.

Conclusions
This study describes the rotation matrix method and the relationship of the geometric model and the noise error for single-frequency and single-epoch GNSS attitude determination. The rotation matrix method reduces the computational complexity compared with the analytical resolution method [5,6]. In RMM, the calculation is reduced greatly and the error angles of different scales in trigonometric functions are effectively avoided. In contrast to the AFM, the geometric correlation methods [10][11][12][13][14][15] get better performance in solving the relationship of the geometric model and the noise error. Although the AFM is more flexible, there is a lack of analysis on this aspect. In the study, this relationship about the AFM is analyzed in detail. The computation error and the noise error are effectively treated. Not only is the flexibility of the AFM inherited, but the success rate is also increased. According to our simulations and real-time experiments, this method is verified as very reliable and effective. The computational complexity is greatly reduced and the success rate is effectively increased. In future studies, we plan to combine the method with an inertial navigation system (INS) for tight combination.