Abstract
Aircraft and spacecraft navigation precision is dependent on the measurement system for position and attitude determination. Rotation of an aircraft can be determined measuring two vectors in two different reference systems. Velocity vector can be determined in the inertial reference frame from a GNSS-based sensor and by integrating the acceleration measurements in the body reference frame. Estimating gravity vector in both reference frames, and combining with velocity vector, determines rotation of the body. A new approach for gravity vector estimations is presented and employed in an attitude determination algorithm. Nonlinear simulations demonstrate that using directly the positioning and velocity outputs of GNSS sensors and strap-down accelerometers, aircraft attitude determination is precise, especially in ballistic projectiles, to substitute precise attitude determination devices, usually expensive and forced to bear high solicitations as for instance G forces.
1. Introduction
Obtaining precise attitude information is essential for navigation and control. The effectiveness of navigation and control is determined by the degree of precision of navigation and control systems, including inertial measurement units [1]. There is an extensive body of literature regarding attitude estimation using various sensor inputs [2].
Traditionally, in order to obtain accurate values for determining attitude, expensive and/or weighty units, such as laser or fiber optic gyroscopes and accelerometers, or their MEMS equivalent must be employed. Moreover, when high-demanding maneuvers are performed, this equipment may become extremely expensive.
Other methods, based on very precise and advanced carrier-phase single-frequency GNSS attitude determination algorithms or multi-GNSS developments of integrating GNSS systems for attitude determination, are not possible to be applied on these artillery devices due the small space available on the vehicle. At present, the best carrier-phase-based method for GNSS attitude determination is the multivariate constrained MC-LAMBDA method. This method is currently the best performing attitude determination method as was demonstrated in land, sea, and air experiments. Example publications describing such actual real-world GNSS-attitude determination results can be found in [3–6]. The underlying theory and advanced algorithms for much of the above-referred integer-ambiguity-resolved carrier-phase GNSS attitude determination methods can be found in [7–9]. Multi-GNSS developments of integrating GNSS systems for attitude determination are stated in [10, 11].
It is well known that the attitude of an aero-vehicle can be determined, starting from an initial condition, by integrating the angular rates (pitch, roll, and yaw rates) of the vehicle and propagating them forward time. Nevertheless, accuracy requirements usually cannot be satisfied by using inexpensive sensors [1]. This problem becomes even more important when the vehicle cannot be reused: low-cost attitude determination systems are of key importance for these applications. For example, Gebre-Egziabher et al. [12] describe an attitude determination system that is based on two measurements of nonzero, noncollinear vectors. Using the Earth’s magnetic field and gravity as the two measured quantities, a low-cost attitude determination system is proposed. Gebre-Egziabher et al. [13] develop an inexpensive attitude heading reference system for general aviation applications by fusing low-cost automotive-grade inertial sensors with GPS. The inertial sensor suit consists of three orthogonally mounted solid-state-rate gyros. Eure et al. [14] describe an attitude estimation algorithm derived by postprocessing data from a small low-cost inertial navigation system recorded during the flight of a subscale commercial off-the-shelf UAV. Estimates of the UAV attitude are based on MEMS gyro, magnetometer, accelerometer, and pitot tube inputs.
Henkel and Iafrancesco [15] state that low-cost GNSS receivers and antennas can provide a precise attitude and drift-free position information, but accuracy is not continuous. Inertial sensors are robust to GNSS signal interruption and very precise over short time frames, which enables a reliable cycle slip correction, but low-cost inertial sensors suffer from a substantial drift. The authors propose a tightly coupled position and attitude determination method for two low-cost GNSS receivers, a gyroscope and an accelerometer, and obtain a heading with an accuracy of 0.25°/baseline length (m) and an absolute position with an accuracy of 1 m. Similar developments may be found within space vehicles, for example, in [16]. In [17], the use of an inertial navigation system (INS) and a multiple GPS antenna system for attitude determination of an off-road vehicle is developed, and in [18], attitude determination using the GPS carrier phase has been applied successfully to aircrafts in experiments by a number of researchers.
Also, improved algorithms for estimating attitude in case of failures have been proposed in the literature. For example, Soken and Hajiyev [19] introduce algorithms with filter gain correction for the case of measurement malfunctions. Two different algorithms are proposed and applied for the attitude estimation process of a pico satellite. The results of these algorithms are compared for different types of measurement faults in different estimation scenarios, and recommendations about their applications are given. Another example is exposed in [20]. This algorithm is based on interacting multiple models (IMM) extended Kalman filters (EKF) using a star sensor and gyroscope data. In Liu et al. 2011, a computationally efficient nonlinear attitude estimation strategy is based on the vector observations for satellites.
However, as stated in [21], many of the presented methods such as those employing local magnetic field vectors are only valid for estimating the orientation of a slow-rotation body. On the other hand, for high spin rate bodies, electromagnetic interactions degrade magnetic measurements. In the research described in this paper, two measured quantities are used to obtain attitude information for high dynamic and high spin rate vehicles: speed and gravity vectors. They are obtained in two different reference frames using a GNSS sensor and a strap-down accelerometer.
1.1. Contributions
The main contributions of this paper are twofold. First is the development of a novel algorithm which avoids the use of gyroscopes for attitude determination. This is purposed to decrease costs in attitude sensors and even to improve attitude determination by applying filtering techniques, especially for artillery device applications, where high acceleration conditions increase the price of precise attitude determination devices such as gyroscopes. Second is the development of an estimation method in order to obtain the gravity vector in body axes. This estimator is motivated by the need of having two vectors expressed in two different triads in order to determine attitude changes. Flight nonlinear simulations are performed to determine real attitude and compare it to the estimated attitude obtained from the algorithms. The applicability of the proposed solution for aircraft flight navigation, guidance, and control or for ballistic rocket terminal guidance, where attack and side-slip angles or total angle of yaw are relatively small, is also demonstrated.
This paper is organized as follows. Section 2 describes the problem in detail. Section 3 is devoted to the estimation of the gravity vector. The attitude determination algorithm is detailed in Section 4. Flight dynamics and sensor models are briefly described in Section 5. Section 6 shows simulation results. Conclusions and references follow.
2. Problem Description
Attitude determination is a fundamental field in aircraft engineering, as maneuvers in order to change vehicle trajectories are governed by aerodynamic forces, which depend directly on ship orientation angles. Furthermore, attitude in the artillery rocket terminal phase is vital as it determines factors as penetrability of payload or avoidability of countermeasures for the artillery rocket. Therefore, developing algorithms to improve attitude determination is a cornerstone in precise guidance, navigation, and control research.
One of the techniques to determine attitude is to calculate it from the director cosine matrix (DCM) which determines completely the rotation of a body between two reference frames. In order to obtain this DCM, a geometrical plane defined in both reference triads must be obtained. Every geometrical plane is defined by two vectors, and therefore, knowing the same two vectors expressed in the two reference frames, DCM may be calculated.
2.1. Triad Definition
Two axis systems are defined in order to express forces and moments: north-east-down (NED) axes and body (B) axes. NED axes are defined by subindex . is pointing north, is perpendicular to and pointing nadir, and is forming a clockwise trihedron. Body axes are defined by subindex . is pointing forward and contained in the plane of symmetry of the rocket, is perpendicular to pointing down and contained in the plane of symmetry of the rocket, and is forming a clockwise trihedron. The origin of body axes is located at the center of mass of the rocket.
2.2. Involved Vector Determination
If a GNSS sensor device is equipped on the aircraft, velocity vector can be directly expressed from the sensor velocity output on the NED triad. Expression for velocity in NED is defined in (1), where , , and are the components of this velocity in NED axes.
The same velocity vector expressed in body triad can be obtained from a set of accelerometers equipped on the ship, one on each of the axes. These devices are able to measure acceleration. After integration of each of the components along time, velocity vector is obtained, which is expressed in (2), where , , and are the components of the acceleration in body axes and describes angular velocity expressed in body axes.
Another vector expressed in two bases is required in order to define the rotation. Gravity vector is very simple to determine in NED triad as it is always parallel to , as it is shown in (3), where is the gravity acceleration, which is fixed constant in this model (9.81 m/s2). Note that precision may be increased using more complicated models, that is, it can be modeled depending on latitude and longitude of the aircraft.
The keystone of the presented attitude determination method is determining gravity vector in body axes. Methods to calculate this gravity vector are developed in [22]. For example, by determining the constant component of the measured acceleration employing a low pass filter, where jerk in body axes is calculated by derivation of acceleration, then, it is integrated in order to obtain the nonconstant component of acceleration, and, by subtracting this nonconstant component from the measured acceleration, gravity vector is estimated. However, this method is not valid when the aircraft rotates. Another method to obtain gravity vector is by integrating the mechanization equations [23] and manipulating consequently the resulting equations. Again, gyros are needed in order to implement this method.
In the next section, a new estimation method, which is valid for nonrotating and rotating aircrafts and which does not need gyros, is presented.
3. Gravity Vector Estimation
The gravity vector estimation method presented here depends on the aerial platform on which it will be employed. This means that the method must be adjusted for the aircraft of interest. Without loss of generality, the estimation method detailed in this section (and also Sections 5 and 6) is particularized for a four-canard-controlled spinning rocket. But, using the applicable flight mechanics equations, it may be applied to other aircraft.
The estimated gravity vector may be expressed as in (4), where its components in body axes are shown.
In order to describe the estimator, the gravity vector must be modeled as a subtraction between accelerations measured by accelerometers in body axes () and an estimation of the aerodynamic and control forces expressed also in body axes. Equation (5) describes this estimator as a particularization of Newton-Euler equations, where describes aerodynamic and inertial forces, describes control forces, and describes aircraft mass.
In order to compute , , , and must be estimated exclusively from data obtained from GNSS sensors and accelerometers. First of all, the Mach number () and angle of attack (), which are shown in (6), are estimated from accelerometer data. This hypothesis is completely valid for small angles of attack, which is the case for the flights described in this paper. The Mach number and angle of attack must be calculated in order to estimate forces and moments involved on aircraft flight. Speed of sound is obtained as a function of altitude () employing the International Standard Atmosphere (ISA) model. This altitude is directly obtained from the GNSS positioning output of the sensor.
Passive aerodynamic and inertial forces () are calculated according to where is the drag force, is the lift force, is the Magnus force, is the pitch damping force, and is the Coriolis force.
Each of these forces is estimated by the expressions in (8), where is the air density as a function of altitude using the ISA model, is the rocket caliber, is the drag force linear coefficient, is the drag force square coefficient, is the lift force linear coefficient, is the lift force cubic coefficient, is the Magnus force coefficient, and are the rocket inertia moments in body axes, is the pitch damping force coefficient, and is the earth angular velocity vector. Note that aerodynamic coefficients are dependent of the Mach number, which was calculated in (6).
In order to estimate the required forces to estimate the gravity vector, angular velocity is required. Newton-Euler equations are used in order to estimate angular velocity (see (9)), where is the rocket inertia tensor defined in (10), are the passive aerodynamic moments which actuate on the projectile, and are the moments generated by control surfaces. Passive aerodynamic moments are defined in (11), where is the overturn moment, is the pitch damping moment, is the Magnus moment, and is the spin damping moment.
Each of the aerodynamic moments is estimated in (12), where is the overturning moment linear coefficient, is the overturning moment cubic coefficient, is the pitch damping moment coefficient, is the Magnus moment coefficient, and is the spin damping moment coefficient.
Similarly, control forces and moments in body axes for each of the four fins must be estimated. For this purpose, the similar model as defined in [24] is employed, but with some modifications in the aerodynamics. The first modification is that the effective incidence aerodynamic speed on each of the four fins must be estimated.
The mathematical expressions for incidence aerodynamic speed decomposition are defined in where denotes the fin and and are the aerodynamic speed minus the fin-span component of the aerodynamic speed and the perpendicular component of to , respectively, , , , and .
The effective angle of attack for each fin () may be expressed for each fin as the sum of the local angle of attack, , and fin deflection, (see the following equation). where , , , and .
And finally, the control force provided by each fin () is given in where is the rocket caliber, is the air density, is the aerodynamic coefficient of the normal force for a fin, and is the characteristic surface of airfoil.
The total control force () is expressed as follows:
The control moment provided by each fin () may be expressed as follows: where is the longitudinal distance, parallel to , of the airfoil center of pressure () to the rocket center of mass (), which depends on the Mach number (). is the lateral distance, which is orthogonal to and parallel to for each fin, from the airfoil center of pressure to the rocket center of mass; it is supposed to be constant in this model.
The total control moment is given as follows:
Finally, note that each of the aerodynamic coefficients defined in this section is dependent on the Mach number; an interpolation (Table 1) based on experimental data for the geometry of the rocket analyzed is provided next.
After this previous mathematical development, the gravity vector can be estimated on body axes () just by accelerometer data and a set of coefficients that are preloaded on the rocket guidance, navigation, and control system.
4. Attitude Determination Algorithm
In order to simplify the calculus, an orthonormal base must be defined in both axis systems, B and NED. This orthonormal base is defined by unitary vectors , , and expressed in both bases and defined by the following expressions:
After defining an orthonormal base in the two systems of axes, the expression to determine the DCM is indicated in (20), where is a 3 × 3 square matrix composed of orthonormal vectors in the body triad, expresses the same concept in the NED triad, and is the director cosine matrix that transforms the NED triad into a body triad.
The DCM matrix can be solved from (20) as it is shown in (21). Note that employing an orthonormal base simplifies the calculation of the inverse matrix as it is the transposed matrix.
After obtaining the director cosine matrix, the rotation must be characterized. The most suitable method to express this rotation is through quaternions, as they avoid any possible singularities on the poles of rotation. The matrix equation that relates and quaternions is shown in (22), where for are the quaternions and for and are the DCM matrix coefficients.
The expressions in (23) may be used in order to solve quaternions. where , , , , and are defined in (24).
It is known that quaternions themselves are enough to express rotations without singularities, but it is also known that conceptually, they are difficult to be visualized. An easier manner to define these rotations is by means of Euler angles. Concretely, the most common aeronautical rotation is defined by the roll (), pitch (), and yaw () angles. The expressions that relate quaternions with Euler angles are defined by (25).
These Euler angles obtained from measurements may be used to characterize rotations and angular velocities in navigation guidance and control algorithms.
5. Flight Dynamics and Sensor Model
The flight dynamics model in [24] is employed. The equations of motion are integrated using a fixed time step Runge-Kutta scheme of fourth order to obtain a single flight trajectory. A set of ballistic simulated shots, with no wind perturbations, is performed varying shot angles from 15° to 75° taking 1° steps and also varying initial azimuth angles from 0° to 360° taking 45° steps. An example of these trajectories for a maximum range of 45°, minimum range of 15°, and maximum height of 75° for the whole set of azimuth shots is shown in Figure 1.
GNSS sensors are modeled, based on a u-blox series 7 sensor, adding a white noise of null average and maximum absolute value of mp = 3 m for position values and for velocity of mv = 0.01 m/s to the simulated position and velocity obtained from the model.
Accelerometers are modeled as second-order systems, based on MTi 1-series accelerometer, with a natural frequency of 300 Hz and a damping factor of 0.7. A white noise of null average and maximum absolute value of ma = 10−3 m/s2 and a bias of 10−3 m/s2 is added. Note that accelerometers will provide acceleration in body axes, which will be employed on the determination of the velocity vector in body axes.
The modeling of the sensors has been based on relatively low-cost equipment. This fact will induce significant errors on measurements, which will determine the suitability of the attitude calculation method. However, as it is shown in the following section, the results are accurate enough for terminal guidance tasks in artillery.
6. Simulation Results
Simulation results are presented using the nonlinear flight dynamics model [24] in order to demonstrate suitability of attitude determination algorithm. A set of ballistic simulated shots was performed as defined in the previous section, in order to compare real (denoted with subindex “real”) and estimated attitudes. MATLAB/Simulink R2016a on a desktop computer with a processor of 2.8 GHz and 8 GB RAM was used.
In order to compare the estimated attitude against the real or estimated attitude, gravity vector, and angular velocity vector, the following expressions are represented in Figures 2, 3, and 4: , , and for the attitude; for the gravity vector; and for the angular velocity. Note that real magnitudes are obtained by simulation results and are compared against the magnitudes obtained by algorithms defined on the previous sections which are represented by the second terms on the last expressions.
All these figures represent results for a shot angle of 45° and three different azimuth angles: 0°, 45°, and 90°. Note that a larger amount of simulations have been performed, as it was defined previously, but they are not shown in figures because the results are virtually the same and in order to simplify the comprehension of the paper.
Note that in the worst of the situations the difference between real and estimated attitude, which is the goal of these algorithms, is very precise for the pitch and yaw angles and precise for the roll angle which assures the applicability of the algorithms for the sensors employed and its associated noises for rocket terminal guidance applications.
In order to quantify the error made in each of the shots along the whole simulated trajectories, the root mean squared errors (RMSEs) for each of the 3 Euler angles (, , and ), for each of the three components of gravity vector in body axes (, , and ), and for each of the three components of angular velocity vector in body axes (, , and ) are calculated in (26), (27), and (28), respectively. This error depends on the shot launch angle () and azimuth ().
Figures 5, 6, and 7 show the previously defined root mean squared errors computed for the whole trajectory (z-axis) depending on shot angle (x-axis) and for the 8 azimuths (y-axis).
Note that the obtained root mean squared errors for all the measures of interest presented in this paper are small enough to consider the proposed algorithm for estimating attitude of high dynamic vehicles at a low cost. First of all is contained between 0° and , is contained between 0° and , and is contained between 0° and . These results verify that the algorithms applied to this concrete set of sensors and attitude determination vectors ( and for ) are more precise in the measurements related with the pitch and yaw angles and less precise for the roll angle. Note that in an axis-symmetric aircraft such as rockets or missiles, the roll angle can be handled by autopilot easily [24]. Root mean squared errors for gravity vector and angular velocity prove that estimations resolved by the methods explained in this paper are precise enough, as RMSEs are small.
Table 2 shows the root mean squared errors for the 3 Euler angles for different shot angles and azimuths, which is basically the numerical values of Figure 5. The first column shows the shot angle and the second one the displayed error (, , and ). The rest of the columns show the numerical values of the errors for different azimuths. In order to present a significant result for the whole method, root mean squared errors have been calculated for the values shown on the surface represented by Figure 5. In order to accomplish this task, the last column, named total, shows the root mean squared errors for the set of azimuths represented on the same row. Similarly, the last three rows represent the root mean squared errors for the set of shot angles represented on the same column. Finally, the three values of the bottom-right corner can be considered as an estimation of the whole error of the attitude calculation algorithms.
7. Conclusions
A novel low-cost approach for estimating attitude of aerial platforms has been proposed. This approach is based on the measurement of two different vectors: speed and gravity vectors. They are measured and estimated in two different reference frames using a GNSS sensor and a set of accelerometers. A new estimation method for obtaining the gravity vector in body axes has also been presented.
The presented algorithm is of special interest for high dynamic vehicles where other low-cost approaches, such as the ones using magnetometers, especially in ballistic rockets or unmanned air vehicles are not suitable due electronic interferences with other electronic equipment.
Nonlinear simulations based on ballistic rocket launches are performed to obtain ideal attitude and compare it to the attitude obtained by the proposed approach. Computational results show that errors are small enough to consider the presented algorithm as a good candidate for use within a control algorithm. Concretely, the total root mean square errors defined in Table 2 are below for the three attitude angles in the totality of tested cases. It is demonstrated that through the attitude determination algorithm proposed, precision is maintained. Note that the needed low-cost components make this approach highly interesting for nonreusable vehicles such as ballistic rockets.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
The authors would like to thank Lieutenant Colonel Jesús Sánchez (NMT) and “Instituto Tecnológico La Marañosa” of the National Institute for Aerospace Technology (INTA) for the solid modeling of the concept.