A Low-Cost INS-Integratable GNSS Ultra-Short Baseline Attitude Determination System

Traditional attitude determination using global navigation satellite system (GNSS) carrier phases is mostly applied on geodetic-grade receivers with sufficient baseline length. However, for some special applications such as mobile communication base station smart antenna attitude determination, only low-cost receivers with ultra-short baselines can be employed, and the environments are more challenging. When solving the ambiguity resolution (AR) problem with low-cost receivers, it is hard for the traditional methods in ambiguity domain to estimate float ambiguities accurately due to the large code pseudorange noises; thus, such systems fail to determine the correct ambiguities. Aiming at improving the AR success rate for ultra-short baselines attitude determination with low-cost receivers, we provide an objective function named Mean Square Residual (MSR) based on the geometrical relationship among the position spherical search space, the fractional carrier phases, and the possible ambiguities. The method can be calculated without code pseudoranges, and thus, can provide a higher AR success rate when using low-cost receivers. The corresponding analysis and acceptance test method are discussed in this contribution, and further, as an extension for more complicated urban dynamic applications, a GNSS/Inertial Navigation System (INS) integrated system is introduced. Several experiments have been conducted to verify performance.


Introduction
Attitude determination with differential global navigation satellite system (GNSS) carrier phases is performed in a wide range of applications, including surveying, vehicle active safety, and driver assistance system and aircraft/ship attitude determination.
Traditional GNSS-based attitude determination applications are mostly implemented with geodetic-grade receivers and antennas with sufficient length of baselines and are generally deployed in open areas that are less affected by multipaths. However, for some special scenarios, due to limits of budget and load capacity, only small and light-weight, low-cost receivers and antennas are applicable. For example, for mobile communication base station smart antenna attitude determination, because of the large application demand, manufacturers prefer the use of low-cost devices, and usually, the GNSS antennas are installed right next to the base station within a small area, bringing problems like signal blockage and strong multipaths. Besides that, some other applications like attitude determination for urban multi-rotor unmanned air vehicle (UAV) have extra requirement of instantaneity. As a common ground for the above-mentioned applications, because of the limited installation space, the antennas are usually placed at a short distance, forming an ultra-short baseline.

Measurement Models
The basic principle of attitude determination with GNSS, as shown in Figure 1, is to measure the precise relative position between multiple antennas fixed on the target surface, and then to determine the attitude of the target. If the baseline length is known, the baseline vector r rb in East-North-Up (ENU) coordinate whose origin point is the base antenna can be expressed as a function of the heading ϕ and pitch θ: where l is the known baseline length, e rb , n rb and u rb are the east, north and up components of r rb respectively.
Sensors 2018, 18, x 3 of 19 system is given, along with an extended system integrated with INS. Finally, the experimental results under different environments are presented to verify performance.

Measurement Models
The basic principle of attitude determination with GNSS, as shown in Figure 1, is to measure the precise relative position between multiple antennas fixed on the target surface, and then to determine the attitude of the target. If the baseline length is known, the baseline vector rb r in East-North-Up (ENU) coordinate whose origin point is the base antenna can be expressed as a function of the heading  and pitch  :  (1) where l is the known baseline length, rb e , rb n and rb u are the east, north and up components of rb r respectively. The undifferenced code pseudorange and carrier phase measurements can be modeled as [16]:  The undifferenced code pseudorange and carrier phase measurements can be modeled as [16]: P i r = ρ i r + cdt r − cdT i + T i r + I i r + E i r + MP i r,P + ε i r,P where P i r and Φ i r are the code pseudorange and carrier phase measurements between satellite i and receiver r, respectively, λ is the wavelength, ρ i r is the geometric range, c is the speed of light, dt r and dT i are the clock biases of the receiver and satellite, respectively, I i r is the ionospheric delay, T i r is the tropospheric delay, E i r is the ephemeris error, N i r is the unknown integer ambiguity, MP i r,P and MP i r,Φ are code and carrier multipath errors, respectively, and ε i r,P and ε i r,Φ represent all other errors that cannot be modeled.
To eliminate the effect of the clock bias of the satellite i, we can do single differencing between two antennas; for short-baseline cases, ephemeris and atmospheric errors can also be eliminated by single differencing, and we can get the single differenced (SD) measurements: ∆P i rb = P i r − P i b = ∆ρ i rb + c∆dt rb + ∆MP i rb,P + ∆ε i rb,P Further, to remove the effects of clock biases of receivers, double differencing between two satellites are done, forming the following double differenced (DD) measurements: For traditional AR methods in the ambiguity domain, the above DD code pseudorange and carrier phase measurement equations are usually combined to obtain the baseline vectorr rb and the float ambiguity estimation ∇∆ N. Then, for the LAMBDA method, the following objective function is used to select the integer ambiguity solution ∇∆ N based on least-squares principle [2].
where Q ∇∆N is the covariance matrix of the estimated float ambiguities. For the above ILS problem, the most common ambiguity acceptance test is the ratio-test [17,18], which is defined as follows Accept ∇∆ N iff : where ∇∆ N is integer ambiguity vector of the second smallest value of q(∇∆N), and thres q is the tolerance threshold. From (8), we may observe that the float ambiguity estimation ∇∆N plays an important role as a referenced center in the objective function q(∇∆N), and thus, determines the success rate of the AR problem. For low-cost receivers under challenging environments, due to the large pseudorange noises, float ambiguities ∇∆N cannot be estimated accurately, thus causing difficulties for the search and latter acceptance test judgment. For attitude determination applications with known baseline length, the objective function (8) can be modified by introducing additional term of tight or soft constraint [4,19], but the float ambiguity estimation still predominates the function.

Basic Idea of Ambiguity Resolution in Angle Domain
To get rid of the effect of float ambiguity estimation, in this contribution, we focus on the single-epoch AR algorithm in position domain, particularly, in angle domain. Furthermore, since the quality of code pseudoranges is hard to determine definitively for low-cost receivers, code pseudoranges are no longer used, which means only carrier phase measurements are employed in our algorithm.
To better understand the relationship among the baseline vector search space, the carrier phase measurements, and ambiguities, a formula variant is made (below). We can reserve only the fractional parts of the carrier phase measurements and rearrange the DD carrier phase measurement (7), as follows: 1 where ∇∆ Φ ij rb is the fractional part of the DD carrier phase measurement whose range is from −0.5 to 0.5 cycle, and ∇∆N ij rb is the DD ambiguity corresponding to ∇∆ Φ ij rb . The SD and DD geometric ranges can be approximated as: 1 where e i is the unit line-of-sight (LOS) vector pointing from the receiver to the satellite i. If we substitute (13) for (11), and ignore the multipath and other noises, we get: The right side of (14) is the baseline vector 1 λ r rb 's projection on e j − e i direction, and the value of the projection is given by the left side, which is composed of the fractional part of the DD carrier phase measurement and the unknown DD ambiguity. Figure 2 shows an error-free example for the relationship between the baseline vector (measured in unit of cycle) and the fractional parts of the DD carrier phase measurements in the ENU coordinate. The baseline length is 0.4 m, and the true pitch and heading are 0 • and 240 • , respectively. The blue, black, and red arrows represent the directions e 1 − e 2 , e 1 − e 3 and e 1 − e 4 , respectively, and the elevations and azimuths of the referenced satellite and other satellites are listed in Table 1. With knowledge of the baseline length, the possible baseline vectors fall into a spherical surface with a radius of the baseline length. For each differenced LOS vector e j − e i , according to (14), a set of parallel circles perpendicular to the vector can be drawn, each of which is determined by the fractional part of the DD carrier phase and possible integer DD ambiguity. If no errors exist, different sets of circles will intersect, as shown in Figure 2, and the intersection point is the true value of the baseline vector. However, affected by the subsistent noises, these circles may not intersect exactly at the same point. Therefore, our solution is to discretize the spherical surface into a set of candidate points by a 2-D traversal in angle space and select an approximate intersection point. The criterion for selecting the approximate intersection point in our algorithm is to see how far the point is from the closest circle of each DD measurement.  (13) where i e is the unit line-of-sight (LOS) vector pointing from the receiver to the satellite i . If we substitute (13) for (11), and ignore the multipath and other noises, we get: The right side of (14) is the baseline vector 1 rb  r 's projection on    j i ee direction, and the value of the projection is given by the left side, which is composed of the fractional part of the DD carrier phase measurement and the unknown DD ambiguity. Figure 2 shows an error-free example for the relationship between the baseline vector (measured in unit of cycle) and the fractional parts of the DD carrier phase measurements in the ENU coordinate.  (14), a set of parallel circles perpendicular to the vector can be drawn, each of which is determined by the fractional part of the DD carrier phase and possible integer DD ambiguity. If no errors exist, different sets of circles will intersect, as shown in Figure 2, and the intersection point is the true value of the baseline vector. However, affected by the subsistent noises, these circles may not intersect exactly at the same point. Therefore, our solution is to discretize the spherical surface into a set of candidate points by a 2-D traversal in angle space and select an approximate intersection point. The criterion for selecting the approximate intersection point in our algorithm is to see how far the point is from the closest circle of each DD measurement.  r , we can calculate its distance to the nearest circles for all the pairs of satellites. However, when calculating the "distances", considering the computational complexity, we employ the vertical distances to the circle planes instead of the direct distances to the circle edges. For each DD geometric range, we can process this in the following way:   Elevation By varying in angle space, we obtain a candidate position point set r rb,k on the spherical surface. For a candidate position point r rb,k , we can calculate its distance to the nearest circles for all the pairs of satellites. However, when calculating the "distances", considering the computational complexity, we employ the vertical distances to the circle planes instead of the direct distances to the circle edges. For each DD geometric range, we can process this in the following way: Again, we substitute (15) into 1 λ r rb,k ' s approximate expression in (13) form, we can get the baseline vector 1 λ r rb,k 's projection expression, with actual measurement noises as follows: The term containing the DD fractional carrier phase and the integer in (16) is the distance from the original point to the nearest circle plane, and then, the vertical distance to the nearest circle plane can be easily obtained by the remaining term: Since the actual carrier phase measurement ∇∆ Φ ij rb contains error, which is scaled by 1/ e j − e i in the above distance calculation formula, we can multiply the calculated distance by a weight e j − e i before making statistics for all the distances. Thus, we can define a so-called mean square residual (MSR) as the objection function. For the candidate position point r rb,k , its MSR value is: where · 2 is the 2-norm operator (mean square function It is reasonable to believe that the best candidate position point is the one with the smallest MSR, and once decided, the ambiguities can be decided by the nearest circles of the point: From the calculation formula, we can find the essence of MSR is to judge the fractional differences between the calculated geometric ranges and the carrier phase measurements. For low-cost receivers, carrier phases are less effected by multipath than code pseudoranges. Compared to traditional methods in ambiguity domain using code pseudoranges, our algorithm uses carrier phases only. The algorithm is more focused on the relationship between the position search space and the satellite geometric distribution, which leads to a higher success rate for AR problem; thus, the main advantage of the proposed method is the applicability to low-cost receivers.

Candidate Position Point Set
In this section, the derivation of the candidate position point set, or specifically, how to discretize the spherical search space more efficiently, is discussed. The aim is to minimize the number of the candidate positions under the precondition of the algorithm stability, so as to lighten the calculation burden.
To discretize the spherical search space more efficiently, we firstly vary θ at a certain search step to get a set of horizontal circles, and then divide each circle by varying ϕ at a search step length determined by the candidate θ of the circle. For attitude determination application with a constant baseline length, the candidate position point set only needs to be calculated once. According to the traversal sequence, the search step length ∆θ for θ needs to be decided before the decision of the search step length ∆ϕ for ϕ. As for the selection of a proper search step, to ensure the validity of the selected approximate intersection point, the searching steps must be small enough that the closest candidate position point to the true position can be selected as the best one, and can provide correct ambiguities. However, on the other hand, excessively small search step lengths may increase the burden of calculation. Therefore, to balance the reliability and efficiency, the selection of ∆θ and ∆ϕ is of great importance. With these considerations in mind, we will discuss the upper bounds of the search step lengths below.
In detail, for the sake of validity, we hope the largest difference for 1 λ ∇∆ρ ij rb caused by angle (θ and ϕ) difference between the true value and its nearest candidate is smaller than 0.25 cycle, which means in extreme situations, the largest difference for SD geometric range 1 λ ∆ρ i rb should be smaller than 0.125 cycle. If we convert both the unit LOS vector e i and the baseline vector r rb into angle forms, the SD geometric range 1 λ ∆ρ i rb can be presented as a function with baseline length and angle variables only. Firstly, to obtain the upper bound of ∆θ, we can take a derivative of 1 λ ∆ρ i rb with respect to θ, and rearrange the expression using trigonometric expressions: where el i and az i are the elevation and azimuth of satellite i, respectively. Since the biggest difference between the true θ and its nearest candidate value is half the search step length, i.e., ∆θ/2, we can get: For the right side of (22), the following inequality always holds, no matter what values the angles are.
Thus, we can set the universal upper bound of ∆θ, as below: For the range of ∆ϕ, similarly, we can take the derivative of 1 λ ∆ρ i rb with respect to ϕ as follows: and get the corresponding range of ∆ϕ as: Then, the upper bound of ∆ϕ can be set as: From (24) and (27), we can find that the minimum total number of the candidate position points has an approximate quadratic polynomial relationship with the baseline length; thus, the algorithm is more suitable for the ultra-short baseline situation. Furthermore, for the attitude determination system, if there are some other aiding sensors like INS or tilt sensors which can provide pitch information, the position search space can be further compressed, from a spherical surface to a circle determined by the pitch. In this way, we can not only improve the calculation efficiency of the algorithm, but also increase the success rate of the AR method, especially under the harsh environment with few available satellites and strong multipath interference.

Algorithm Performance Analysis and Acceptance Test Method
We have introduced the principle of the proposed AR algorithm in angle domain. In this part, the algorithm performance and the corresponding acceptance test method will be discussed.
The performance of the algorithm, or whether the AR problem can be solved successfully, depends on whether the selection of the approximate intersection point is correct. Figure 3 gives a MSR distribution example in natural logarithmic form. According to the best candidate point selection criterion, we can get the approximate intersection point at the global minimum point. As we can see from the figure, in addition to the global minimum point, there are other valleys, namely local minimums. The MSR distribution can be affected by several factors such as satellite geometry distribution, baseline length and carrier phase noises; in some cases, the global minimum may appear at a wrong point.
From (24) and (27), we can find that the minimum total number of the candidate position points has an approximate quadratic polynomial relationship with the baseline length; thus, the algorithm is more suitable for the ultra-short baseline situation. Furthermore, for the attitude determination system, if there are some other aiding sensors like INS or tilt sensors which can provide pitch information, the position search space can be further compressed, from a spherical surface to a circle determined by the pitch. In this way, we can not only improve the calculation efficiency of the algorithm, but also increase the success rate of the AR method, especially under the harsh environment with few available satellites and strong multipath interference.

Algorithm Performance Analysis and Acceptance Test Method
We have introduced the principle of the proposed AR algorithm in angle domain. In this part, the algorithm performance and the corresponding acceptance test method will be discussed.
The performance of the algorithm, or whether the AR problem can be solved successfully, depends on whether the selection of the approximate intersection point is correct. Figure 3 gives a MSR distribution example in natural logarithmic form. According to the best candidate point selection criterion, we can get the approximate intersection point at the global minimum point. As we can see from the figure, in addition to the global minimum point, there are other valleys, namely local minimums. The MSR distribution can be affected by several factors such as satellite geometry distribution, baseline length and carrier phase noises; in some cases, the global minimum may appear at a wrong point. Even for error-free situations, under some special satellite geometry distributions, the phenomenon of multiple intersection points will occur. An example of multiple intersection points is given in Figure 4., where there are three pairs of DD satellites. The baseline length is 0.25 m, and the true pitch and heading are 0° and 90°, respectively. The blue, black and red arrows represent the Even for error-free situations, under some special satellite geometry distributions, the phenomenon of multiple intersection points will occur. An example of multiple intersection points is given in Figure 4, where there are three pairs of DD satellites. The baseline length is 0.25 m, and the true pitch and heading are 0 • and 90 • , respectively. The blue, black and red arrows represent the directions of e 1 − e 2 , e 1 − e 3 and e 1 − e 4 , respectively, and the elevations and azimuths of the referenced satellite and other satellites are listed in Table 2. The true position point is marked by green, but from the figure we can find that in addition to the true position point, another two intersection points exist. Actually, such a multi-intersection case requires a special integer linear dependence property of the differenced LOS vector set e j − e i i , which happens mostly under the situation where only 4 or 5 satellites are available. referenced satellite and other satellites are listed in Table 2. The true position point is marked by green, but from the figure we can find that in addition to the true position point, another two intersection points exist. Actually, such a multi-intersection case requires a special integer linear   Table 2. Satellite geometry distribution for Figure 4.

Sat 1 (Ref) Sat 2 Sat 3 Sat 4 Elevation
To illustrate the effects of different factors on the MSR method, a success rate simulation is conducted, as shown in Figure 5b, with different baseline lengths l , 0.2 m, 0.4 m and 0.8 m respectively, and different carrier phase noise standard deviations (STD) ranging from 0 to 0.1 cycle. For each simulation set, the number of the test samples is 5,000, and the success rate is defined as the percentage of the results that get the correct integer ambiguities. The satellite distribution is plotted in Figure 5a. As we can see from the figure, the success rate drops as the carrier phase noises rise, and the algorithm achieves better performance with short baselines. Also, we conduct a similar simulation under a poor satellite geometry distribution with only five satellites in view, all of which are in the northern sky. The success rate and the satellite geometry distribution are plotted in Figure 6. With such satellite geometry distribution, the success rate drops greatly in comparison with the situation under normal geometry distribution.  Elevation To illustrate the effects of different factors on the MSR method, a success rate simulation is conducted, as shown in Figure 5b, with different baseline lengths l, 0.2 m, 0.4 m and 0.8 m respectively, and different carrier phase noise standard deviations (STD) ranging from 0 to 0.1 cycle. For each simulation set, the number of the test samples is 5,000, and the success rate is defined as the percentage of the results that get the correct integer ambiguities. The satellite distribution is plotted in Figure 5a. As we can see from the figure, the success rate drops as the carrier phase noises rise, and the algorithm achieves better performance with short baselines.  Table 2. The true position point is marked by green, but from the figure we can find that in addition to the true position point, another two intersection points exist. Actually, such a multi-intersection case requires a special integer linear   Table 2. Satellite geometry distribution for Figure 4.

Sat 1 (Ref) Sat 2 Sat 3 Sat 4 Elevation
To illustrate the effects of different factors on the MSR method, a success rate simulation is conducted, as shown in Figure 5b, with different baseline lengths l , 0.2 m, 0.4 m and 0.8 m respectively, and different carrier phase noise standard deviations (STD) ranging from 0 to 0.1 cycle. For each simulation set, the number of the test samples is 5,000, and the success rate is defined as the percentage of the results that get the correct integer ambiguities. The satellite distribution is plotted in Figure 5a. As we can see from the figure, the success rate drops as the carrier phase noises rise, and the algorithm achieves better performance with short baselines. Also, we conduct a similar simulation under a poor satellite geometry distribution with only five satellites in view, all of which are in the northern sky. The success rate and the satellite geometry distribution are plotted in Figure 6. With such satellite geometry distribution, the success rate drops greatly in comparison with the situation under normal geometry distribution. Normally, the shorter the baseline, the higher success rate; but sometimes this is not the case (for example, in Figure 6b in two similar directions are less distinguishable for shorter baselines, and thus, increase the difficulty in finding the correct approximate intersection point. Therefore, with challenges such as poor satellite geometry distribution and large carrier phase noises, the MSR distinction among the local minimums tends to be less obvious, and it is hard to find the correct approximate intersection point. To judge whether the result is valid or not, or strictly speaking, to assess the reliability of the selected approximate intersection point, we propose an acceptance test method based on MSR distinction. It is reasonable to believe that the more obvious the distinction between the smallest valley and others, the higher the probability of finding the correct intersection point.

Flow Chart of the MSR-Based GNSS/INS Integrated Attitude Determination System
As analyzed in the previous section, the success rate of the MSR method with GNSS is vulnerable to the poor satellite geometry and multipath effects. For attitude determination applications in urban area, challenging environments are inevitable, which makes the availability, precision, and instantaneity of the attitude results hard to guarantee. To settle the problem, our solution is to integrate the MSR-based GNSS attitude determination system with INS.
INS can provide higher precision attitude results at a higher output rate and keep on working for a short while when the GNSS is unavailable. What's more, INS can output pitch information with Normally, the shorter the baseline, the higher success rate; but sometimes this is not the case (for example, in Figure 6b, the success rate of l = 0.2 m is lower than the one of l = 0.4 m), because under extreme poor satellite geometry with few visible satellites, the circles decided by e j − e i i in two similar directions are less distinguishable for shorter baselines, and thus, increase the difficulty in finding the correct approximate intersection point.
Therefore, with challenges such as poor satellite geometry distribution and large carrier phase noises, the MSR distinction among the local minimums tends to be less obvious, and it is hard to find the correct approximate intersection point. To judge whether the result is valid or not, or strictly speaking, to assess the reliability of the selected approximate intersection point, we propose an acceptance test method based on MSR distinction. It is reasonable to believe that the more obvious the distinction between the smallest valley and others, the higher the probability of finding the correct intersection point. In detail, if MSR 1st and MSR 2st represent the MSR values of the smallest and second-smallest valleys: Accept r rb,1st iff : where r rb,1st is the precise baseline vector re-calculated by the integer ambiguities of the smallest valley, and thres MSR is the acceptance test threshold. One important thing to note is that the above MSR ratio is the ratio of the second smallest valley to the smallest valley, not that of the second smallest MSR to the smallest MSR. Sometimes, the points with smallest and second smallest MSRs are neighbors of the true position, which means that the distinction between them may be very small, and both can give the same correct integer ambiguities and be regarded as the correct approximate intersection points. For the selection of the threshold, thres MSR = 1.3 is recommended, with which the judgment success rate can reach 90% in most ultra-short baseline attitude determination tests we have conducted (an example is given by the multipath experiment in Section 7.1.2).

Flow Chart of the MSR-Based GNSS/INS Integrated Attitude Determination System
As analyzed in the previous section, the success rate of the MSR method with GNSS is vulnerable to the poor satellite geometry and multipath effects. For attitude determination applications in urban area, challenging environments are inevitable, which makes the availability, precision, and instantaneity of the attitude results hard to guarantee. To settle the problem, our solution is to integrate the MSR-based GNSS attitude determination system with INS.
INS can provide higher precision attitude results at a higher output rate and keep on working for a short while when the GNSS is unavailable. What's more, INS can output pitch information with which only the traversal in heading is needed for the derivation of the candidate position point set.
For the accuracy of the pitch provided by INS, according to the derivation process of ∆θ, the absolute error must be less than 0.125λ/l.
A flow chart of the attitude determination system is depicted in Figure 7. The proposed system contains two main processing modules, GNSS processing module and INS processing module, and supports two attitude determination modes, GNSS only mode and GNSS/INS integration mode. which only the traversal in heading is needed for the derivation of the candidate position point set.
For the accuracy of the pitch provided by INS, according to the derivation process of   , the absolute error must be less than  0.125 / l . A flow chart of the attitude determination system is depicted in Figure 7. The proposed system contains two main processing modules, GNSS processing module and INS processing module, and supports two attitude determination modes, GNSS only mode and GNSS/INS integration mode. The main function of the GNSS processing module is to provide various GNSS results, including the positioning, velocity and attitude results. Every epoch, before calculation, satellite selection is made based on satellite elevation mask angle. Then, according to the method described previously, the MSR values are calculated for all the candidate position points which are derived with the baseline length and the pitch if available. Next, we pick up the candidate position points of the smallest and second-smallest MSR valleys, then calculate the MSR ratio between them. To judge whether the GNSS attitude result is valid, the following judging condition (JC), which includes the acceptance test method of (28) and the tests of the prior knowledge, is used: thres are the thresholds for baseline length and pitch, respectively.
In our system, the traditional loose coupling method [20] is adopted to calibrate INS with the results obtained from the GNSS processing module. To enhance the reliability of the INS output, the validity of the GNSS results needs to be judged. Although the acceptance test method based on the MSR ratio has a low missing detection rate, the false alarm rate is relatively high. Sometimes, with a poor satellite geometry distribution and strong multipath effects, the candidate position point of the second-smallest valley, not the smallest one, is the correct, needed intersection point. Thus, to get The main function of the GNSS processing module is to provide various GNSS results, including the positioning, velocity and attitude results. Every epoch, before calculation, satellite selection is made based on satellite elevation mask angle. Then, according to the method described previously, the MSR values are calculated for all the candidate position points which are derived with the baseline length and the pitch if available. Next, we pick up the candidate position points of the smallest and second-smallest MSR valleys, then calculate the MSR ratio between them. To judge whether the GNSS attitude result is valid, the following judging condition (JC), which includes the acceptance test method of (28) and the tests of the prior knowledge, is used: where θ rb,1st is the pitch calculated by r rb,1st , θ is the pitch provided by INS, and, thres l and thres θ are the thresholds for baseline length and pitch, respectively. In our system, the traditional loose coupling method [20] is adopted to calibrate INS with the results obtained from the GNSS processing module. To enhance the reliability of the INS output, the validity of the GNSS results needs to be judged. Although the acceptance test method based on the MSR ratio has a low missing detection rate, the false alarm rate is relatively high. Sometimes, with a poor satellite geometry distribution and strong multipath effects, the candidate position point of the second-smallest valley, not the smallest one, is the correct, needed intersection point. Thus, to get more valid GNSS attitude results, we leave two candidate position points, namely, those of the smallest and second-smallest valley, for selection, and calculate their precise headings. Then, we adopt the Innovation Pre-filtering based Extended Kalman Filter (IPEKF) [21] to judge the validity of these two candidate headings, along with the other GNSS results which will be used in the subsequent filtering process.
The state and observation equation are given as follows [20]: represents the INS state of epoch k, composed of the position vector r I,k , velocity vector v I,k , attitude vectors α k , accelerometer biases x a,k and gyroscope biases x g,k , y k = r T G,k v T G,k ϕ k T represents measurements of epoch k provided by the GNSS processing module, composed of the position vector r G,k , velocity vector v G,k and heading ϕ k , ω k−1 and ν k represent system and measurement noises, respectively, F k−1 represents the transition matrix whose expression can be found in Section 11.6 of [20], and H k represents the measurement matrix.
The basic idea of IPEKF-based judgment is to test the consistency between the measurements and the state estimators with the use of the following measurement innovation, which is given by: wherex − k is the EKF predicted state estimator derived fromx k−1 . For the two candidate headings, we will select the one with a smaller innovation value for the subsequent validity judgment and EKF process. Before we judge the validity of each measurement innovation component with a certain threshold, the following normalization needs to be done: where δy − k,j is the jth element of δy − k and c − k,j,j is the jth diagonal element of the innovation covariance matrix C − k which is calculated by: where P − k and R k is the covariance matrixes of the predicted state estimators and the measurement noises in EKF process, respectively. If the normalized innovations exceed the thresholds, the corresponding measurements will be judged invalid, and the corresponding state components won't be filtered at this epoch.
Note, for the initialization stage without reliable INS predicted information, JC of (29) is adopted to judge the validity of the GNSS results. Once calibrated, INS can keep providing valid results for a short time, even losing the aid of the GNSS; however, errors will increase as time goes by. Therefore, in our designed system, if the time of losing valid GNSS results exceeds 30 s, INS will not output heading results any longer; after that, once there are new GNSS results (JC passed) available, the initialization process needs to be re-performed.

Experiment Results
Several static/dynamic experiments in GNSS only or GNSS/INS integration mode under different environments were conducted to evaluate the practical performance of the proposed method. In our system, to realize GNSS attitude determination, two Unicore ® UM220-III H GNSS modules were employed, along with two SenseStone ® GNSS built-in active antennas with double feed points. For INS part, Xsens ® MTi-3-8A7G6 INS module was employed. The raw data can be transmitted in RTCM 3.2 format by the GNSS modules. The message types of the measurements and ephemeris we used are listed in Table 3. In the experiments, the satellite elevation masks were set as 15 • . To better verify the correctness of the results and make an analysis on the measurement performance of the low-cost antennas and receivers under different environments, reference baseline values (for static experiments only) was obtained based on the long-term observation results from the TRIMBLE ® NETR9 multi-frequency geodetic-grade receivers in advance. For short baselines, the measurements in DD form are considered to be free of the ephemeris, atmospheric, and clock errors. Hence, if we calculate the DD geometric ranges based on the precise reference baseline and subtract them from the DD measurements, the obtained residuals can be regarded as noises from multipath and receivers only. With the DD residuals, we can evaluate the measurement qualities of the patch antennas and low-cost receivers that we used under different environments.

Experiments with GNSS Only
In this part, the performance of the proposed system with GNSS only was tested.

Experiments with GNSS Only in Open Area
Firstly, the static experiment was performed in an open area. The antennas were placed on the arm of a turntable, as shown in Figure 8a, with a baseline length of 0.267 m and a certain pitch. The DD residuals of the measurements were calculated to assess the code pseudorange and carrier phase quality of the low-cost antennas and receivers under the open area, which are given in Figure 8b. The number of the available GPS satellites was 6, and by statistical analysis, the standard deviations of DD code pseudorange noises range from 0.70 m to 1.07 m and the ones of DD carrier phase noises were less than 0.013 cycle.  Table 3. In the experiments, the satellite elevation masks were set as 15°. To better verify the correctness of the results and make an analysis on the measurement performance of the low-cost antennas and receivers under different environments, reference baseline values (for static experiments only) was obtained based on the long-term observation results from the TRIMBLE ® NETR9 multi-frequency geodetic-grade receivers in advance. For short baselines, the measurements in DD form are considered to be free of the ephemeris, atmospheric, and clock errors. Hence, if we calculate the DD geometric ranges based on the precise reference baseline and subtract them from the DD measurements, the obtained residuals can be regarded as noises from multipath and receivers only. With the DD residuals, we can evaluate the measurement qualities of the patch antennas and low-cost receivers that we used under different environments.

Experiments with GNSS Only
In this part, the performance of the proposed system with GNSS only was tested.

Experiments with GNSS Only in Open Area
Firstly, the static experiment was performed in an open area. The antennas were placed on the arm of a turntable, as shown in Figure 8a, with a baseline length of 0.267 m and a certain pitch. The DD residuals of the measurements were calculated to assess the code pseudorange and carrier phase quality of the low-cost antennas and receivers under the open area, which are given in Figure 8b. The number of the available GPS satellites was 6, and by statistical analysis, the standard deviations of DD code pseudorange noises range from 0.70 m to 1.07 m and the ones of DD carrier phase noises were less than 0.013 cycle. The heading and pitch results of the proposed MSR based algorithm are plotted in Figure 9, in which the green/yellow points present the results that pass/fail the acceptance test method (AT) of (28). The red dashed lines in Figure 9 are the averages results of TRIMBLE ® NETR9 multi-frequency The heading and pitch results of the proposed MSR based algorithm are plotted in Figure 9, in which the green/yellow points present the results that pass/fail the acceptance test method (AT) of (28). The red dashed lines in Figure 9 are the averages results of TRIMBLE ® NETR9 multi-frequency geodetic-grade receivers, which are used as reference value here. As we can see from Figure 8, the pseudorange and carrier phase measurement noises of the low-cost receivers in open area are also quite small. The success rate of the MSR method was 100%, and the standard deviations of the headings and pitches are 0.42 • and 0.68 • , respectively. geodetic-grade receivers, which are used as reference value here. As we can see from Figure 8, the pseudorange and carrier phase measurement noises of the low-cost receivers in open area are also quite small. The success rate of the MSR method was 100%, and the standard deviations of the headings and pitches are 0.42°and 0.68°, respectively. Then, a dynamic experiment was conducted at the same place. This time, the turntable was moving at a rotational speed of about 36°/s. The attitude results are plotted in Figure 10. According to the angular rate converted from the heading results, the obtained attitude results were consistent with what we set.

Experiment with GNSS Only under Multipath Environment
As one of the applications for which we want to adopt the proposed system, the situation of the attitude determination application for the base station smart antennas was simulated. We conducted a static experiment under multipath environment and placed the GNSS antennas next to a pylon, as shown in Figure 11a. The DD residuals of the measurements under such situations were given in Figure 11b. As we can see from the figure, under such an environment, both code pseudorange and carrier phase measurement qualities deteriorated dramatically, but the carrier phase accuracy was still within the tolerable range. The standard deviations of the DD pseudorange measurements ranged from 3.8 m to 6.7 m, while the ones of the DD carrier phase measurements ranged from 0.07 cycle to 0.12 cycle. Then, a dynamic experiment was conducted at the same place. This time, the turntable was moving at a rotational speed of about 36 • /s. The attitude results are plotted in Figure 10. According to the angular rate converted from the heading results, the obtained attitude results were consistent with what we set. geodetic-grade receivers, which are used as reference value here. As we can see from Figure 8, the pseudorange and carrier phase measurement noises of the low-cost receivers in open area are also quite small. The success rate of the MSR method was 100%, and the standard deviations of the headings and pitches are 0.42°and 0.68°, respectively. Then, a dynamic experiment was conducted at the same place. This time, the turntable was moving at a rotational speed of about 36°/s. The attitude results are plotted in Figure 10. According to the angular rate converted from the heading results, the obtained attitude results were consistent with what we set.

Experiment with GNSS Only under Multipath Environment
As one of the applications for which we want to adopt the proposed system, the situation of the attitude determination application for the base station smart antennas was simulated. We conducted a static experiment under multipath environment and placed the GNSS antennas next to a pylon, as shown in Figure 11a. The DD residuals of the measurements under such situations were given in Figure 11b. As we can see from the figure, under such an environment, both code pseudorange and carrier phase measurement qualities deteriorated dramatically, but the carrier phase accuracy was still within the tolerable range. The standard deviations of the DD pseudorange measurements ranged from 3.8 m to 6.7 m, while the ones of the DD carrier phase measurements ranged from 0.07 cycle to 0.12 cycle.

Experiment with GNSS Only under Multipath Environment
As one of the applications for which we want to adopt the proposed system, the situation of the attitude determination application for the base station smart antennas was simulated. We conducted a static experiment under multipath environment and placed the GNSS antennas next to a pylon, as shown in Figure 11a. The DD residuals of the measurements under such situations were given in Figure 11b. As we can see from the figure, under such an environment, both code pseudorange and carrier phase measurement qualities deteriorated dramatically, but the carrier phase accuracy was still In addition to the MSR method, the LAMBDA-based methods were also tested this time. The attitude results of the MSR method and the LAMBDA method that uses EKF to estimate float ambiguities [22] are plotted in Figure 12a. For the LAMBDA method, the traditional ratio-test method (10) was adopted as the acceptance test method with a threshold of 2 [17]. As we can see from the figure, affected by the large pseudorange noises, the float ambiguity estimator could hardly converge, even with the use of EKF, which made the traditional LAMBDA method unable to fix the integer ambiguities. In contrast, our single-epoch MSR method relies on carrier phase only, and thus, could still provide correct results for most of the test time, despite of the attitude precision decline due to the multipath noises and the limit of the ultra-short baseline.
Also, we made a comparison with the results of the single-epoch LAMBDA method, along with single-epoch LAMBDA methods with tight/soft constraints (TC/SC-LAMBDA) [4,19]. The results are listed in Table 4, including the total success rates of all results, the fix rates based on the corresponding test methods, and the success rates of the fix results. As shown by the table, compared with the LAMBDA-based methods, our MSR method could provide a much higher single-epoch success rate under the multipath environment, and the judgment of the corresponding acceptance test method was also more reliable. In addition to the MSR method, the LAMBDA-based methods were also tested this time. The attitude results of the MSR method and the LAMBDA method that uses EKF to estimate float ambiguities [22] are plotted in Figure 12a. For the LAMBDA method, the traditional ratio-test method (10) was adopted as the acceptance test method with a threshold of 2 [17]. In addition to the MSR method, the LAMBDA-based methods were also tested this time. The attitude results of the MSR method and the LAMBDA method that uses EKF to estimate float ambiguities [22] are plotted in Figure 12a. For the LAMBDA method, the traditional ratio-test method (10) was adopted as the acceptance test method with a threshold of 2 [17]. As we can see from the figure, affected by the large pseudorange noises, the float ambiguity estimator could hardly converge, even with the use of EKF, which made the traditional LAMBDA method unable to fix the integer ambiguities. In contrast, our single-epoch MSR method relies on carrier phase only, and thus, could still provide correct results for most of the test time, despite of the attitude precision decline due to the multipath noises and the limit of the ultra-short baseline.
Also, we made a comparison with the results of the single-epoch LAMBDA method, along with single-epoch LAMBDA methods with tight/soft constraints (TC/SC-LAMBDA) [4,19]. The results are listed in Table 4, including the total success rates of all results, the fix rates based on the corresponding test methods, and the success rates of the fix results. As shown by the table, compared with the LAMBDA-based methods, our MSR method could provide a much higher single-epoch success rate under the multipath environment, and the judgment of the corresponding acceptance test method was also more reliable. As we can see from the figure, affected by the large pseudorange noises, the float ambiguity estimator could hardly converge, even with the use of EKF, which made the traditional LAMBDA method unable to fix the integer ambiguities. In contrast, our single-epoch MSR method relies on carrier phase only, and thus, could still provide correct results for most of the test time, despite of the attitude precision decline due to the multipath noises and the limit of the ultra-short baseline.
Also, we made a comparison with the results of the single-epoch LAMBDA method, along with single-epoch LAMBDA methods with tight/soft constraints (TC/SC-LAMBDA) [4,19]. The results are listed in Table 4, including the total success rates of all results, the fix rates based on the corresponding test methods, and the success rates of the fix results. As shown by the table, compared with the LAMBDA-based methods, our MSR method could provide a much higher single-epoch success rate under the multipath environment, and the judgment of the corresponding acceptance test method was also more reliable.

Experiment with GNSS/INS Integration
The GNSS-only mode is more suitable to the applications that have low requirement of instantaneity. However, considering the challenges like multipath effects and satellite availability, the requirement of instantaneous output is hard to guarantee with GNSS only for the dynamic applications in urban environments; thus, the GNSS/INS integration mode is more recommended.
A dynamic experiment with GNSS/INS integration was conducted in an actual urban environment with surrounding trees and buildings. The moving route is given by Figure 13. To better evaluate the performance of the proposed system, the attitude results of NovAtel ® SPAN-CPT high precision GNSS/INS receiver were also collected at the same time. The available satellite numbers and attitude results of SPAN-CPT and the proposed system are plotted in Figures 14 and 15, respectively.

Experiment with GNSS/INS Integration
The GNSS-only mode is more suitable to the applications that have low requirement of instantaneity. However, considering the challenges like multipath effects and satellite availability, the requirement of instantaneous output is hard to guarantee with GNSS only for the dynamic applications in urban environments; thus, the GNSS/INS integration mode is more recommended.
A dynamic experiment with GNSS/INS integration was conducted in an actual urban environment with surrounding trees and buildings. The moving route is given by Figure 13. To better evaluate the performance of the proposed system, the attitude results of NovAtel ® SPAN-CPT high precision GNSS/INS receiver were also collected at the same time. The available satellite numbers and attitude results of SPAN-CPT and the proposed system are plotted in Figures 14 and  15, respectively.

Experiment with GNSS/INS Integration
The GNSS-only mode is more suitable to the applications that have low requirement of instantaneity. However, considering the challenges like multipath effects and satellite availability, the requirement of instantaneous output is hard to guarantee with GNSS only for the dynamic applications in urban environments; thus, the GNSS/INS integration mode is more recommended.
A dynamic experiment with GNSS/INS integration was conducted in an actual urban environment with surrounding trees and buildings. The moving route is given by Figure 13. To better evaluate the performance of the proposed system, the attitude results of NovAtel ® SPAN-CPT high precision GNSS/INS receiver were also collected at the same time. The available satellite numbers and attitude results of SPAN-CPT and the proposed system are plotted in Figures 14 and  15, respectively.    As we can see from the figures, for the dynamic situation in urban environment, the availabilities of both systems declined to a certain extent. Still, the advantages of the proposed system can be observed. For comparison purposes, three route segments have been marked in the figures, and we can find some obvious mistakes in the attitude results of SPAN-CPT, including a two-minute sustained wrong output in Route Segment C. In contrast, our proposed system provided a more stable output with fewer anomalies, even though we realized the system with the use of the low-cost antennas and receivers, and had fewer available satellites for calculation. From Route Segment B, we can find that the INS part was re-initialized quickly and accurately with the MSR-based GNSS attitude determination method, which has a higher success rate for single-epoch solutions.
The introduction of the IPEKF further enhances the robustness of the whole system, and reduces the adverse effects of the instability of the GNSS output. Figure 16 gives the IPEKF attitude selection results under different satellite availability conditions during the test. The total number of the epochs that can provide more than three available satellites is 601. For the situation with sufficient satellites, the attitude results with the smallest MSR value still had high reliability, but for the situation with less satellites, sometimes the attitude results of the second-smallest MSR valley were correct. With the aid of the INS and IPEKF, it is much easier to select the correct results or give a validity judgment for the results, and thus, enhance the robustness of the whole system.   As we can see from the figures, for the dynamic situation in urban environment, the availabilities of both systems declined to a certain extent. Still, the advantages of the proposed system can be observed. For comparison purposes, three route segments have been marked in the figures, and we can find some obvious mistakes in the attitude results of SPAN-CPT, including a two-minute sustained wrong output in Route Segment C. In contrast, our proposed system provided a more stable output with fewer anomalies, even though we realized the system with the use of the low-cost antennas and receivers, and had fewer available satellites for calculation. From Route Segment B, we can find that the INS part was re-initialized quickly and accurately with the MSR-based GNSS attitude determination method, which has a higher success rate for single-epoch solutions.
The introduction of the IPEKF further enhances the robustness of the whole system, and reduces the adverse effects of the instability of the GNSS output. Figure 16 gives the IPEKF attitude selection results under different satellite availability conditions during the test. The total number of the epochs that can provide more than three available satellites is 601. For the situation with sufficient satellites, the attitude results with the smallest MSR value still had high reliability, but for the situation with less satellites, sometimes the attitude results of the second-smallest MSR valley were correct. With the aid of the INS and IPEKF, it is much easier to select the correct results or give a validity judgment for the results, and thus, enhance the robustness of the whole system. As we can see from the figures, for the dynamic situation in urban environment, the availabilities of both systems declined to a certain extent. Still, the advantages of the proposed system can be observed. For comparison purposes, three route segments have been marked in the figures, and we can find some obvious mistakes in the attitude results of SPAN-CPT, including a two-minute sustained wrong output in Route Segment C. In contrast, our proposed system provided a more stable output with fewer anomalies, even though we realized the system with the use of the low-cost antennas and receivers, and had fewer available satellites for calculation. From Route Segment B, we can find that the INS part was re-initialized quickly and accurately with the MSR-based GNSS attitude determination method, which has a higher success rate for single-epoch solutions.
The introduction of the IPEKF further enhances the robustness of the whole system, and reduces the adverse effects of the instability of the GNSS output. Figure 16 gives the IPEKF attitude selection results under different satellite availability conditions during the test. The total number of the epochs that can provide more than three available satellites is 601. For the situation with sufficient satellites, the attitude results with the smallest MSR value still had high reliability, but for the situation with less satellites, sometimes the attitude results of the second-smallest MSR valley were correct. With the aid of the INS and IPEKF, it is much easier to select the correct results or give a validity judgment for the results, and thus, enhance the robustness of the whole system.

Conclusions
To settle the AR problem for ultra-short baseline attitude determination with low-cost GNSS receivers under challenging environment, we proposed a single-epoch AR algorithm, namely MSR, in the angle domain. Unlike traditional algorithms in the ambiguity domain, MSR method neither estimates the float ambiguities nor uses the code pseudoranges, relying rather on carrier phases only which are less affected by multipath; thus, it has a certain adaptation to challenging environments. The performance of the MSR method and the extension MSR-based GNSS/INS integrated attitude determination system for urban dynamic applications under different environments with the low-cost devices has been validated. The experiment results show that the proposed systems could provide a higher success rate and reliability when using low-cost receivers even under the challenging environments.
The above discussion and experiments only involves two degrees of freedom of attitude with one baseline, but the application can be easily extended to three-degree-of-freedom situation with two baselines, which can be used for UAV attitude determination, and will be tried in our future works.