GNSS Carrier Phase-Based Attitude Determination

The GNSS (Global Navigation Satellite Systems) are a valid aid in support of the aeronautic science. GNSS technology has been successfully implemented in aircraft design, in order to provide accurate position, velocity and heading estimations. Although it does not yet comply with aviation integrity requirements, GNSS-based aircraft navigation is one of the alternative means to traditional dead-reckoning systems. It can provide fast, accurate, and driftless positioning solutions. Additionally, ground-based GNSS receivers may be employed to aid navigation in critical applications, such as precision approaches and landings.

Introduction
The GNSS (Global Navigation Satellite Systems) are a valid aid in support of the aeronautic science. GNSS technology has been successfully implemented in aircraft design, in order to provide accurate position, velocity and heading estimations. Although it does not yet comply with aviation integrity requirements, GNSS-based aircraft navigation is one of the alternative means to traditional dead-reckoning systems. It can provide fast, accurate, and driftless positioning solutions. Additionally, ground-based GNSS receivers may be employed to aid navigation in critical applications, such as precision approaches and landings.
One of the main issues in airborne navigation is the determination of the aircraft attitude, i.e., the orientation of the aircraft with respect to a defined reference system. Many sensors and technologies are available to estimate the attitude of a aircraft, but there is a growing interest in GNSS-based attitude determination (AD), often integrated at various levels of tightness to other types of sensors, typically Inertial Measurements Units (IMU). Although the accuracy of a stand-alone GNSS attitude system might not be comparable with the one obtainable with other modern attitude sensors, a GNSS-based system presents several advantages. It is inherently driftless, a GNSS receiver has low power consumption, it requires minor maintenance, and it is not as expensive as other high-precision systems, such as laser gyroscopes.
GNSS-based AD employs a number of antennas rigidly mounted on the aircraft's structure, as depicted in Figure 1. The orientation of each of the baselines formed between the antennas is determined by computing their relative positions. The use of GNSS carrier phase signals enables very precise range measurements, which can then be related to angular estimations. However, carrier phase measurements are affected by unknown integer ambiguities, since only their fractional part is measured by the receiver. The process of reconstructing the number of whole cycles from a set of measurements affected by errors goes under the name of ambiguity resolution (AR). Only after these ambiguities are correctly resolved to their correct integer values, will reliable baseline measurements and attitude estimations become available. This chapter focuses on novel AR and AD methods. Recent advances in GNSS-based attitude determination have demonstrated that the two problems can be formulated in an integrated manner, i.e., aircraft attitude and the phase ambiguities can be considered as the unknown parameters of a common ambiguity-attitude estimation method. In this integrated approach, the AR and AD problems are solved together by means of the theory of Constrained Integer Least-Squares (C-ILS). This theory extends the well-known least-squares theory (LS), by having geometrical constraints as well as integer constraints imposed on parameter subsets. The novel AR-AD estimation problem is discussed and its various properties are analyzed. The method's complexity is addressed by presenting new numerical algorithms that largely reduce the required processing load. The main objective of this chapter is to provide evidence that: • GNSS carrier-phase based attitude determination is a viable alternative to existing attitude sensors • Employing the new ambiguity-attitude estimation method enhances ambiguity resolution performance • The new method can be implemented such that it is suitable for real-time applications The structure of this contribution is as follows. Section 2 gives the observation and stochastic model which cast the set of GNSS observations, with special focus on the derivation of the GNSS-based attitude model. Section 3 reviews the most common attitude parameterization and estimation methods, mainly focusing on those widely used in aviation applications. Section 4 introduces a new ambiguity-attitude estimation method, which enhances the existing approach for attitude determination using GNSS signals. Section 5 presents flight-test results, which provide practical evidence of the novel method's performance. Finally, section 6 draws several conclusions.

The GNSS-based attitude model
A GNSS receiver works by tracking satellites in view and storing the data received. Each GNSS satellite broadcasts a coded message with information about its orbit, the time of transmission, and few other parameters necessary for the correct processing at receiver side (Misra & Enge, 2001). By collecting signals from three or more satellites a GNSS receiver determines its own position with a triangulation procedure, exploiting the knowledge about both the satellites positions and the slant distance (range) by each satellite in view. The range measurements are obtained by detecting the time of arrival of the signal, from which the range can be inferred. This measurement is affected by several error sources: the satellite and receiver clocks are not perfectly synchronized; the signal travels through the atmosphere, which causes delays; the direct signal may be affected by unwanted reflections (multipath) that cannot be perfectly eliminated by careful antenna design. If not properly modeled, each of these effects will limit the achievable GNSS accuracy. The observed pseudorange or code observable is therefore modeled as where the superscript s indicates the satellite and the subscripts r and f indicate the receiver and the frequency, respectively. The different terms are: The magnitude of errors involved in these observations -decimeter or meter level -would not allow high-precision applications, such AD, which require cm-or mm-level accuracy in the final positioning product. Therefore, another set of observations is considered: the phase of the tracked signal, modeled as with ϕ the phase of the generated carrier signal (original or replica) in cycles, t 0 t h et i m eo f reference for phase synchronization, and λ f the wavelength of frequency f . The phase reading is characterized by different atmospheric delays (the ionosphere causes an anticipation of phase instead of a delay), different instrumental biases (indicated with δ), different multipath and an additional bias which is represented by the unknown number of whole cycles that cannot be detected by the tracking loop, since only the fractional part is measured. These are the integer ambiguities z. In case of GNSS, the precision of the phase measurements far exceeds the one of code observations: typically the phase observable is two orders of magnitude more accurate than the code measurement.
The many sources of error in (1) and (2) can be mitigated in relative positioning models. First, we form the so-called single difference (SD) code and carrier phase observations by taking the differences between observations simultaneously collected at two antennas tracking the same satellite: (3) where subscript r 12 indicates the difference between two antennas: r 12 = r 2 − r 1 .T h e phase value ϕ s f (t 0 ), relative to the common satellite, is eliminated. The instrumental delays and clock errors of the satellite are usually considered constant over short time spans, since the travel time difference with respect to any two points on the Earth surface is small (Teunissen & Kleusberg, 1998).
The terms cdt r 12 , cd r 12 , f and δ r 12 , f refer to the relative clock errors and relative instrumental delays between the two receivers. A perfect synchronization between receivers implies the cancellation of the clock biases, and a correct calibration would reduce the impact of instrumental delays. In the case of a single receiver connected to two antennas, these two sources of relative error could cancel out with a proper calibration.
The receiver clock errors and hardware delays in the single difference equations (3) are common for all the satellites tracked at the same frequency. Therefore these terms can be eliminated by forming a double difference (DD) combination, obtained by subtracting two SD measurements from two different satellites: P s 12 r 12 , f = ρ s 12 r 12 , f + I s 12 r 12 , f + T s 12 r 12 + dm s 12 r 12 , f + ε s 12 P,r 12 , f Φ s 12 r 12 , f = ρ s 12 r 12 , f − I s 12 r 12 , f + T s 12 r 12 + δm s 12 r 12 , f + λ f z s 12 r 12 , f + ε s 12 It has been assumed that the real-valued initial phase of the receiver replica does not vary for different tracked GNSS satellites.
The differential atmospheric delays depend on the distance between antennas. For sufficiently short baselines -typically shorter than a kilometer -the signals received by the antennas have traveled approximately the same path, thus the atmospheric delays becomes highly correlated. The differencing operation makes these errors negligible with respect to the measurement white noise for the baselines typically employed in AD applications, which rarely exceeds a few hundred meters.
Note that the relation between observations and baseline coordinates is nonlinear, since these are contained in the range term with r s and r r the satellite and receiver antenna position vectors, respectively. By assuming the atmospheric delays negligible and applying the Taylor expansion to expression (4) one obtains the linearized relations △P s 12 r 12 , f = −(u s 12 r ) T △r r 12 + ε s 12 P,r 12 , f △Φ s 12 r 12 , f = −(u s 12 r ) T △r r 12 + λ f z s 12 where the observables are now 'observed minus computed' terms, and the unknowns are expressed as increments with respect to a computed approximate value. △r r 12 is the baseline vector -the difference between the absolute antennas positions -whereas u s 12 r = u s 2 r − u s 1 r is the difference between unit line-of-sight vectors of different satellites. Also note that the multipath terms have been lumped into the remaining unmodeled errors ε s 12 P,r 12 , f and ε s 12 Φ,r 12 , f . Consider now two antennas simultaneously tracking the same m + 1 satellites at N frequencies. The vector of DD observations of type (6) are cast in the linear(ized) functional model (Teunissen & Kleusberg, 1998) with y the 2mN-vector of code and carrier phase observations, z the unknown integer-valued ambiguities and b the vector of real-valued baseline coordinates. A and G are the design matrices with Λ the diagonal matrix of N carrier wavelengths and U the m × 3 matrix of DD unit line-of-sight vectors. Symbol ⊗ denotes the Kronecker product (Van Loan, 2000). Model (7) describes the linear relationship between GNSS observables and the parameters of the two antennas. However, a single baseline is generally not sufficient to estimate the full orientation of an aircraft with respect to a given reference frame. At least three non-aligned antennas are necessary to guarantee that each rotation of the aircraft can be tracked unambiguously. It is straightforward to generalize the model formulation (7) to cast n DD baseline observations, obtained with n + 1 GNSS antennas (Teunissen, 2007a): This formulation is obtained by casting the observations at each baseline in the columns of the 2mN × n matrix Y.C o n s e q u e n t l y ,Z = [z 1 ,...,z n ] is the matrix whose n columns are the integer ambiguity mN-vectors, and B = [b 1 ,...,b n ] is the 3 × n matrix that contains the n real-valued baseline vectors. We exploited here once again the short baseline hypothesis: the same matrix of line-of-sight vectors U is used for all baselines.
Besides describing the functional relationship between observables and unknowns, a proper modeling should also capture the observation noise, i.e., the measurement error. The error is relative to the receiver, to the satellite, to the frequency and to the type of observations (code or phase). The variance-covariance (v-c) matrix of a vector of DD observations y collected at baseline b will be denoted as D(y)=Q yy ,w i t hD(·) the dispersion operator. For the multibaseline model (9), the description of measurement errors requires a further step: the observations are cast into a 2mNn vector by applying the vec operator, which stacks the columns of a matrix. The v-c matrix Q YY that characterizes the error statistic of vec(Y) is A simple expression for Q YY is obtained by assuming that each of the baselines is characterized by the same v-c matrix Q yy : with P n the n × n matrix that takes care of the correlation which is introduced by having a common antenna: Expressions (9) and (11) define the GNSS multibaseline model that we use in this contribution as the foundation of our GNSS-based attitude estimation theory.
With the available code and phase observations it is possible to estimate the set of baseline coordinates. These can then be used to provide the aircraft attitude, but only when a further condition is realized: the positions of the antennas installed aboard the given aircraft are known, rigid and do not change over time (or, if change occurs, it is perfectly known and predictable). This is so because it is necessary to have a one-to-one relationship between aircraft attitude and baselines attitude. As an example, consider two antennas mounted on the two extremities of a flexible mast: it is not possible to separate the rotations of the mast from its deformations by only observing the variations of the mutual position between the two antennas.
The rigidity assumption is formalized in the following way. Consider two orthonormal frames, defined by the basis {u 1 , u 2 , u 3 } and {u ′ 1 , u ′ 2 , u ′ 3 }. Let us assume that the second frame is integrally fixed with the aircraft. An arbitrary vector x can be equivalently described by using either reference system: The relation between the components of vectors x and x ′ is completely defined by the mutual orientation of the two reference systems. The linear transformation x = Rx ′ allows for a one-to-one relationship. Matrix R, hereafter referred to as rotation matrix or attitude matrix, belongs to the class of orthonormal matrices O: its column vectors r i are normal and their product null: r T i r j = δ i,j ,w i t hδ i,j the Kronecker's delta (δ i,j = 1i fi = j,0o t h e r w i s e ) . These constraints are necessary for the admissibility of transformation x = Rx ′ .I na b s e n c eo f deformations, the scalar product between any two vectors should be invariant with respect to the transformation: whereas the vectorial product is invariant under rotations about the axis defined by x ′ × y ′ : These conditions are fulfilled for orthonormal rotation matrices with determinant equal to one.
Model (9) can then be reformulated by means of the linear transformation B = RF,whereF is used to cast the set of known local baseline coordinates and R is the orthonormal (R T R = I q ) matrix that rotates B into F. The complete GNSS attitude model reads then (Teunissen, 2007a;) Parameter q is introduced in order to make model (16) of general applicability. The n baselines may be aligned or coplanar, impeding the estimation of a full 3 × 3 matrix R. Therefore, q defines the span of matrix F. For baseline sets formed by aligning n + 1 antennas we set q = 1, whereas configurations of coplanar antennas are defined by q = 2. With four or more non-coplanar antennas, q = 3.
The GNSS attitude model (16) is a nonlinear model. Although the relation between observables and unknowns remain linear, the orthonormal constraint is of a nonlinear nature, and profoundly affects the estimation process. This is investigated in section 4. First, the following section gives an overview of common attitude parameterization and estimation methods.

Attitude parameterization and estimation
The orthonormality of R (R T R = I q )imposes q(q+1) 2 constraints on its components r ij .Th efu ll matrix R can then be parameterized with a properly chosen set of variables, whose number c a nb ea sl i t t l ea st w o( i fq = 1) or three (if q ≤ 2). To this purpose, several representations may be used, and few are briefly reviewed in the following.
From a set of code and phase observations cast as in (16), the problem of extracting the components of the attitude representation involves, as shown in section 4, a nonlinear least squares problem. Its formulation and solution are the second topic discussed in this section.

Attitude parameterization
Several attitude parameterizations are available in the literature, see e.g., Shuster (1993) and references therein. The most common parameterizations are briefly reviewed in the following.

Direction cosine matrix
The transformation between two basis of orthonormal frames reads with r ij the entries of R. The scalar product between any two unit vectors of the two frames is Hence, the attitude matrix can be expressed by nine direction cosines, i.e., the nine cosines of the angles formed by the three unit vectors of the first frame and the three unit vectors of the second frame: This representation fully defines the mutual orientation of the two frames, by using a set of nine parameters (see Figure 2). Each configuration can be described without incurring any singularity, at the cost of having a larger number of parameters than other representations. The main axis u ′ 1 is completely defined by the knowledge of the three direction cosines u T 1 u ′ 1 , u T 2 u ′ 1 and u T 3 u ′ 1 .

Euler angles
Consider counterclockwise rotations about one of the main axis of a frame {u 1 , u 2 , u 3 }.T h e n the rotation matrix R is obtained through one of the following expressions: Any arbitrary rotation can always be decomposed as a combination of three consecutive rotations about the main axis u 1 , u 2 or u 3 , represented by one of the relations in (20). Figure  3 shows the example of a 321 rotation: the first rotation is about the third main axis u 3 with magnitude ψ, the second is about the (new) second main axis u ′ 2 with magnitude θ,t h el a s t about the (new) first main axis u ′′ 1 with magnitude φ. The rotation matrix that defines the transformation between the frames {u 1 , u 2 , . Twelve combinations of rotations are possible, whose choice depends on the application. As an example, the sequence 321 is commonly used to describe the orientation of an aircraft, where the angles ψ, θ, φ are named heading, elevation and bank, respectively. It is easy to see that the Euler angles representation is not unique: e.g, the combination 321 is equivalently expressed as The main advantage of the Euler angles representation is its straightforward physical interpretation, of importance for human-machine interfaces. The disadvantage lies in fact that the construction of the attitude matrix requires the evaluation of trigonometric functions, of higher computational load than other parameterizations. Also, the derivatives of the components of the rotation matrix are nonlinear (trigonometric), and affected by singularities. The three consecutive rotations that rotate the frame {u 1 , u 2 , u 3 } into the frame {u ′′′ 1 , u ′′′ 2 , u ′′′ 3 }.Thefirstoneisaboutthemainaxisu 3 and magnitude ψ, the second is about the main axis u ′ 2 and magnitude θ and the third is about the main axis u ′′ 1 and magnitude φ.

Quaternions
A quaternion is an order-4 vector whose components can be used to define the mutual rotation between reference systems:q = (q 1 , q 2 , q 3 , q 4 ) T q 4 is named the scalar (real) component of the quaternion, whereas (q 1 , q 2 , q 3 ) T forms the imaginary (or vectorial) part. The components of a quaternion must respect the constraint q Tq = 1. Physically, the four components ofq define the magnitude and axis of the rotation necessary to rotate one reference system into the other, see Figure 4. The attitude matrix R is parameterized in terms of quaternions as This parameterization is non ambiguous and it does not involve any trigonometric function, so that the computational burden is lower than with other representations. The quaternion representation is of common use in attitude estimation and control applications, since it guarantees high numerical robustness.

Attitude estimation
As it will be shown in the next sections, the least-squares solution of model (16), requires the solution of a constrained least-squares problem of the type: with · 2 Q = (·) T Q −1 (·).T h es h a p eo fQ drives the choice of the solution technique to be adopted for solving (25).
If Q is a diagonal matrix, problem (25) becomes an Orthogonal Procustes Problem (OPP), see Schonemann (1966). This class of constrained least-squares problem have been thoroughly analyzed, and fast algorithms have been devised to quickly extract the minimizerŘ,s e e (Davenport, 1968;Shuster & Oh, 1981). Various fast methods for the solution of an OPP have been introduced -and widely used in practice -based on the Singular Value Decomposition (SVD) or the EIGenvalues decomposition (EIG), such as the QUaternion ESTimator (QUEST) (Shuster, 1978;Shuster & Oh, 1981), the Fast Optimal Attitude Matrix (FOAM) (Markley & Landis, 1993), the EStimator of the Optimal Quaternion (ESOQ) (Mortari, 1997) or the Second ESOQ (ESOQ2) (Mortari, 2000) algorithms, which have been extensively compared in Markley & Mortari (1999; and Cheng & Shuster (2007).
For nondiagonal matrices Q, the extraction of the orthonormal attitude matrixŘ has to be performed through nonlinear estimation techniques. A first numerical scheme for the solution of (25) is derived by applying the Lagrangian multipliers method. The Lagrangian function is with [λ] q the q by q matrix of Lagrangian multipliers: The last term of (26) gives the q(q+1) 2 constraining functions that follows from the orthonormality of R: q constraints are given by the normality (unit length) of the columns of R,whereas q(q−1) 2 constraints are given by the orthogonality of the columns of R.
The gradient of the Lagrangian function (26), together with the q(q+1) 2 constraining functions, defines the nonlinear system to be solved: Due to the symmetry of matrix R T R − I q , only its upper (or lower) triangular part has to be considered in (28). The Newton-Raphson method can then be applied to iteratively converge to the sought orthonormal matrix of rotations.
This method is computationally heavier than other iterative schemes, since it requires the explicit computation of larger-sized matrices than other methods given in the following.
A second viable solution scheme is obtained by re-parameterizing the attitude matrix with the vector of Euler angles μ = (ψ, θ, φ) T . Following the reparameterization, matrix R (μ) implicitly fulfills the constraint R T R = I q , and problem (25) is rewritten aš with h(μ)=Q − 1 2 vec R − R(μ) . The nonlinear least-squares problem (29) is solved by applying iterative methods, e.g., the Newton method. This approach (Euler angles parameterization) works with a minimal set of unknowns -the Euler angles -and it can quickly converge to the sought minimizer if an accurate initial guess is used. The disadvantage is that trigonometric functions have to be evaluated, increasing the computational load.
A third viable approach is devised by employing the quaternions parameterization of R and to solve for (25):q = arg min The orthonormality of R is guaranteed by the normality of the quaternion: this introduces a single constraint in the minimization problem (30). A Lagrangian function is formed as and the (nonlinear) system to be solved is with J R(q) the Jacobian of vec(R(q)).
The three iterative solutions given above rigorously solve for problem (25), but are generally slower than the methods available for diagonal Q matrices (SVD, EIG, QUEST, FOAM, ESOQ, and ESOQ2). Figure 5a illustrates the mean number of floating-point operations for different attitude estimation methods, per number of baselines employed. 10 4 samplesR have been generated via Monte Carlo simulations for a given fully-populated Q matrix. The gray bars span between the maximum and minimum numbers obtained for each algorithm. The off-diagonal elements of Q a r ed i s r e g a r d e dw h e na p p l y i n gt h eS V D ,E I G ,Q U E S T ,F O A M , ESOQ, and ESOQ2 methods. These techniques outperform each iterative method: the number of required floating-point operations is generally two to three orders of magnitude lower. Among the iterative methods, the Lagrangian multiplier technique generally requires the highest number of operations, making it the least efficient method, while the Euler angle method and the Quaternion parameterization provide better overall results. Figure 5b shows the corresponding mean, maximum and minimum computational times marked during the simulations. The Lagrangian parameterization method generally takes the longest time to converge, whereas the quaternion and Euler angle methods show better results. Note that higher number of floating operations does not directly translate into longer computational times, because modern processor architectures efficiently operate by means of multi-threading and parallel processing.

Reliable attitude-ambiguity estimation methods
This section reviews the solution of the GNSS attitude model (16). This can be presented by addressing two consecutive steps: float estimation and ambiguity resolution.

Float ambiguity-attitude solution
We indicate with float the solution of (16) obtained by disregarding the whole set of constraints, i.e., the integerness of Z and the orthonormality of R:   This float solution follows from solving the system of normal equations where the v-c matrix Q YY is written as in (11). Inversion of the normal matrix M gives the v-c matrix of the float estimatorsR andẐ: The float estimators are explicitly derived aŝ Next to the above float solution, we can also define the following conditional float solution for the attitude matrix: In this case the ambiguity matrix is assumed completely known. The solutionR(Z) can be computed form the float solutionsR andẐ as: Application of the variance propagation law gives There is a very large difference in the precision of the float solutionR and the precision of the conditional float solutionR(Z). This can be demonstrated by comparing expression (39) with QˆRˆR in (35), whose relation with the design matrices can be made explicit as is characterized by much smaller entries than G T Q −1 yy G −1

.T h i s i s demonstrated as follows. Matrices A, G and Q yy may be partitioned as
where we assumed, for simplicity, the same code and phase standard deviations for each observation, independent from the combination of satellites, receivers and frequency. Λ is the diagonal matrix of carrier wavelengths, whereas Q is the matrix that introduces correlation due to the DD operation.
It follows that

206
Recent Advances in Aircraft Technology

www.intechopen.com
The ratio between the entries of matrix QRR and QR (Z)R(Z) is then proportional to the ratio In GNSS applications, this phase-code variance ratio is in the order of 10 −4 . This clearly demonstrates the importance of ambiguity resolution: if we can integer-estimate Z with sufficiently high probability, then the attitude matrix R can be estimated with a precision that is comparable with the high precision ofR(Z).

Ambiguity resolution
The second step consists of the resolution of the carrier phase integer ambiguities. The solution of model (16) is obtained through the following C-ILS minimization problem: Both sets of constraints are now imposed: the matrix of ambiguitiesŽ is integer valued and the matrixŘ belongs to the class of 3 × q orthonormal matrices O 3×q .T h eC -I L Ss o l u t i o nŽ can be computed from the float solutions as (Teunissen & Kleusberg, 1998): The cost function C(Z) is the sum of two terms. The first weighs the distance between a candidate integer matrix Z and the float solutionẐ, weighted by the v-c matrix QẐẐ.T h e second weighs the distance between the conditional (on the candidate Z) attitude matrix R(Z) and the orthonormal matrixŘ(Z) that follows from the solution of (45). Therefore, the computation of cost function C(Z) also involves a term that weighs the distance of the conditional attitude matrix from its orthogonal projection. This second term greatly aids the search for the correct ambiguities: integer candidates Z that produce matricesR(Z) too far from their orthonormal projection contribute to a much higher value of the cost function.
Since the minimization problem (44) is not solvable analytically due to the integer nature of the parameter involved, an extensive search in a subset of the space of integer matrices Z mN×n has to be performed. The definition of an efficient and fast solution scheme for problem (44) is not a trivial task. In order to highlight the intricacies of such formulation, we first give an approximate solution, obtained by neglecting the orthonormal constraint.

The LAMBDA method
Consider first the integer minimization problem (44) without the orthonormality constraint on R. Then the second term of C(Z) reduces to zero and the integer minimization problem becomesŽ This is the usual approach of doing GNSS integer ambiguity resolution. Due to the absence of the orthonormality constraint on R one may expect lower success rates, i.e., lower probability of identifying the correct ambiguity matrix Z. However, the ILS problem (46) is of lower complexity than (44) and a very fast implementation of it is available: the LAMBDA (Least-squares AMBiguity Decorrelation Adjustment) (Teunissen, 1995) method, see, e.g., Boon & Ambrosius (1997); Cox & Brading (2000); Huang et al. (2009);Ji et al. (2007); Kroes et al. (2005). It consists of two steps, namely decorrelation and search.
The integer minimizer has to be extensively searched within a subset of the whole space of integers: Ω U is the so-called search space, a region of the space of integer matrices that contains only those candidates Z for which the squared norm (46)  , is an option, as well as bootstrapping an integer matrix, as in Teunissen (2000;2007b).
Searching for the integer minimizer in Ω U proves inefficient due to the weight matrix QẐẐ. Geometrically, the search space defines a hyperellipsoid centered inẐ and whose shape and orientation are driven by the entries of matrix QẐẐ. The difficulty of the search lies in the fact that the search space is highly elongated, as detailed in Teunissen & Kleusberg (1998). The reason is that the ambiguities are highly correlated. While the set wherein the independent ambiguities (e.g., three ambiguities for a single baseline scenario) can be chosen is rather large, the set of admissible values for the remaining ambiguities is very small. This causes major halting problems during the search, since many times the selected subset of independent ambiguities does not yield admissible integer matrix candidates. This issue is tackled and solved in the LAMBDA method with a decorrelation step. The decorrelation of matrix QẐẐ is achieved by an admissible transformation matrix T. In order to preserve the integerness, such matrix has to fulfill the following two conditions: T as well as its inverse T −1 need to have integer entries. The matrix of transformed ambiguities Z ′ and corresponding v-c matrix are then obtained as The decorrelation procedure is described in Teunissen & Kleusberg (1998). The v-c matrix is iteratively decorrelated by a sequence of admissible transformations T i , until matrix cannot be further decorrelated. Note that due to the integer conditions on T,af u l l decorrelation cannot generally be achieved. Figure 6 shows three steps of the decorrelation process for a two-dimensional example. Figure 6a shows the original (elongated) ellipse associated to QẐẐ, Figure 6b shows an intermediate decorrelation step, and Figure 6c shows the final decorrelated search space.
After the decorrelation step, the actual search is performed by operating the LDL T factorization of matrix QẐ ′Ẑ′ , so that the quadratic form in (46) can be written as a summation: where the scalarsẑ ′ i|I and σ 2 i|I are the conditional float ambiguity estimator and its corresponding conditional variance, respectively. These are conditioned to the previous I = 1,...,i − 1 values, and directly follow from the entries of matrices L and D. More details on the way the search is actually performed can be found in de Jonge & Tiberius (1996).
Due to the decorrelation step, the extensive search for the integer minimizerŽ U is performed quickly and efficiently, making the LAMBDA method perfectly suitable for real-time applications.

The MC-LAMBDA method
The MC-LAMBDA method is an extension of the LAMBDA method that applies to the geometrically-constrained problem (44). The MC-LAMBDA method shares the same working principle of the LAMBDA method: first the search space is decorrelated, then the search for the integer minimizer is performed. However, an extensive search within a (decorrelated) search space is generally not efficient as it is with the LAMBDA method, as explained in the following.
The search space is now defined as The cost function C(Z) takes, for the same candidate Z, much larger values than the first quadratic term in (44), due to the matrix QR (Z)R(Z) , whose inverse has entries two orders of magnitude larger than the entries of Q −1 ZẐ Teunissen, 2007a). For this reason it is not trivial to set a proper value of χ 2 , since the cost function C(Z) is highly sensitive to the choice of Z (Giorgi, 2011;Giorgi et al., 2011). This problem becomes more marked for weaker models (single frequency, low number of satellites tracked, high noise levels). Obviously, larger values of χ 2 imply longer computational times due to the larger number of candidates to be evaluated. Also, the constrained least-squares problem (45) has to be solved for each of the integer candidates in Ω C χ 2 , thereby further increasing the computational load.
The aforementioned issues are solved with a novel numerical efficient search scheme for the solution of (44). This is achieved by employing easier-to-evaluate bounding functions and introducing new search algorithms.
First, consider two functions, C 1 (Z) and C 2 (Z), that satisfy the following inequalities: These functions provide a lower and an upper bound for the cost function C(Z). The choice for these bounding functions is driven by two requirements: their evaluation should be less time consuming than the evaluation of C(Z), and each bound should be sufficiently tight. Several alternatives have been studied in (Giorgi, 2011;Giorgi et al., 2012;Nadarajah et al., 2011;Teunissen, 2007a;c), based on -the eigenvalues of matrix Q −1

R(Z)R(Z)
-the analytical solution of Wahba's problem (Wahba, 1965) -a tighter geometrical bound based on Procustes problem (Schonemann, 1966) -a QR factorization (Gram-Schmidt process) For example, the first method listed exploits the inequalities ξ m · 2 I ≤ · 2 Q ≤ ξ M · 2 I , with ξ m and ξ M the smallest and largest eigenvalues of Q −1 R(Z)R(Z) , respectively. After some manipulation, the two bounding functions read wherer i (Z) are the column vectors ofR(Z).
Two efficient search methods have been developed to reduce the computational burden associated to an extensive search. Independently from the bounding functions used, these novel search schemes allow for a quick minimization of C(Z).
Consider first the lower bound C 1 (Z). The search space associated to C 1 (Z) is Obviously, the search space Ω C χ 2 is contained within Ω 1 χ 2 . One may proceed, for example, by choosing χ 2 = vec(Ẑ − Z ′ ) 2 QẐẐ with Z ′ a given integer matrix (both rounding the float solution and bootstrapping an integer matrix are viable choices). Then, we can enumerate all the integers matrices contained in Ω 1 χ 2 and compute C(Z) for each candidate (if any, since set Ω 1 χ 2 may also turn out empty), in order to also evaluate Ω C χ 2 . If this set turns out non-empty, then one has simply to extract the minimizerŽ by sorting the integer matrices according to the values of C(Z). However, there is no guarantee that Ω C χ 2 is non-empty. If the search space Ω C χ 2 is empty, the size of Ω 1 χ 2 is increased and the process repeated iteratively until the minimizerŽ is found. This search scheme, illustrated with the flow chart in Figure 7, is named Expansion approach, since the size of the search space is iteratively 'expanded'.
An alternative approach is devised by considering the upper bound C 2 (Z). Its search space is which is contained in the set Ω C χ 2 . Consider the following iterative procedure. First, the scalar χ 2 is set such that it guarantees the non-emptiness of Ω 2 χ 2 , and therefore Ω C χ 2 is non-empty either. This can be done by choosing χ 2 = C 2 (Z ′ ) for an integer matrix Z ′ ,which can be the rounded float solution, a bootstrapped solution, or an integer matrix obtained by other means (see for further options Giorgi et al. (2008)). Then, the search proceeds by looking for an integer candidate in the set Ω 2 χ 2 , aiming to find a matrix Z 1 that provides a smaller valuefortheupperboundC 2 (Z 1 )=χ 2 1 < χ 2 . When it is found, the set is shrunk to Ω 2 χ 2 1 and the search continues by looking for another integer candidate Z 2 capable of reducing the value C 2 (Z 2 )=χ 2 2 < χ 2 1 . This process is repeated until the minimizer of C 2 (Z),sa yŽ 2 , is found. Since this may differ from the minimizer of C(Z), the search space Ω C χ 2 ,w i t h χ 2 = C 2 (Ž 2 ), is evaluated and the sought-for integer minimizerŽ extracted. This iterative search scheme is named Search and Shrink approach, and it is detailed in the flow chart of Figure 8.
Both the Expansion and the Search and Shrink approaches implement the search for integer minimizer (44) in a fast and efficient way, such that the algorithm can be used for real-time applications.
The MC-LAMBDA method achieves very high success rates. The success rate is defined as the probability of providing the correct set of integer ambiguities. The inclusion of geometrical constraints, which follow from the a priori knowledge of the antennas relative positions aboard the aircraft, largely aids the ambiguity resolution process, allowing for higher success rates in weaker models, such as with the single-frequency and/or high measurement noise scenarios. These performance improvements associated to the MC-LAMBDA method with respect to classical methods (such as the LAMBDA) are analyzed in the following section with actual data collected during two different flights tests.

Flight test results
The performance of the MC-LAMBDA method is analyzed with data collected on two flight-tests performed with a Cessna Citation jet aircraft. The aircraft attitude is extracted from unaided, single-epoch, single-frequency (N = 1) GNSS observations, in order to demonstrate the method capabilities in the most challenging scenario, i.e., stand-alone, high observation noise and low measurements redundancy. Also, single-epoch performance is extremely important for dynamic platforms, where a quick recovery from changes of tracked satellites, cycle slips and losses of lock is necessary to avoid undesired loss of guidance. The  single-frequency case is of interest for many aerospace applications, where limits on weight and power consumption must often be respected.
In both tests the same receiver (Septentrio PolaRx2@) was connected to three antennas, placed on the middle of the fuselage, on the wing and on the nose (see Figure 9). In the first  Fig. 8. The Search and Shrink approach: flow chart.
test analyzed (T − I) the nose antenna was placed on the extremity of a boom, whereas in second test (T − II) it was directly placed on the aircraft body. The two tests largely differ by the flight dynamic. Test T − I was conducted with aggressive maneuvering and few zero-gravity parabolas, whereas T − II was performed as part of a gravimetry campaign, with very few smooth maneuvers, as shown in Figure 10. During test T − II,t h ea i r c r a f t was also equipped with an Inertial Navigation System (INS), whose output is used to test the GNSS-based attitude estimation accuracy. Figure 11 reports the number of tracked satellites for the duration of the two tests. The PDOP (Precision Dilution of Precision) is also shown.

Instantaneous ambiguity resolution
The success rate marked by the LAMBDA and MC-LAMBDA methods applied to both flight tests is reported in Instead, application of the MC-LAMBDA method yields a strong performance improvement. The constrained method is capable of providing the correct integer solution for more than 80% of the epochs in test T − I and more than 88% in test T − II.
Both airborne tests confirm the very large improvement that is obtained by strengthening the underlying model with the inclusion of geometrical constraints. It is stressed that all the ambiguity resolution performance reported are obtained by processing the GNSS signals without any a priori information or assumption about the attitude or the aircraft motion. Also, mask angles, elevation-dependent models, dynamic models or any kind of filtering are not applied.  Figure 12. The high dynamics of the flight is evident from the steep variations of the attitude. In particular, Figure 13 shows a zero-gravity maneuver: the aircraft promptly pitched up, gained some altitude, and performed an ample arc to create a virtual absence of gravity on board. Figure 14 shows the GNSS-based attitude angles for the test T − II. The INS solutions are also reported in the figures, in order to provide a comparison between the two systems. Table 2 reports the standard deviations of the differences between the INS and GNSS-based attitude estimations. Taking the precise INS output as benchmark solution, it can be inferred that the accuracy obtained is within the expected range, given the baseline lengths employed. The heading angle is estimated with the highest precision, whereas the elevation estimation is characterized by the highest noise levels. This is due to the relative geometry of the antennas and to the fact that the vertical components of the GNSS-based baseline estimations are inherently less accurate that the horizontal components. The bank angle is estimated with higher precision than the elevation angle, being driven by the longer baseline Body − Wing.   (b) Altitude profile during the zero-gravity maneuver. Fig. 13. T − I test: zero-gravity maneuver.

Summary and conclusions
Ambiguity resolution can be effectively enhanced by means of a rigorous formulation of the ambiguity-attitude estimation problem. In order to infer the aircraft's orientation from the GNSS antenna positions, each antenna location on the aircraft body has to be precisely known. This geometrical information can be embedded in the ambiguity resolution step, thus strengthening the underlying functional model -i.e., additional information is added to the functional model -and enhancing the whole estimation process. The higher ambiguity resolution performance comes at the cost of an increased computational complexity. In order to overcome the issue, a number of solutions are presented, which allow for fast and reliable solutions without requiring extensive computational loads. A fast implementation of the geometrically constrained problem is obtained by modifying a well-known method for ambiguity resolution: the LAMBDA (Least-squares AMBiguity Decorrelation Adjustment) method. This method is nowadays the standard for carrier-phase based applications, and it is being implemented in a number of receivers employed for high-precision navigation applications. The complexity of the constrained estimation method requires the development of novel strategies to extract the solution in a timely manner. This is achieved by properly modifying the LAMBDA method to address the specific ambiguity-attitude estimation problem: the Multivariate Constrained (MC)-LAMBDA method. Through the use of two novel search schemes the sought-for set of carrier phase ambiguities can be efficiently estimated.
The method is tested on actual data collected on two different flight tests. Each test indicates the feasibility of employing GNSS as attitude sensor, an application that might be increasingly adopted in the aviation industry, either stand-alone for non-critical applications, or in combinations with other sensors for safety-critical applications.


extfieldsofstudy:
Engineering
Computer Science

text:
Development of an EKF-based localization algorithm using compass sensor and LRF

This paper presents the implementation of a perceptual system for a mobile robot. The system is designed and installed with modern sensors and multi-point communication channels. The goal is to equip the robot with a high level of perception to support a wide range of navigating applications including Internet-based telecontrol, semi-autonomy, and autonomy. Due to uncertainties of acquiring data, a sensor fusion model is developed, in which heterogeneous measured data including odometry, compass heading and laser range is combined to get an optimal estimate in a statistical sense. The combination is carried out by an extended Kalman filter. Experimental results indicate that based on the system, the robot localization is enhanced and sufficient for navigation tasks.

I. INTRODUCTION
Perceptual system is the key to the intelligence of a mobile robot. When a mobile robot travels in an unknown or partially known environment, it must understand the environment for obstacle avoidance or path planning. The perception of robot is done by taking measurements using various sensors and then extracting meaningful information from those measurements.
Based on advanced material technologies, modern sensors can be nowadays equipped for the mobile robot such as optical incremental encoders, heading sensors, ultrasonic sensors, infrared sensors, laser range finders and vision systems. These sensors are selected so as to accord with the goal of application, the specific constraints of the working environment, and the individual properties of the sensors themselves. Nevertheless, there is no single sensor which can adequately capture all relevant features of a real environment. It is necessary to combine the data from different sensors into a process known as sensor fusion. The expectation is that the fused data is more informative and synthetic than the original.
Several methods have been reported to cope with this trend. Durrant-Whyte has developed a multi-Bayesian estimation technique for combining touch and stereo sensing [1], [2]. Tang and Lee proposed a generic framework that employed a sensorindependent, feature-based relational model to represent information acquired by various sensors [3]. In [4], a Kalman filter update equation was developed to obtain the correspondence of a line segment to a model, and this correspondence was then used to correct position estimation. In [5], an extended Kalman filter was conducted to manipulate image and spatial uncertainties.
In this work, we develop a multi-sensor perceptual system for the mobile robot. Sensors include but not limit to optical quadrature encoders, compass sensors, ultrasonic sensors, laser range finders, global positioning systems (GPS) and vision systems. The goal is to equip the robot with diverse levels of perception to support a wide range of navigating applications including Internet-based telecontrol, semi-autonomy, and autonomy. At this stage of research, we use the optical quadrature encoders for position measurement, compass sensor for deflect angle calculation and laser range finder (LRF) for object-boundary detection. These types of data have their strengths and are often used in robotic applications, but combining them provides even more useful information. This combined approach deserves further investigation. In our system, an extended Kalman filter (EKF) has been designed. It fuses the raw measurement data of optical quadrature encoders, the compass sensor, and the LRF to obtain optimal estimation of robot positions while reduces uncertainties in measurements. The proposed fusing algorithm also provides constraints that filter out unreliable data from the sensor readings. Outputs of EKF combined with the boundaries of objects detected from LRF ensure the success navigation of the mobile robot in indoor environments.
The paper is arranged as follows. Details of the perceptual system are described in Section II. The algorithm for sensor fusion using EKF is explained in Section III. Section IV introduces experimental results. The paper concludes with an evaluation of the system, with respect to its strengths and weaknesses, and with suggestions of possible future developments.

II. SENSOR SYSTEM DESIGN
The sensor system consists of an Omni-directional camera, a GPS, a compass sensor, a LRF, three quadrature encoders and eight ultrasonic sensors. Fig.1 describes the sensors in relation with actuators and communication channels in a mobile robot. The communication is performed via low-rate and high-rate channels. The low-rate channel with standards of RS-485 is developed by an on-board 60MHz Microchip dsPIC30F4011based micro-controller and MODBUS protocol for multi-point interface. The high-rate channels use the USB-to-COM and IEEE-1394/firewire available commercial ports. In the system, optical quadrature encoders are used. An optical encoder is basically a mechanical light chopper that produces a certain number of sine or square wave pulses for each shaft revolution. As the diameter of wheel and the gear coefficient are known, the angular position and speed of wheel can be calculated. In the quadrature encoder, a second light chopper is placed 90 degrees shifted with respect to the original resulting in twin square waves as shown in fig.2. Observed phase relationship between waves is employed to determine the direction of the rotation. In the system, measurements from encoders are used as input data for positioning and feedback for a closed-loop motor speed controller.  The GPS is mainly applied for positioning in the outdoor environment. A HOLUX GPS UB-93 module is used [7]. Due to the presence of networking in our system, an Assisted GPS (A-GPS) can be also used in order to locate and utilize satellite information of the network in the poor signal condition.
The system provides eight SFR-05 ultrasonic sensors split into four arrays, two on each, arranged at four sides of the robot. The measuring range is from 0.04m to 4m.
The vision system is detachable and mounted on the top of the robot. It mainly consists of an Omni-directional digital camera Hyper-Omni Vision SOIOS 55 with a high-rate IEEE-1394 data transfer line. With a 360-degree field of view, the Omni-directional camera is a good vision sensor for landmarkbased navigation [8].
Last but not least, a 3D laser range finder (LRF) with a range from 0.04m to 80m is developed for the system. Its operation is based on the time-of-flight measurement principle. A single laser pulse is emitted out and reflected by an object surface within the range of the sensor. The lapsed time between emission and reception of the laser pulse is used to calculate the distance between the object and the sensor. By an integrated rotating mirror, the laser pulses sweep a radial range in its front so that a 2D field/area is defined ( fig.4). Due to the fixation of the pitching angle in the scanning plane, the information of 2D image may be, in certain cases, insufficient for obstacle detection ( fig.4). In those situations, a 3D image is necessary. As the 2D scanner is popular and lowcost, building a 3D LRF from the 2D is usually a prior choice [9][10]. In our system, a 3D LRF is developed based on the 2D SICK-LMS 221 [9]. It has the view angle of 180° in horizontal and 25° in vertical. During the measurement time, a set of deflect angle β and distance ρ values are received. The set (β, ρ) is then combined with the pitching angle α to create the 3D data. Based on these data, we can define the Cartesian coordinates of one point in the space.   The proposed sensor platform equips the robot with the ability to perceive many parameters of the environment. Their combination, however, presents even more useful information. In this work, the raw data of three different types of sensors including the compass sensor, the quadrature encoder and the LRF is syndicated inside an EKF. The aim is to determine the robot position during operation as accurately as possible. We start with the kinematic model of the mobile robot. The two wheeled, differential-drive mobile robot with non-slipping and pure rolling is considered. Fig.7 shows the coordinate systems and notations for the robot, where (X G , Y G ) is the global coordinate, (X R , Y R ) is the local coordinate relative to the robot chassis. R denotes the radius of driven wheels, and L denotes the distance between the wheels.
During one sampling period ∆t, the rotational speed of the left and right wheels ω L and ω R create corresponding increment distances ∆s L and ∆s R traveled by the left and right wheels respectively: These can be translated to the linear incremental displacement of the robot's center ∆s and the robot's orientation angle ∆θ : The coordinates of the robot at time k+1 in the global coordinate frame can be then updated by: In practice, (3) is not really accurate due to unavoidable errors appeared in the system. Errors can be both systematic such as the imperfectness of robot model and nonsystematic such as the slip of wheels. These errors have accumulative characteristic so that they can break the stability of the system if appropriate compensation is not considered. In our system, the compensation is carried out by the EKF.

Let [
] T x yθ = x be the state vector. This state can be observed by some absolute measurements, z. These measurements are described by a nonlinear function, h, of the robot coordinates and an independent Gaussian noise process, v. Denoting the function (3) is f, with an input vector u, the robot can be described by: where the random variables w k and v k represent the process and measurement noise respectively. They are assumed to be independent to each other, white, and with normal probability distributions: The steps to calculate the EKF are then realized as below: 1. Prediction step with time update equations: where ˆn k − ∈ ℜ x is the priori state estimate at step k given knowledge of the process prior to step k, ˆk − P denotes the covariance matrix of the state-prediction error, A k is the Jacobian matrix of partial derivates of f to x: W is the Jacobian matrix of partial derivates of f to w: and 1 k − Q is the input-noise covariance matrix which depends on the standard deviations of noise of the right-wheel rotational speed and the left-wheel rotational speed. They are modeled as being proportional to the rotational speed ω R,k and ω L,k of these wheels at step k. This results in the variances equal to 2 R δω and 2 L δω , where δ is a constant determined by experiments. The input-noise covariance matrix Q k is defined as: 2. Correction step with measurement update equations: where ˆn k ∈ ℜ x is the posteriori state estimate at step k given measurement k z , K k is the Kalman gain, H is the Jacobian matrix of partial derivates of h to x: R k is the covariance matrix of noises estimated from the measurements of compass sensor and LRF as follows. x G cos β j + y G sin β j = ρ j When the robot moves, a new scan of LRF is performed and a new map of environment, namely local map, is constructed which also consists of a set of line segments described by the equation: where ψ i and r i are the parameters of lines ( fig.8).
The line segments of the global and local map are then matched using weighted line fitting algorithm [11]. The matching line parameters i r and i ψ from the local map are collected in the vector z k , which is used as the input for the correction step of the EKF to update the robot's state: From (16) and (19), the data of the compass sensor is utilized to correct the robot orientation. At step k, this data is employed as the absolute measurement for the element θ k of z k .
The noise of this measurement is achieved from the manufactory specification. The accuracy of 0.1 0 corresponds to the noise variance of 0.01. This collected into the covariance matrix R i,k (19) gives R k for the correction step of EKF.

IV. EXPERIMENTS
The setup of the sensor system has been described in Section II. In this section, experiments are conducted to evaluate the efficiency of the fusion algorithm.

A. Experimental Setup
A rectangular shaped flat-wall environment constructed from several wooden plates surrounded by a cement wall is setup. The robot is a two wheeled, differential-drive mobile robot. Its wheel diameter is 10 cm and the distance between two drive wheels is 60cm.
The speed stability of the motor during each sampling time step is an important factor to maintain the efficiency of the EKF. For that reason, motors are controlled by microprocessorbased electronic circuits which carry out a PID algorithm. The stability checked by a measuring program written in LABVIEW is ±5%.
The compass sensor has the accuracy of 0.1 0 . The LRF has the accuracy of 30mm in distance and 0.25 0 in deflect angle.
The sampling time ∆t of the EKF is 100ms. The proportional factor δ of the input-noise covariance matrix Q k is experimentally estimated as follows. The deviation between the true and the position estimated by the kinematic model when driving the robot straight forwards several times (from the minimum to the maximum tangential speed of the robot) is determined. The deviation between the true robot orientation and the orientation estimated by the kinematic model when rotating the robot around its own axis several times (from the minimum to the maximum angular speed of the robot) is also determined. Based on the error in position and orientation, the parameter δ is calculated. In our system, δ is estimated to be the value 0.01.

B. Sensor Fusion Algorithm Evaluation
In order to evaluate the efficiency of the fusion algorithm, different configurations of the EKF were implemented. Fig.9 describes the trajectories of a robot movement estimated by four methods: the odometry, the EKF with compass sensor, the EKF with LRF, and the EKF with the combination of LRF and compass sensor. The deviations between each trajectory and the real one are represented in fig.10.

a)
In another experiment, the robot is programmed to follow predefined paths under two different scenarios: with and without the EKF. Fig.11a represents the trajectories of the robot moving along a rounded rectangular path in which the one with dots corresponds to the non-existence of EKF and the one with pluses corresponds to the existence of EKF in the implementation. The trajectories in case the robot follows a rounded triangular path are shown in fig.11b.
It is concluded that the EKF algorithm improves the robot localization and its combination with LRF and compass sensor gives the optimal result.

C. Autonomous Navigation
In this experiment, we evaluate the applicability of the proposed fusion approach in a real autonomous navigation application. The goal is to navigate the mobile robot to go a round closed to the line segments of the global map extracted from the first scan of the LRF. Fig.12 shows the extracted map and the success trajectory of the robot during the navigation. A sequence of images showing the motion of the robot in the experiment is described in fig.13. It is necessary to develop a perceptual system for the mobile robot. The system is required to not only be wellworking but also critically support various levels of perception. In this work, many types of sensors including position speed encoders, integrated sonar ranging sensors, compass and laser finder sensors, the global positioning system (GPS) and the visual system have been implemented in a real mobile robot platform. An EKF has been designed to fuse the raw data of compass sensor and LRF. Experiments show that this novel combination significantly improves the accuracy of robot localization and is sufficient to ensure the success of robot navigation. Further investigation will be continued with more sensor combination to better support the localization in outdoor environments.
ACKNOWLEDGMENT This work was supported by Vietnam National Foundation for Science and Technology Development (NAFOSTED).