Real-Time Single-Frequency GPS/MEMS-IMU Attitude Determination of Lightweight UAVs

In this paper, a newly-developed direct georeferencing system for the guidance, navigation and control of lightweight unmanned aerial vehicles (UAVs), having a weight limit of 5 kg and a size limit of 1.5 m, and for UAV-based surveying and remote sensing applications is presented. The system is intended to provide highly accurate positions and attitudes (better than 5 cm and 0.5∘) in real time, using lightweight components. The main focus of this paper is on the attitude determination with the system. This attitude determination is based on an onboard single-frequency GPS baseline, MEMS (micro-electro-mechanical systems) inertial sensor readings, magnetic field observations and a 3D position measurement. All of this information is integrated in a sixteen-state error space Kalman filter. Special attention in the algorithm development is paid to the carrier phase ambiguity resolution of the single-frequency GPS baseline observations. We aim at a reliable and instantaneous ambiguity resolution, since the system is used in urban areas, where frequent losses of the GPS signal lock occur and the GPS measurement conditions are challenging. Flight tests and a comparison to a navigation-grade inertial navigation system illustrate the performance of the developed system in dynamic situations. Evaluations show that the accuracies of the system are 0.05∘ for the roll and the pitch angle and 0.2∘ for the yaw angle. The ambiguities of the single-frequency GPS baseline can be resolved instantaneously in more than 90% of the cases.


Introduction
The recent interest in unmanned aerial vehicles (UAVs) necessitates the development of small and lightweight direct georeferencing systems for the guidance, navigation and control (GNC) of UAVs and UAV-based surveying or remote sensing applications. Direct georeferencing comprises the determination of the position (e.g., X, Y , Z) and the attitude (e.g., α, β, γ) of a sensor or the vehicle relative to the Earth, represented in a predefined coordinate frame. The attitude is usually defined as the orientation between the body fixed coordinate frame and the navigation coordinate frame.
Although UAVs exist in many different size classes, this paper focuses on micro-and mini-sized UAVs, having a weight limit of 5 kg and a size limit of 1.5 m. Thus, the term UAV in this paper always relates to micro-and mini-sized UAVs.
In recent years, direct georeferencing systems have extensively been researched and used, especially in airborne applications [1,2]. However, these systems cannot easily be adopted for small UAVs, which is mainly due to weight and size constraints. Examples of the first approaches to the direct georeferencing of UAVs can be found in [3][4][5][6][7][8]. Nevertheless, in [3][4][5][6][7], direct georeferencing is not performed in hard real-time onboard the UAV. In particular, for the GNC and the realization of automatic operations, the real-time capability is very important. Falco et al. [8] describe the development of a real-time-capable low-cost system for the position and attitude determination of small UAVs. However, the system accuracies are only <2 m for the horizontal positions and <2 • for the yaw angle, which is not sufficient for many applications, e.g., surveying.
In this paper, a newly-developed direct georeferencing system is presented, which is intended to provide highly accurate positions and attitudes (σ pos < 5 cm, σ att < 0.5 • ) in real-time onboard a UAV, using lightweight components. Even if the system estimates both positions and attitudes [9], this paper is focused on the attitude determination with the developed system. The algorithms are described in detail, and the accuracy and the performance of the attitude determination are comprehensively evaluated.
The challenge in the attitude determination of UAVs comes with the usage of small and lightweight IMUs (inertial measurement unit), such as micro-electro-mechanical system (MEMS) IMUs, since they have a significantly lower quality than the navigation-grade IMUs, which are usually applied in direct georeferencing systems. Without the frequent availability of absolute attitude information, the output of the gyroscopes of MEMS IMUs would be significantly wrong after a few seconds. This applies in particular to the yaw angle. While an IMU itself measures the gravity, which in the long term can be used to bound the roll and the pitch error, the yaw angle error is not bounded, if only IMU readings are available and no assumption about the motion behavior is made. To overcome this problem, there are multiple options, which may serve as a yaw reference for a UAV attitude determination.
One option comes with the availability of position information in situations where the UAV's horizontal acceleration is not zero [10][11][12]. In this case, the yaw angle is correlated to the position observation. A drawback of this approach is clearly the need for non-zero accelerations.
Another option, which is often used in MEMS IMU-based attitude determination, is the integration of a magnetometer [12,13]. However, the main problem when using a magnetometer on a UAV is the distortion of the magnetic field measurements, caused by various time-constant and time-varying effects. Although most of the constant effects, coming from ferromagnetic material in the vicinity of the sensor, can be compensated by appropriate calibration methods [14,15], the non-constant effects, e.g., high electrical currents during the operation of rotary-wing UAVs, inhibit sub-degree attitude accuracies in most cases.
A third option is the realization of a GNSS multi-antenna system on the UAV to provide a direct yaw measurement. Examples for the use of multi-antenna GPS systems on lightweight UAVs can be found in [3,8,16]. Due to the low price and the small weight of some currently available single-frequency GPS receivers, this option is becoming more and more attractive. However, for the attitude determination, GPS carrier phases have to be observed, and a resolution of the carrier phase ambiguities is required. This is a challenging task, especially for low-cost single-frequency GPS observations [17]. During UAV applications in urban areas, frequent losses of the GPS signal lock can occur due to obstacles in the signal path. Since every loss of lock necessitates a re-initialization of the ambiguities, the ambiguity resolution process should be as fast as possible.
The contribution of this paper is the presentation and the evaluation of a small and lightweight real-time-capable direct georeferencing multisensor system, which also acts as an attitude and heading reference system (AHRS) for micro-and mini-sized UAVs. The desired attitude accuracy of the system is less than 0.5 • for all attitude angles. To achieve this and to overcome the disadvantages of the different yaw determination options, the following information is used and combined in a suitable way: • 3D position measurement: This can be a code-based GPS position with accuracies of several meters or a cm-accurate RTK (real-time kinematic) GPS position. We use an RTK GPS position in this paper. • MEMS IMU readings: The IMU readings include angular rates and accelerations, measured by three-axis gyroscopes and accelerometers. • Single-frequency GPS observations: The GPS observations, measured by two onboard GPS receivers, enable the determination of an onboard GPS baseline, which contributes to the yaw determination. • Magnetic field observations: The magnetic field observations come from a calibrated magnetometer, which is placed as far away as possible from the electric currents on the UAV platform.
In the following sections of this paper, we first explain the GPS attitude baseline determination algorithms, including the concept of an instantaneous ambiguity resolution for the single-frequency onboard GPS baseline. Afterwards, the details of the sensor integration and the system design will be presented. The evaluation of different field tests enables the assessment of the attitude accuracy and the ambiguity resolution performance of the system.

GPS Attitude Baseline Determination
It has been known for many years that GPS multi-antenna systems are well suited for the attitude determination of mobile objects [18,19]. The accuracy of a GPS attitude angle depends on the baseline length, which is limited to ca. 1 m for small UAVs, and the accuracy of the baseline coordinates. To achieve sub-cm accuracies of the baseline coordinates, which lead to sub-degree attitude accuracies for a 1-m baseline, carrier phases have to be observed, and the carrier phase ambiguities have to be fixed to their correct integer value. For kinematic applications, this ambiguity resolution should be fast and reliable at the same time.
In recent years, many different approaches have been investigated to allow for the ambiguity resolution in the attitude determination. For example, motion-based methods rely on the motion of the vehicle or changes in the receiver-satellite-geometry [19,20]. Many other procedures make use of the LAMBDA (Least-squares AMBiguity Decorrelation Adjustment) method [21]. Due to the integration of the known baseline length or geometry in the ambiguity objective function, this procedure has been improved for the attitude determination [22][23][24][25]. Further improvements are possible, when also inertial sensor readings are used to reduce the number of candidates in the ambiguity search space [16,[26][27][28][29][30][31]. For example, in Eling et al. [29], an instantaneous ambiguity resolution procedure based on the ambiguity function method (AFM) [32], aided by MEMS gyroscopes, is presented. Roth et al. [28] applied an enhanced approach of the extended LAMBDA method, which has been developed by Mönikes et al. [16]. To accelerate the ambiguity resolution, they use the baseline length and magnetic field observations as constraints. However, as described in [31], a fast or even instantaneous ambiguity resolution with a constrained LAMBDA method can be difficult, when the ambiguity float solution is systematically wrong, due to poor GPS measurements under challenging GPS measurement conditions in urban areas.
In our approach, we aim for an instantaneous and reliable single-frequency GPS ambiguity resolution during UAV flights. Therefore, we use the information from additional sensors (GPS position, inertial sensors and magnetic field observations) and the known baseline length to improve the ambiguity float solution. Afterwards, we perform a search in the ambiguity domain, using the modified LAMBDA method (MLAMBDA) [33]. For a further search and also a validation step, we perform a transformation into the coordinate domain, where the AFM and further checks are also applied. Figure 1 shows the overall process of the attitude determination with the developed direct georeferencing system. To aid the ambiguity estimation of the GPS baseline, first, an approximate attitude is estimated in a Kalman filter (Pos/IMU/Mag), using a position measurement (Pos) (in our case, an RTK GPS position), IMU readings and magnetic field observations (Mag) (Section 3). Afterwards, these attitude results serve as an input to the GPS baseline determination, which consists of the base float solution and the integer ambiguity resolution (Sections 2.2 and 2.3). If the ambiguities can be resolved successfully, the base fixed solution (Base) can also be integrated in the Kalman filter (Pos/IMU/Mag/Base).

Overview
The whole strategy assumes the availability of position measurements and GPS observations. If they are not available (no GPS observations), the sensor integration switches to a second mode (IMU/Mag), which only estimates the attitude based on inertial sensor and magnetometer readings (Section 4).
The coupling of the approximate information about the attitude and the GPS baseline raw data could generally also be realized in a tightly-coupled approach, but for reliability and practical reasons in the real-time programming, we decided to separate these steps here.  [33], and the AFM is the ambiguity function method [32].

GPS Baseline Float Solution
In the first step of the ambiguity resolution, the ambiguities are estimated as real values in the float solution. This estimation is realized in a total state space extended Kalman filter (EKF). Beside the single-difference (SD) [34] ambiguities N j AB , the state vector x SD of this filter also contains the baseline parameters ∆x e b and the first time derivative of the baseline parameters ∆ẋ e b : where A and B are the notations for both GPS antennas of the baseline and j is the notation for the satellites. The reason for estimating SD instead of double-difference (DD) [34] ambiguities is to avoid the hand over problem that would arise for DD ambiguities, when the reference satellite changes [35]. Measurement model: As observations DD carrier phases φ jk AB and DD pseudoranges P jk AB on the GPS L1 frequency (L1=1575.42 MHz) are used: where ρ jk AB (t k ) is the geometric range DD, λ is the GPS L1 wavelength and ϑ is the noise.
To improve the ambiguity resolution, the prior information about the attitude, known from the Pos/IMU/Mag solution (Section 3), is added to the observation vector. For this, the known body-frame (b-frame) baseline parameters ∆x b b are transformed into the e -frame (earth centered, earth fixed), using the C e b rotation matrix, which represents the attitudes from the Pos/IMU/Mag solution: Finally, the known baseline length s can be added as a constraint to the observation vector. Accordingly, the observation vector reads: A linearization of the measurement equations yields the measurement matrix: Therein, the matrix A includes the satellite/receiver geometry [34], and the matrix D is the double differencing matrix. n is the number of visible DD observations, and m = n + 1 is the number of visible satellites.
As the stochastic model for the GPS observations, the following equation is used: where f is a factor, which is f P = 100 for the pseudoranges and f φ = 1 for the carrier phases. el is the satellite elevation angle, and a and b are free parameters, weighting the influence of the elevation angle on the measurement noise. The choice of the values for a and b depends on various effects (antennas, receivers, antenna near-fields, etc.), which may be different for every measurement setup. In order to find realistic values representing the measurement noise of the GPS observations on our system during kinematic applications, we analyzed different flights. Based on the covariance matrices of the parameters and the residuals of the measurements, we found that a = 2 mm and b = 2 mm lead to the best results for the parameter estimation with our setup. Finally, in the complete measurement noise matrix R k , also the DD correlations are regarded: where σ ∆x b is the noise of the prior information and σ s the noise of the baseline length constraint. System dynamics model: The standard total state space EKF time update is the following: We use a simple random walk model as the system dynamics model. Thus, the state transition matrix is composed of: and G k is given by: The process noise matrix Q k is a diagonal matrix consisting of two types of process noise. The process noise for the first time derivative of the baseline parameters is set to σ ∆ẋ = 1 m/s. This choice reflects the approximate dynamic capabilities of the UAV, but does not smooth the baseline parameters too much. In contrast, the ambiguity parameters are assumed to be constant. However, to allow for a better reaction of the ambiguity parameters on changing satellite configurations, the noise parameter for the ambiguities is not set to zero, but to a very small value: σ N = 1 · 10 −4 cycles.
Since the initialization of the ambiguity parameters, which is based on the difference between the SD pseudoranges and the SD carrier phases, is very imprecise, the process noise for new initialized ambiguities is set to a very large value of σ N = 100 cycles. In this way, new initialized ambiguities do not have any influence on the estimation of the other parameters, but benefit from the current float solution.
In the final step of the float solution, the estimated SD float ambiguities and their covariance matrix have to be transformed to DD values, applying the D-matrix again.

Integer Ambiguity Resolution
The results of the float solution are real-valued DD ambiguities and their covariance matrix. However, the true ambiguities are integer valued, and they need to be found within the solution space, which is spanned by the estimated float ambiguities and their covariances. This process is basically a search problem, called integer ambiguity resolution, and it is the key to enable mm-accuracies for the baseline parameters. In our approach, three steps are performed (see also Figure 1) to resolve the ambiguities correctly: Step (i): In the first step of the integer ambiguity resolution, the MLAMBDA method [33] is applied. The MLAMBDA method is an ambiguity search technique in the ambiguity domain and an extension of the well-known LAMBDA method [21], leading to faster computation times [33]. The MLAMBDA method is based on an integer least squares problem to obtain the estimates of the DD integer ambiguities solving the minimization problem [21]: where N is the integer ambiguity vector,N is the float ambiguity vector and QN represents its covariance matrix. The inputs for the MLAMBDA method are the DD float ambiguities and their covariance matrix.
Due to the aid of the additional information about the attitude, known from the Pos/IMU/Mag solution (Section 3), the baseline float solution is more accurate, compared to a GPS-only float solution. The more accurate the float ambiguities are, the more likely is a fast and correct ambiguity resolution with the MLAMBDA method. Thus, there is a good chance that the ambiguities can already be fixed at this step of the ambiguity resolution process in our approach. As outputs of the MLAMBDA method, we consider ten sets of ambiguity candidates. To make a decision if any of these sets can be fixed, we first apply the simple ratio test [36], where the squared norm of the ambiguity residuals of the best set of ambiguities (R 1 ) and the second best set of ambiguities (R 2 ) is compared. Only if the quotient R 2 /R 1 exceeds a threshold of three, the ambiguities of the best set will be fixed to integer values.
Step (ii): In case the ambiguities could not be fixed within Step (i), further effort is needed. The ten best solutions from the MLAMBDA method are considered further. In the first instance, we try to perform an elimination of single ambiguity parameters, which probably prevent the fixing of all of the other ambiguities. For example, in the case that all ten sets of ambiguities are equal except the ambiguity of one DD observation, this one varying ambiguity parameter is omitted. In the case that more than one ambiguity parameter is varying in the ten sets of ambiguities, always one ambiguity parameter is eliminated in a random order, and the MLAMBDA method is applied again, until the ratio test is positive or a maximum number of iterations is reached. Due to the elimination of single ambiguity parameters, some of the ten sets of ambiguity candidates can be omitted now, since they probably became equal to other ambiguity sets.
Step (iii): In the last step of the ambiguity estimation, a transformation into the coordinate domain is performed. This means that all of the remaining sets of ambiguity candidates are used to determine baseline parameters. Afterwards, these baseline parameters are tested in the ambiguity function method (AFM) [32], which is a search technique based on a trigonometric cost function, the ambiguity resolution function (ARF): By use of the determined baseline parameters, all observed DD carrier phases φ jk obs can also be calculated (φ jk calc ). If the determined baseline parameters are correct, the differences between the observed and the calculated DD carrier phases are the integer ambiguities. Therefore, the result of the ARF should be close to the number of the available DD carrier phases (n) if the correct candidate is tested. Beside the AFM, we also consider the lengths and the residuals of the determined baseline parameters to check the correctness of the ambiguities.
As can be seen in Figure 1 Step (iii) of the ambiguity resolution process is also applied to fixed solutions from Step (i) or (ii), to validate these results. Thus, due to the combination of Steps (i), (ii) and (iii), a fast and reliable ambiguity resolution can be performed at the same time.

Fixed Solution
In the fixed solution, the final baseline parameters ∆x e b are estimated. Once the ambiguities have been fixed successfully, they can generally be held fixed for the further epochs, as long as no signal interruption or cycle slip occurs. To ensure that the ambiguities are held fixed correctly, they are resolved independently, epoch by epoch, until the same ambiguity resolution has been determined in ten consecutive epochs. The risk that the same incorrect ambiguity resolution is found in ten consecutive epochs during dynamic applications is negligible.

Cycle Slip Detection
A loss of a satellite signal lock leads to a reinitialization of the corresponding integer counter in the GPS receiver. As a result, the previously known ambiguity parameter of the affected satellite is not valid any longer. These cycle slips have to be detected, if the ambiguities are held fixed. The cycle slip detection can be done using a Kalman filter [37,38], where an exceptional change in the carrier phases between two measurement epochs will be conspicuous, when a good prediction of the baseline parameters is available. These requirements are met here, since the motion behavior of the baseline parameters between two epochs is well known from the sensor integration.
The transfer of the information from the sensor integration to the cycle slip detection is realized using the attitude change, which can be converted into the first deviation of the baseline parameters ∆x e b . Finally, this information serves as control input, to check the consistency of the ambiguity parameters in a second GPS Kalman filter, which is decoupled from the float solution.
The state vector x C of this simple filter includes the baseline and the ambiguity parameters. The system dynamics model is a random walk model: The observation vector includes the carrier phases and the pseudoranges: and the measurement matrix represents the nonlinear functional model for a single baseline: Cycle slips can now be detected considering the innovation: In case no cycle slips exist, the innovation of the carrier phases remains close to zero. When the innovation of one DD carrier phase exceeds a threshold d > 0.1 m, a cycle slip is probably the reason, and the corresponding ambiguity parameter is re-initialized. If at least three DD carrier phases show inaccuracies, the complete ambiguity set is rejected, and the ambiguity resolution process is re-initialized.

The Sensor Integration
The sensor integration can be separated into the strapdown algorithm (SDA) and the Kalman filter update. In the SDA, the highly dynamic movement of the system is determined integrating the angular rates and accelerations of the MEMS IMU in real time. As a result of various errors (initial alignment, inertial sensors, computational, etc.), the SDA solution is drifting over time. Therefore, the additional sensors, providing long-term, stable measurements, are needed to correct and bound the drift of the inertial sensor integration [39]. The estimation of the SDA errors is realized in an error state space Kalman filter. Some advantages for this approach can, for example, be found in [40,41].
In the following, the SDA and the Kalman filter update will be summarized. More details can be found in [42][43][44]. Here, the navigation equations of the b-frame are expressed in the e-frame. The b-frame to e-frame rotation C e b is split into a rotation of the b-frame with respect to the local-level navigation frame (n-frame) C n b and a rotation of the n-frame with respect to the e-frame C e n . Since the axes of the n-frame are aligned with the directions east, north and up, the C e n rotation is: where λ o and ϕ o are the ellipsoidal coordinates of the b-frame origin. The rotation from the b-frame to the n-frame may be expressed as the product of three separate rotations: using the Euler angles roll α, pitch β and yaw γ.
In the full state vector x, not the Euler angles, but a quaternion q is used for the attitude representation. Further states are the position x e p , the velocity v e eb , the accelerometer bias b b a and the gyro bias b b ω :

Strapdown Algorithm
Within the SDA, the MEMS IMU readings are used to propagate the attitude, the velocity and the position of the system.
Attitude: The gyroscopes of the IMU measure the attitude change ω b ib of the b-frame with respect to an inertial frame represented in the b-frame. Since the Earth rotation rate ω e ie is small compared to the measurement noise of the MEMS gyroscopes, the approximation: is valid. Assuming that the direction of the angular rate vector ω b eb remains fixed in space during the time interval ∆t = t k − t k−1 , the rotation angle ∆σ k can be determined with the trapezoidal rule: With the quaternion ∆q k : the current attitude q k results from: where q k represents the quaternion relating the b-frame to e-frame axes at the epoch k.

Velocity: The velocity propagation is based on the integration of the measured accelerations a b ib
: The assumption of a constant direction of the acceleration vector during the time interval ∆t is also valid, with the result that Equation (27) can be simplified: Then, the velocity at time k can be determined applying Equation (29).
v e eb,k = v e eb,k−1 + C e b a b ib + The vector g e l is the local gravitation vector represented in the e-frame. Position: Finally, the current position may be derived integrating the velocity. For this integration, again, the trapezoidal rule is used, since the integration time is short here:

Kalman Filter Update
The Kalman filter update is based on an error state extended Kalman filter, which means that the estimated state is not the full state vector itself, as described in Equation (22), but the error vector between the SDA predicted state and the full state. To estimate this error, the following measurements are used: System dynamics model: The derivation of the error dynamics equations can, for example, be found in [42] or [45]. Hence, the linearized continuous system dynamics model specialized to the e-frame is: where 0 is always a 3 × 3 zero matrix and I a 3 × 3 identity matrix. [â e ib ×] denotes a 3 × 3 skew symmetric matrix with elements of the bias-corrected acceleration vector a b ib , which has been rotated into the e-frame. The vector n includes the stochastic variables.
Measurement model: For the integration of the GPS positions the lever arm l b between the GPS antenna reference point x e a and the system reference point x e p , which can be identical to the IMU reference point, has to be regarded [45]: x e a = x e p + C e b l b + ϑ xa = x e p + (I + Ψ)Ĉ e b l b + ϑ xa (32) = x e p + (I + Ψ)l e + ϑ xa = x e p +l e − l e × ψ + ϑ xa where ϑ xa denotes the GPS position measurement noise. Due to the coupling between the GPS position and the measured accelerations from the IMU, the yaw angle error is bounded in the presence of horizontal accelerations. However, in the absence of horizontal accelerations, for example when the UAV is hovering, the yaw angle is still unobserved [12].
To overcome the problem of an undetermined yaw angle, the GPS attitude baseline is arranged on the UAV. In the measurement model, the link between the GPS baseline and the attitude determination is given by the transformation of the b-frame GPS baseline coordinates ∆x b b , which can be determined by calibration measurements, into the e-frame: where ϑ x b is the measurement noise of the GPS attitude baseline. Finally, as long as the ambiguities of the GPS baseline are unresolved, the baseline parameters do not have the quality to correct the yaw error adequately. To be able to also bound the yaw error in these cases, we integrated magnetic field sensors into our system.
Magnetic field sensors, when used as a compass, measure two or three components of the Earth's magnetic field in the sensor system. Therefore, they enable the determination of the orientation of the sensor system relative to the Magnetic North Pole. However, the main difficulties with magnetic field sensors are their sensitivities to ferromagnetic material in the vicinity of the sensor. In case this material is part of the sensor platform, disturbances can be compensated by calibration [15]. Nevertheless, the changing currents of the rotors or ferromagnetic material in the environment can still lead to distortions of the magnetic field observations. This is why the attitude accuracy, based on magnetometer measurements, is usually only in the range of a few degrees, which is still good enough to bound the gyro drift.
The measurement equation for a magnetic field measurementh b is the following [12]: where h e is the known local magnetic field vector, represented in the e-frame, and ϑ m is the magnetometer measurement noise.
Since erroneous magnetic field observations can lead to significant systematic attitude errors in all three axis, Equation (35) is rearranged, with the result that the magnetic field observations only have an impact on the yaw angle determination: With the ambiguities fixed, the parameters of the baseline can be determined with a mm-accuracy. For a baseline with a length of 1 m, this leads to attitude accuracies of approximately 0.2-0.5 • for the yaw angle. Since the roll and pitch angle can be determined more reliably regarding the relation to the gravitation vector, the GPS baseline should also only be used to estimate the yaw angle error: In case the ambiguities of the GPS baseline are unresolved so far, only the GPS positions and the magnetic field observations are used in the Pos/IMU/Mag integration (Figure 1), which serves as the prior information for the GPS baseline ambiguity resolution. If the ambiguities are resolved successfully, the Pos/IMU/Mag/Base solution can be determined, where all available sensory input is integrated.
Kalman filter formulation: In the prediction step of the Kalman filter, only the error covariance matrix is propagated: Therein, the matrices Φ k and G k are the discretizations of the matrices F and G at the time k. The matrix Q k−1 includes the discrete time noise.
The filter update is computed using the standard Kalman filter equations: where K k is the Kalman gain matrix, R k is the measurement noise matrix, H k is the measurement matrix, z k is the observation vector and y k is the vector of the predicted observations, which are based on the results of the SDA. For the update of the full state vector, the error state vector is simply added to the SDA-derived state vector. For the attitude update, Equations (25) and (26) have to be applied using δψ instead of ∆σ.

Special Case: GPS Unavailability
When GPS observations have not been available for a few seconds, the Kalman filter, as it was presented in Section 3, cannot be applied any longer. The reasons for this are that the position and velocity errors, resulting from the inertial sensor integration, cannot be bounded without the GPS observations. Furthermore, the accelerometer biases are no longer observable longer. On this account, the system is intended to stop providing position and velocity estimates, and the Kalman filter update switches to a second mode [12], which is referred to as the IMU/Mag solution in Figure 1.
The Kalman filter in the second mode is a subset of the previously-presented filter (Section 3). The attitude prediction is still determined integrating the angular rates, but the state vector only consists of the attitude and the gyro bias errors. Since the magnetic field observations are not enough information to observe the full orientation (the rotation around the Earth's magnetic field axis would be undetermined), the accelerometer readings are added to the measurement model. The measurement equation for the accelerations is the following [45]: It is assumed that the main force measured by the accelerometers is the local Earth gravitation vector g e l . The influence of the trajectory dynamics on the accelerometer measurements is modeled by a large measurement noise with ϑ a = 3 m/s 2 .

System Realization
The developed direct georeferencing system is shown in Figure 2. Its dimensions are 11.0 cm × 10.2 cm × 4.5 cm, and its weight is 240 g, without the GPS antennas.
The system consists of a dual-frequency geodetic-grade GPS receiver (Novatel OEM 615), which acts as the main positioning device. Together with the GPS raw data of a master station, which is transmitted via a radio module (XBee Pro 868), it allows for an RTK GPS (real-time kinematic) positioning with a rate of 10 Hz, leading to centimeter position accuracies.
The dual-frequency GPS receiver also serves as one receiver of the onboard GPS baseline for the GPS attitude determination (1 Hz). However, for the attitude determination, only the L1 GPS frequency is used, since the second onboard GPS device is a low-cost single-frequency GPS receiver (Ublox LEA6T).
The tactical grade IMU (Analog Devices Adis16488), including three-axis gyroscopes and accelerometers, provides high rate inertial data (100 Hz), which is used for the position and attitude determination. A magnetometer (Honeywell HMC5883L) is placed at the outer end of one of the rotor-free UAV arms, where it is not affected by electric currents.
A real-time processing unit (National Instruments sbRIO 9606), which is a reconfigurable IO board, including an FPGA (field programmable gate array) and a 400-MHz processor, enables the onboard position and attitude determination in real time. Figure 3 shows the current version of the UAV platform, which has been developed within a project called Mapping on Demand [46]. In this project, the UAV is intended to fly fully autonomously. Therefore, additional sensors, such as two stereo camera pairs and a 5 MPixel camera, are also fixed on the platform.  The GPS antennas of the attitude baseline can be seen on the left and the right side of Figure 3. The baseline length is 92 cm. Both antennas are geodetic-grade antennas (3G + C, navXperience), which are originally too heavy for UAV applications. In order to reduce weight, we omitted all unnecessary components of the antennas, such as the antenna housings and the 5/8" screw threads. Due to the changes in the antenna design, we had to recalibrate the antennas in our anechoic chamber [9,47]. By comparison to the original antenna, the antenna modifications led to significant changes in the phase center offsets (approximately 4 cm in the up, <1 mm in the north and east component) and in the phase center variations (<5 mm) of the antennas.

Results
In this section, the ambiguity resolution performance and the attitude accuracy of the direct georeferencing system will be evaluated. For the evaluation of the ambiguity resolution, real flight test data are used. Since it is difficult to get a highly accurate attitude reference for UAV flights, a test setup has been applied for the accuracy evaluation. On this setup (Figure 4), the direct georeferencing system, the GPS antennas and the magnetometer are mounted in the same configuration as on the UAV platform. As the reference system, a navigation-grade INS (inertial navigation system) from Imar has been used, which consists of fiber optical gyros (FOG) and provides attitudes with a rate of 1000 Hz. The specifications of the Imar INS are shown in Table 1, together with the specifications of the MEMS IMU (Adis16488), which is used on the direct georeferencing system. Due to its high level of bias stability (<0.003-0.01 • /h) and the low influence of the angular random walk (0.001 • / √ h), the attitudes of the Imar INS serve as an excellent attitude reference for about one hour after the measurement start.

Ambiguity Resolution
Concerning the evaluation of the ambiguity resolution, there are two questions we would like to answer. To make this possible, data of three flight tests (Flights 1-3) and three measurements with the previously-presented test setup (Tests 1-3) are used. All of these measurements were performed on different days and at different locations.
1. How often does the ambiguity resolution algorithm succeed in fixing the ambiguities in the first epoch (instantaneously) after a loss of lock? Table 2 shows single-epoch ambiguity resolution success rates. For calculating these values, the complete integer ambiguity estimation process was re-initialized in every single epoch of the datasets. In the presentation of these values, a distinction is made between the three steps (Steps i-iii) of our ambiguity resolution algorithm (Section 2.3), on the one hand, and a standard ambiguity resolution algorithm, on the other hand. This standard ambiguity resolution algorithm includes a float solution without any prior information, an integer ambiguity estimation with the MLAMBDA method and a ratio test with a threshold of three. It can be seen that our algorithm leads to instantaneous ambiguity resolutions in more than 90% of the epochs. Compared to the standard approach, this is a significant improvement, especially during the UAV flights, where the GPS measurement conditions have been more challenging, due to buildings and vegetation in the surrounding areas.
Within the different steps of our algorithm, a significant increase of the success rate can be seen between Step (i) and Step (ii), whereas the difference between Step (ii) and Step (iii) is smaller for the single-epoch ambiguity resolutions. In Flight 1, the success rate after Step (iii) even is smaller than after Step (ii). The reason for this is that some fixed solutions from Step (i) or (ii) have been rejected in the validation of Step (iii).

What is the average time to fix the ambiguities?
To answer the second question, the integer ambiguity estimation process was restarted in every epoch again, but now, the algorithm always proceeded for as many epochs as needed to resolve the ambiguities. In this way, each time to first fix could be measured. Figure 5 presents the percentage of cases in which the ambiguities could be resolved after a certain time. Here, also, a comparison is made between the standard and our improved ambiguity resolution algorithm. Accordingly, with our algorithm, the ambiguities could at least after 5 s be resolved in more than 99% of the cases during UAV flights. In contrast, with the standard algorithm, the ambiguities could after 5 s only be resolved in 20% of the cases during the UAV flights. Figure 5. Ambiguity resolution (Amb. res.) success rates after a certain time, determined with a standard ambiguity resolution algorithm (standard) and our ambiguity resolution algorithm (improved) for flight test data (Flight) and test setup data (Test).

Accuracy
For the accuracy evaluation, the dataset "Test 2" is used in the following, which includes static and kinematic parts. During the kinematic part, the measurement setup was carried by two persons, trying to simulate typical UAV dynamics.
In Figure 6, the differences between the method described here and the reference attitude (Imar INS) are shown. It can be seen that the maximum difference is less than 0.5 • for all attitudes, where the roll and the pitch angle seem to be more accurate than the yaw angle. The mean, the root mean square (rms) values and the maximum deviations (max), which were determined on the basis of the shown differences, are presented in Table 3. The rms for the roll and the pitch are about 0.03 • in the static mode and 0.05 • in the kinematic mode. The rms for the yaw angle is about 0.2 • for a static and a kinematic mode, where the deviations of the yaw angle are mainly caused by the systematic errors of the GPS attitude baseline.
In Figure 7, the Pos/IMU/Mag/Base yaw angle differences from Figure 6 are compared to the GPS baseline attitude differences and the differences of the Pos/IMU/Mag solution. The maximum deviation of the Pos/IMU/Mag solution lies in the first static part of the measurement, where the yaw angle determination has mainly been influenced by the magnetic field observations. Even if these deviations are in the order of a few degrees, the solution still serves as good prior information for the ambiguity search. From experience, this is the same for UAV flights, where the deviations between the prior information and the GPS attitudes usually have a comparable magnitude to the deviations seen here. Due to the coupling with the GPS positions (we used RTK GPS positions), in combination with horizontal accelerations, the deviations decreased during the kinematic part of the measurement.
Although the ambiguities were held fixed for the GPS base solution (Figure 7), the cycle slip detection led to a re-initialization of the ambiguities in two epochs during this test. Furthermore, some GPS base outliers are also conspicuous (e.g., after ca. 4.2 min). However, the outliers do not deteriorate the attitudes from the Pos/IMU/Mag/Base integration.

Conclusions
In this paper, a direct georeferencing system has been presented, which is intended to provide highly accurate positions and attitudes (better than 5 cm and 0.5 • ) of a micro-and mini-sized UAV in real time. The paper has mainly been focused on the attitude determination with this system. For the attitude determination MEMS IMU readings, magnetic field observations, a GPS baseline and a GPS position measurement are integrated in a sixteen-state Kalman filter. To be able to resolve the ambiguities of the single-frequency GPS baseline reliably and instantaneously, the float solution has been improved. For the integer ambiguity estimation, the MLAMBDA method and the AFM are combined in a suitable way.
With a test measurement, where a navigation-grade INS served as a reference, it could be shown that the implemented algorithms allow for attitude accuracies of about 0.05 • for the roll and the pitch and 0.2 • for the yaw angle during kinematic applications. With the developed ambiguity resolution algorithms, the ambiguities of the single-frequency GPS attitude baseline could be resolved instantaneously in more than 90% of the epochs. It could also be shown that short-term losses of the GPS baseline can be bridged successfully by the inertial sensors without a deterioration of the attitude accuracy.
As a final conclusion, the use of a low-cost single-frequency GPS baseline can significantly improve the attitude determination of micro-and mini-sized UAVs, even in non-ideal GPS measurement environments.

Acknowledgments
This work was funded by the DFG (Deutsche Forschungsgemeinschaft) under Project Number 1505 "Mapping on Demand".

Author Contributions
Christian Eling developed the algorithms, implemented the software, wrote the manuscript and contributed to the development of the direct georeferencing system. Lasse Klingbeil assisted in the algorithm development and the writing of the manuscript. He also contributed to the development of