Next Article in Journal
Progressive Temporal-Spatial-Semantic Analysis of Driving Anomaly Detection and Recounting
Previous Article in Journal
A Sensor for Spirometric Feedback in Ventilation Maneuvers during Cardiopulmonary Resuscitation Training
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Novel Method for Estimating Pitch and Yaw of Rotating Projectiles Based on Dynamic Constraints

1
School of Energy and Power Engineering, Nanjing University of Science & Technology, Nanjing 210094, China
2
Beijing Key Laboratory of High Dynamic Navigation Technology, Beijing Information Science & Technological University, Beijing 100101, China
*
Author to whom correspondence should be addressed.
Sensors 2019, 19(23), 5096; https://doi.org/10.3390/s19235096
Submission received: 15 September 2019 / Revised: 13 November 2019 / Accepted: 19 November 2019 / Published: 21 November 2019
(This article belongs to the Section Physical Sensors)

Abstract

:
This paper addresses the difficult problem of measuring the attitude of a high-spinning projectile and presents a novel method for estimating the pitch and yaw angles of the projectile in flight. The method is based on analysis of the external moment of the rotating projectile during flight and theoretical derivations obtained from the dynamics’ equations. First, the principle of zero-crossing method is introduced, which explains the process of geomagnetic azimuth and roll measurements by the non-orthogonal geomagnetic sensor combination. Then, the dynamics constraint equations between the Euler angles and flight-path angle, trajectory deflection angle of the projectile are derived using the dynamics equations of the projectile rotating around the centroid, and analysis of the flight characteristics of the projectile in stable flight. Next, the spatial orientation relationship between pitch, yaw angles and magnetic azimuth is established based on the physical principle of geomagnetic azimuth. Finally, the pitch and yaw angles are estimated using the unscented Kalman filter (UKF), with the dynamics constraint equations serving as the driving equations. In the UKF prediction stage, the Runge-Kutta method is used to discretize the state equation that improves the prediction accuracy. Simulation results show that the proposed method can be used to accurately calculate the pitch and yaw angles, and results of experimental data processing also verify the feasibility of the proposed method for real-world applications.

1. Introduction

Due to ever-increasing accuracy requirements for precision-guided weapons, acquisition of accurate flight attitude information of projectiles has become crucially important for analyzing their flight dynamics, as well as providing support for the navigation & guidance system. At present, the most commonly used attitude measurement methods rely on solar sensors [1,2], angular rate gyros [3,4,5], inertial measurement units (IMU) [6,7,8] and magnetometers [9,10,11]. Among these methods, the solar sensors work effectively only under good weather conditions, the angular velocity gyros have an upper limit on the rotational speed of the projectile, and the IMU suffer from error accumulation. Therefore, for measuring attitude of high-speed rotating projectiles, special working conditions, i.e., high temperature, high pressure, high overload, and high speed, as well as the requirements of low cost and small size preclude the use of many sensors. The magnetometer can be widely used in attitude estimation of rotating objects [12,13,14,15,16,17,18] after undergoing a calibration and compensation process [19,20,21], thanks to its features of reliable performance, low cost, and no error accumulation.
Many domestic and international scholars have carried out research on measuring the attitude of rotating projectiles using magnetometers. Wilson M. described several attitude measurement solutions based on combinations of multiple low-cost sensors [22] and conducted an in-depth research on the application of magnetometers in smart bomb [23]. Changey S. et al. conducted lab and flight experiments to verify the effectiveness of an algorithm that estimates roll of the projectile based on data acquired by two magnetometers [24,25]. Maley J. proposed a full attitude estimation method for spin-stabilized projectiles based on steady-state Kalman filtering [26]. Rogers J. et al. developed a low-cost orientation estimator for smart bombs equipped with magnetometers and thermopiles. This orientation estimator, with the advantage of not relying on GPS and other state feedbacks, estimates Euler angles and rotation rates using extended Kalman filtering (EKF) [27]. Most of these conventional methods calculate the Euler angles using the transfer matrix between the relevant coordinate systems [28,29]. Since the solution of the three Euler angles is non-independent, data fusion with devices such as built-in sensors and thermopiles is required in the calculation process [30,31,32]. However, the methods dependent on other built-in devices require sufficient internal space and specific working conditions. Apart from increasing production cost, adding extra electronic devices in a fixed space also increases the errors and noise of the output signal. Researchers in the United States first proposed a method for measuring the magnetic azimuth and rotational speed of a rotating projectile using two uniaxial magnetometers with a specific angle to the projectile axis [33,34]. This low cost and high precision method is called the ‘zero-crossing method’. Based on this method, experts at the Nanjing University of Science and Technology proposed an extreme value ratio method that combined non-orthogonal magnetic sensors, and conducted research in related areas [35,36]. Researchers at the Beijing Information Science and Technology University proposed a novel phase shift ratio method based on the extreme value ratio method [37]. In this study, a novel technique was developed based on the zero-crossing method. The moment applied to the rotating projectile flying in the air was analyzed and then the angular relationship contained in the external moment was extracted based on the ballistic characteristics of the projectile in stable flight. Subsequently, the constraint relationship between pitch, yaw angles and flight-path angle, trajectory deflection angle were deduced.
The EKF has been used in several studies for attitude estimation [38,39,40]. As it requires linearization of a nonlinear system before performing Kalman filtering, it is suitable for linear or weakly nonlinear systems. The EKF also has stringent requirements on the accuracy of the filter parameters, and involves calculation of the Jacobian matrix that is cumbersome. Compared with the EKF, the unscented Kalman filter (UKF) exhibits good robustness in the presence of nonlinearity and uncertainty [41], therefore, it is better at dealing with complex models with high nonlinearity [42,43,44,45,46] and has been used widely in recent years. It is necessary to discretize a continuous system when the filtering algorithm is applied in computers, and the discretization method and discrete step-size directly affect the filtering accuracy. When the step-size is large, the discrete models processed by the conventional methods such as the Euler method significantly differ from the continuous models. On the other hand, reducing the discrete step-size increases the computational complexity. When the fourth-order classical Runge-Kutta method [47,48,49] is used as the discretization method, the reliance on discrete step-size is reduced greatly. Consequently, the discrete models become closer to the theoretical continuous models, and the filtering accuracy is improved.
The method proposed in this paper works as follows: First, the dynamic constraint equations between pitch, yaw angles and flight-path angle, trajectory deflection angle are derived and used as the state model. Then, the geomagnetic vector and the projectile axis vector are simultaneously projected onto the reference coordinate system to obtain the spatial orientation relationship between the pitch, yaw angles and magnetic azimuth, and a measurement model based on geomagnetic azimuth is constructed. Finally, the pitch and yaw angles of the rotating projectile are estimated using the UKF algorithm, which utilizes the fourth-order classical Runge-Kutta method as the discretization method. The effectiveness of the proposed method is verified through simulations and processing of experimental data.

2. Definition of Coordinate System and Principle of Zero-Crossing Method

2.1. Coordinate Systems

To establish the differential equations of projectile dynamics, we use the approach described in [50] to introduce several basic coordinate systems: the reference coordinate system O X Y Z , the ballistic coordinate system O X 2 Y 2 Z 2 , the projectile axis coordinate system O ξ η ζ and the second projectile axis coordinate system O ξ η 2 ζ 2 . Figure 1 shows the angular relationships between these coordinate systems. In the figure, the angles θ and ψ are the Euler angles of the pitch and yaw, respectively, the angle θ a is the angle between the velocity vector and the horizontal plane, the angle ψ 2 is the angle between the velocity vector and the vertical plane, respectively, i.e., flight-path angle and trajectory deflection angle, and δ is the total attack angle of the projectile. Figure 2 further illustrates the pitch component δ 1 and the yaw component δ 2 of the total attack angle.
Both O ξ η ζ and O ξ η 2 ζ 2 are non-rolling coordinate systems that do not roll with the projectile. The axis O ξ of each coordinate system is the vertical axis of the projectile and the only difference between the coordinate planes O η ζ and O η 2 ζ 2 is a turning angle β [50].

2.2. Principle of Zero-Crossing Method

When a projectile and a uniaxial magnetometer with an angle of λ to the projectile axis rotate together in the earth’s magnetic field, the instantaneous field strength along the sensitive axis of the magnetometer is as follows [34]:
M S = c o s ( λ ) | M | c o s ( σ M ) + s i n ( λ ) | M | s i n ( σ M ) s i n ( ϕ )
where M is the geomagnetic field vector, σ M is the magnetic azimuth, i.e., the angle between the projectile axis and the direction of geomagnetic field, ϕ is the roll of the projectile and λ is the angle between the sensitive axis of the uniaxial magnetometer and the projectile axis. The rotating projectile is considered to be flying steadily in the air, and the output signal of the magnetometer changes periodically. When the sensitive axis is orthogonal to the direction of the earth’s magnetic field, the output signal of the magnetometer is zero, and the roll phase of the projectile represents the zero-crossing. There are two zero-crossings in a single cycle.
The zero-crossing method uses two uniaxial magnetometers ( S 1 and S 2 ) with different mounting angles, as shown in Figure 3. With the mounting angles of 90° and 60°, the two magnetometers are in a coplanar relation with the projectile axis, i.e., they have equal initial roll phases. Four zero crossings can be extracted using the output signals of the two magnetometers, which results in two pairs of rolls given as ( φ S 1 A , φ S 1 B ) and ( φ S 2 A , φ S 2 B ). The ratio Φ can be calculated as
Φ = | φ S 2 B φ S 2 A φ S 1 B φ S 1 A |
The magnetic azimuth of the projectile relative to the geomagnetic field during the flight can be determined based on the magnetic azimuth-ratio diagram plotted beforehand, and the roll angular rate and the roll phase angle of the projectile can be obtained using the recorded zero-crossing time.
The three-element attitude information, i.e., roll, pitch, and yaw angles, is converted into two-element attitude information, i.e., roll angle and geomagnetic azimuth. Subsequently, the roll angle information is separated to provide the possibility for secondary processing of the pitch and yaw angles, which is a strength of the zero-crossing method.

3. Method for Estimating Pitch and Yaw

The movement of the projectile in air consists of two parts: the centroid motion and the around- centroid motion. The former is mainly characterized by the position and velocity of the projectile, and is governed by the law of centroid movement. The latter is characterized by the attitude of the projectile, and is governed by the theorem of angular momentum [50].

3.1. Dynamics Constraint Equations

It is necessary to analyze the moment of external forces relative to the center of mass during the flight of the projectile. When there is no wind and the projectile shape does not cause any aerodynamic eccentricity, only the static and equatorial damping moments need to be considered. References [50,51] provide the dynamics equations of the projectile undergoing around-centroid motion. A new set of dynamic equations based on the specific problem are obtained as follows:
{ ω η ˙ = 1 A M η C A ω ξ ω ζ ω ζ ˙ = 1 A M ζ + C A ω ξ ω η θ ˙ = ω ζ c o s ( ψ ) ψ ˙ = ω η m z ˙ = 0 m z z ˙ = 0
where M η and M ζ are components of the external moments in the projectile axis coordinate system O ξ η ζ , A and C are coefficients of the moment of inertia, ω ξ , ω η and ω ζ are projection components of the angular velocity on coordinate system O ξ η ζ and m z and m z z are derivatives of the static and equatorial damping moment coefficients.
The external moments include static and equatorial damping moments, and their vector forms are
M z = ρ S l 2 v r m z 1 s i n δ r ( v r × ξ )
M z z = ρ v r S l d m z z ω / 2
where M z and M zz are the static moment vector and equatorial damping moment vector, ρ is the air density, S is the cross-sectional area of the projectile, l is the projectile length, d is the projectile diameter, v r is the velocity vector of the projectile relative to the wind, m z is the static moment coefficient, δ r is the relative attack angle, ξ is the unit vector of the axis O ξ of the coordinate system O ξ η ζ and ω is the projectile oscillation angular velocity. When there is no wind, v r is equal to v , and δ r is the attack angle δ .
For a small attack angle, m z = m z δ r , and the form of the static moment vector in Equation (4) can be rewritten as follows:
M z = ρ S l 2 v r m z ( v r × ξ )
The component form of the static moment in the projectile axis coordinate system O ξ η ζ is given as
M z ξ = 0 M z η = ρ S l 2 v r m z v r ζ M z ζ = ρ S l 2 v r m z v r η
where v r η and v r ζ are components of the relative velocity v r in the coordinate system O ξ η ζ . Let the components of the relative velocity v r in the coordinate system O ξ η 2 ζ 2 be denoted as v r η 2 and v r ζ 2 , the relationship between two components is as follows:
v r η = v r η 2 c o s β + v r ζ 2 s i n β v r ζ = v r η 2 s i n β + v r ζ 2 c o s β
For a normally flying projectile, as the attack angle and the ballistic deflection are small, δ 1 , δ 2 , ψ , ψ 2 and θ θ a have small values. Thus, the following relationship holds [50]:
β 0 ;   δ 1 θ θ a ;   δ 2 ψ ψ 2
As shown in Figure 2, the rotation relationship between the ballistic coordinate system O X 2 Y 2 Z 2 and the second projectile axis coordinate system O ξ η 2 ζ 2 leads to v r η 2 = v δ 1 and v r ζ 2 = v δ 2 . Consequently, Equation (8) can be further written as
v r η = v r η 2 = v δ 1 v r ζ = v r ζ 2 = v δ 2
If the influence of wind is ignored, the static moment components of the O η and O ζ axes in Equation (7) can be written as
M z η = ρ S l 2 v m z ( v δ 2 ) M z ζ = ρ S l 2 v m z ( v δ 1 )
Substituting Equation (9) into the above equation, we obtain
M z η = ρ S l 2 v 2 m z ( ψ ψ 2 ) M z ζ = ρ S l 2 v 2 m z ( θ θ a )
Defining A m = ρ S l 2 A m z , the static moment components M z η and M z ζ can be rewritten as
M z η = A A m v 2 ( ψ ψ 2 ) M z ζ = A A m v 2 ( θ θ a )
In the same way, the equatorial damping moment in Equation (5) can be written in component form in the coordinate system O ξ η ζ as follows:
M z z ξ = ρ v r 2 S l d m z z ω ξ 0 M z z η = ρ v r 2 S l d m z z ω η M z z ζ = ρ v r 2 S l d m z z ω ζ
Similarly, under the condition of no wind, defining C m = ρ S l d 2 A m z z , the components M z z η and M z z ζ of the equator damping moment can be rewritten as
M z z η = A C m v ω η M z z ζ = A C m v ω ζ
Considering both Equations (13) and (15), the components M η and M ζ of the total external moment can be rewritten as
M η = M z η + M z z η = A A m υ 2 ( ψ ψ 2 ) + A C m υ ω η M ζ = M z ζ + M z z ζ = A A m υ 2 ( ϑ θ a ) + A C m υ ω ζ
Substituting Equation (16) into (3), the dynamics constraint equations including flight-path angle, trajectory deflection angle and two Euler angles are obtained as follows:
{ ω η ˙ = A m v 2 ( ψ ψ 2 ) + C m v ω η C A ω ξ ω ζ ω ζ ˙ = A m v 2 ( ϑ θ a ) + C m v ω ζ + C A ω ξ ω η ϑ ˙ = ω ζ cos ( ψ ) ψ ˙ = ω η A m ˙ = 0 C m ˙ = 0
where the projectile’s flight speed v , flight-path angle θ a and trajectory deflection angle ψ 2 can be calculated using the ballistic radar data. The axial angular velocity of the projectile is denoted by ω ξ .

3.2. Relationship between Pitch, Yaw Angles and Magnetic Azimuth

When the projectile is flying in the air, its instantaneous attitude relative to the earth’s magnetic field can be represented by the pitch, yaw, magnetic dip and magnetic declination. Since the projectile’s position information can be detected by the radar, the geomagnetic field information, including magnetic dip and magnetic declination of the projectile’s position can be calculated based on the geomagnetic field model. Therefore, the magnetic azimuth σ M only contains two pieces of unknown information, i.e., the pitch and yaw.
The shooting direction is denoted as α N , and the influence of the meridional convergence angle is ignored. The unit vectors on the three axes of the reference coordinate system O X Y Z are denoted as i , j and k . The geomagnetic field vector is projected onto the reference coordinate system.
As shown in Figure 4, the geomagnetic vector is described by the north-east-down (NED) coordinate system. Taking the northern hemisphere as an example, the geomagnetic unit vector and its horizontal projection are M and M N , respectively, the magnetic declination is D , north to the east is positive, the magnetic dip is I and the downward direction is positive.
Thus, the geomagnetic unit vector can be expressed as
M = c o s ( I ) c o s ( a N D ) i s i n ( I ) j c o s ( I ) s i n ( a N D ) k
where both the magnetic declination D and the magnetic dip I can be calculated based on the spherical harmonics model of the geomagnetic field, and the shooting direction α N is known before the experiment.
The projectile axis unit vector ξ can be obtained by projecting the projectile axis vector onto the reference coordinate system O X Y Z as shown in Figure 1.
ξ = c o s ( θ ) c o s ( ψ ) i + s i n ( θ ) j + c o s ( θ ) s i n ( ψ ) k
The magnetic azimuth σ M is the angle between the geomagnetic unit vector and the first projectile axis unit vector. It can be calculated as follows using the vector included angle cosine formula:
c o s ( σ M ) = M ξ | M | | ξ | = c o s ( I ) c o s ( α N D ) c o s ( θ ) c o s ( ψ ) s i n ( I ) s i n ( θ ) c o s ( I ) s i n ( α N D ) c o s ( θ ) s i n ( ψ )

3.3. Estimation of Dip and Yaw

When the derived dynamics constraint equations are used as the driving equations, the pitch and yaw angles of the rotating projectile can be estimated based on the spatial relationship between the magnetic azimuth and the Euler angles. A block diagram of the attitude estimation method is shown in Figure 5. The built-in magnetometer provides accurate geomagnetic signals. The magnetic azimuth and rotational speed serve as the inputs to the estimation algorithm and can be calculated using the zero-crossing method. The rotational speed can be used to obtain the roll. The radar collects the velocity and position information, and calculates the geomagnetic field information of the entire trajectory based on the geomagnetic field model to provide support to filtering. The initial firing elements are used to simulate the magnetic azimuth of the initial section of the trajectory and perform initial filtering calibration using the calculated magnetic azimuth as a reference. Finally, the pitch and yaw angles are estimated using the improved UKF algorithm.

4. Design of UKF

The UKF mainly consists of two phases: the prediction phase and the correction phase. In the prediction phase, a set of state prediction points based on the sigma points should be generated. Since the state equation is a continuous model, discretization needs to be carried out that can directly affect the accuracy of the filtering results. The Runge-Kutta method is often used in ballistic calculations as it outperforms other methods in terms of discretization accuracy under the same step-size. In this paper, the fourth-order classical Runge-Kutta method is used for state estimation in the prediction stage. Therefore, the filtering algorithm used in this paper is called the RK4-UKF algorithm.
Assume that the state equation of a continuous nonlinear system is
X ˙ k = f [ X k 1 , k 1 ] + W k 1
The measurement equation is
Y k = h [ X k , k ] + V k
The workflow of the RK4-UKF algorithm is as follows:
  • Calculation of the sigma point set
    X k 1 0 = x ^ k 1 ;   X k 1 i = { x ^ k 1 + ( n + λ ) P x i = 1 , 2 , , n x ^ k 1 ( n + λ ) P x i = n + 1 , , 2 n
  • Prediction phase
    k 1 = f ( X k 1 i ) ;   k 2 = f ( X k 1 i + h 2 k 1 ) ;   k 3 = f ( X k 1 i + h 2 k 2 ) ;   k 4 = f ( X k 1 i + h k 3 ) X k / k 1 i = X k 1 i + h 6 ( k 1 + 2 k 2 + 2 k 3 + k 4 ) x ^ k / k 1 = i = 0 2 n W i m X k / k 1 i ;   P k / k 1 = i = 0 2 n W i c [ X k / k 1 i x ^ k / k 1 ] [ X k / k 1 i x ^ k / k 1 ] T + Q k
  • Correction phase
    Y   k / k 1 i = h ( X k / k 1 i ) ;   y ^ k / k 1 = i = 0 2 n W i m Y   k / k 1 i
    P ( Y Y ) k / k 1 = i = 0 2 n W i c [ Y   k / k 1 i y ^ k / k 1 ] [ Y   k / k 1 i y ^ k / k 1 ] T + R k
    P ( X Y ) k / k 1 = i = 0 2 n W i c [ X k / k 1 i x ^ k / k 1 ] [ Y k / k 1 i y ^ k / k 1 ] T
    K k = P ( X Y ) k / k 1 P ( Y Y ) k / k 1 1
    X ^ k = X ^ k / k 1 + K k ( Y k y ^ k / k 1 )
    P k = P k / k 1 K k P ( Y Y ) k / k 1 K k T

4.1. State Equation

Given the continuous nonlinear state equations in Equation (17), the state variables are written as
x = [ x 1    x 2    x 3    x 4    x 5    x 6 ] = [ ω η    ω ζ    θ    ψ    A m    C m ]
Then, Equation (17) can be written as
X ˙ = f ( x ) = ( x 5 v 2 ( x 4 ψ 2 ) + x 1 x 6 C A ω ξ x 2 x 5 v 2 ( x 3 θ a ) + x 2 x 6 + C A ω ξ x 1 x 2 c o s x 4 x 1 0 0 ) + W
As the nonlinear equations given in Equation (31) only approximately describe the around-centroid motion of the projectile, there will be certain errors. Therefore, Gaussian white noise W N ( 0 , Q ) is introduced to model these errors.

4.2. Measurement Equation

The magnetic azimuth is represented by a measured variable y = ( σ M ) . The measurement equation can be constructed as follows, based on Equation (20):
y = h ( x ) + V = arccos ( c o s ( I ) c o s ( α N D ) c o s ( x 3 ) c o s ( x 4 ) s i n ( I ) s i n ( x 3 ) c o s ( I ) s i n ( α N D ) c o s ( x 3 ) s i n ( x 4 ) ) + V
where measurement noise V is the Gaussian white noise, given as V N ( 0 , R ) , and R = ( σ σ M 2 ) .

5. Simulation and Experimental Results

5.1. Simulation and Analysis

5.1.1. Simulation

The calculation steps for the magnetic azimuth and roll are described in detail in [35]. Therefore, these steps will not be repeated here and instead, only the simulation results will be given. The focus of simulations in this study is the estimation of pitch and yaw angles. Assume that the projectile is launched from a location of (E100°, N39°) with an initial velocity of 800 m/s, a shooting angle of 60° and a shooting direction of 100°. The pitch and yaw components of the angular velocity equal to 2 rad/s are added to simulate the initial disturbance at the time of launch, and the ballistic data are simulated using the 6D ballistic equations. Then, the geomagnetic signal output information of the trajectory is simulated through conversion between the relevant coordinate systems. Finally, the magnetic azimuth is calculated using the zero-crossing method and serves as the true value.
A Gaussian white noise d ~ N(0,0.5°) is added to the true value of the magnetic azimuth to serve as measurement value. Figure 6a shows the simulated and true values of the initial 1 s of the ballistics, and Figure 6b shows the discrepancy between the true and simulated measured values. It can be observed that the maximum error of the simulated measured value is about ±1.6°, which is much larger than the measurement error described in [34].
The pitch angle and yaw angle are estimated using the RK4-UKF algorithm and compared with the corresponding true values, as shown in Figure 7. Figure 7a shows the estimation of the pitch and yaw angles of the entire ballistic. It can be seen from the figure that the projectile flied for more than 100 s. The estimated value is close to the true value in the entire ballistic, and the estimated effect is satisfactory. From the law of yaw angle movement, it can be seen that the existence of dynamic equilibrium angle causes the yaw angle to deviate to the right from the trajectory deflection angle (the yaw angle was defined as positive to the right, so the yaw angle is positive) in the midcourse because of the right-hand twist of projectile. Figure 7b shows the estimation error of the pitch and yaw angles, both of which are mostly in the range of (−0.2° ~ 0.2°). The errors are slightly larger in the beginning of the trajectory phase, and fluctuate in the midcourse. Figure 7c,d show the estimated and true values of the pitch and yaw angles within 1 s of the initial phase of ballistic. The dual-circular motion law of the projectile can be seen clearly from the figures. The pitch and yaw angles oscillate periodically around flight-path angle and trajectory deflection angle, respectively. This oscillation is a slow-circular motion, with a low frequency and continuously diminishing amplitude. At the same time, the projectile axis oscillates periodically around the dynamic balance axis in a fast-circular fashion, with a continuously diminishing nutation amplitude.
To show movement of the body axis, the oscillation trajectory of the projectile is constructed based on the pitch and yaw angles, as shown in Figure 8. The projectile starts from a position of (60°,0°) and undergoes a counterclockwise dual-circular motion. The estimated oscillation trajectory of the projectile obtained using the RK4-UKF algorithm is consistent with the true trajectory. For the sake of clarity, only five cycles of the fast-circular motion in the initial phase and a slow circular motion cycle consisting of their centers are shown in the figure. Each red dot in the figure shows the approximate center of the fast-circular motion and the dotted line passing through the center represents the slow-circular motion, i.e., the motion trajectory of the dynamic balance axis of the projectile. The slow-circular motion is also in the counterclockwise direction. The dual-circular motion of the projectile is prominent.

5.1.2. Monte Carlo Simulation

In order to avoid the contingency of the RK4-UKF simulation results, a Monte Carlo simulation was performed. The RK4-UKF was run 1000 times, and the initial value of the state variable and the noise of the measurement value were randomly changed within a reasonable range before each run. After 1000 runs, the error of pitch and yaw angles estimated by RK4-UKF are tested in terms of the mean, mean square error and maximum of the absolute value, respectively. At the same time, the normal curve fitting was performed on the Monte Carlo test results. The test results are shown in Figure 9.
For the estimation error of the pitch angle, the expectation of mean is about 9.587 × 10−5 degree, the expectation of mean square error is about 0.065°, and the expectation of maximum of the absolute value is about 0.2644°. For the estimation error of the yaw angle, the expectation of mean is about −1.675 × 10−4 degree, the expectation of mean square error is about 0.05984°, and the expectation of maximum of the absolute value is about 0.1985°. Figure 9c is the maximum of the absolute value of the estimation error, i.e., the maximum estimation error. According to the principle of normal distribution, the maximum estimation error of the pitch angle does not exceed 0.457°, and the maximum estimation error of the yaw angle does not exceed 0.297°.

5.1.3. Analysis of simulation results

The simulation results show that the filtering result of the RK4-UKF algorithm is consistent with the true value, with a small error on the order of 10−1 degrees. The Monte Carlo simulation also shows that the pitch angle error estimated by this method does not exceed 0.457°, and the yaw angle error does not exceed 0.297°. The following conclusions were obtained based on the simulation results:
  • When analyzing the moment applied to the projectile during flight, it is assumed that the projectile shape has no eccentricity and there is no wind. Among external moments, only the static and equatorial damping moments are considered, while the smaller Magnus moment is ignored. Moreover, since the attack angle is small, the resulting small ballistic deviation from the firing surface allows approximations of δ 1 θ θ a ;   δ 2 ψ ψ 2 during the state equation derivation process. Thus, there is uncertainty in the adjustment of the state noise parameters.
  • Under the same sampling step-size, the fourth-order classical Runge-Kutta discretization method results in smaller discretization errors compared to other method such as the Euler method, and its discrete equations are close to the continuous model.

5.2. Experiment and Analysis

5.2.1. Experiment

It is impossible to observe directly and record the attitude of a flying projectile using current technology when the range reaches tens of kilometers or hundreds of kilometers. However, based on the pattern of projectile motion and the Lyapunov stability principle [50,51], the stability of the flying projectile can be maintained only when the following two conditions are met: a) The directions of the slow circular motions of the lateral attitude parameters including pitch and yaw are consistent with the velocity direction during the flight of projectile; b) The projectile axis undergoes periodic nutation around the velocity direction and the amplitude diminishes continuously. Therefore, it is feasible to verify the effectiveness of the attitude estimation method using the flight-path angle and trajectory deflection angle measured by a radar or a GPS device.
The experimental verification was conducted at a shooting range. During the experiment, the weather was good and windless, with the presence of a few clouds. The field layout of the verification experiment is shown in Figure 10. The reference coordinate system is the North-Up-East coordinate. The elevation angle θ 0 was 15.3° and the direction of fire α N was 103.3555°. A velocity radar was set up near the artillery location to measure the projectile velocity and an air balloon was launched to collect the meteorological data.
The measurement system used in the experiment consisted of a CPU, a magnetometer unit consisting of two single-axis magnetometers, a data acquisition unit, a signal processing unit, a communication unit, a power supply and other auxiliary unit. Figure 11 shows the block diagram of the measurement system. The geomagnetic unit collected the original voltage signal, the signal processing unit carried out signal conversion and processing, the CPU was responsible for signal processing, and the communication unit was used for transmitting and receiving instructions. Figure 12 shows the photos of the measurement system. The system was mounted inside the standard projectile to form the assembly, as shown in Figure 13a.
The assembly was mounted at the front section of a standard warhead. Its interior was fixed using solid glue and protected with a non-magnetic cover. This arrangement enabled it to withstand the strong impact and large overload during the launch stage, ensuring normal operation of the measurement components. The assembly was recycled after launching, as shown in Figure 13b.
The initial velocity of the projectile measured by the radar was 744.4 m/s. During the flight, the magnetometer recorded the variation of geomagnetic intensity in each axial direction, and the magnetic azimuth and rotational speed of the projectile were calculated using the zero-crossing method. The calculation results are shown in Figure 14.
Figure 14a shows the variation of magnetic azimuth along the entire trajectory. In the initial section of the trajectory, the projectile’s oscillation amplitude is relatively large due to the influence of initial disturbances. Then, the oscillation amplitude decreases gradually. Towards the end section of the trajectory, the projectile begins to oscillate again, and the oscillation amplitude increases continuously. There are two reasons behind this phenomenon: First, as the rotational speed of the projectile decreases, the gyroscopic effect of the projectile’s rotation diminishes gradually. Consequently, the dynamic stability of the projectile reduces gradually, eventually causing oscillations. Second, the change in the projectile’s velocity from supersonic to subsonic also causes oscillations. Radar data show that the projectile’s velocity was equal to 338.69 m/s (about Mach 1) around 24 s, which is in the transonic region.
Figure 14b shows the variation of magnetic azimuth during the initial 2 s section of the trajectory. The initial value of the magnetic azimuth is 117.6°, the minimum and maximum values are 112.4° and 127° in the first cycle of the slow circular motion, respectively, and the oscillation amplitude is about 7.3°. The pattern of dual-circular motion of the projectile can also be clearly seen from the figure. The oscillation amplitude diminishes continuously irrespective of fast or slow circular motion.
Figure 15 shows the variation of rotational speed of the projectile along the entire trajectory. The initial rotational speed of the projectile is 1486 rad/s, which finally drops to about 900 rad/s. This speed drop is fast at first and then gradually slows down.

5.2.2. Initial Alignment

Initial alignment needs to be performed at first to determine the initial value of the filter. To carry out this alignment, the theoretical trajectory is simulated using the 6D rigid body ballistic equations based on the initial firing elements. Then, the geomagnetic information of the entire trajectory is obtained through conversion between the relevant coordinate systems. With this information, the theoretical magnetic azimuth angle is calculated using the zero-crossing method and compared with the measured value. The initial firing conditions and moment coefficients should be adjusted continuously until the theoretical magnetic azimuth becomes roughly consistent with the measured value.
Among the initial launch conditions, the position and velocity of the projectile are provided by the radar, and the elevation angle and direction of fire are known in advance. The initial rotational speed w ξ 0 can be calculated either using the zero-crossing method, or based on the initial velocity as follows [50,51]:
ω ξ 0 = 2 π v 0 η d
where η the twist pitch of is rifling, d is the diameter of the projectile and v 0 is the initial velocity of the projectile after leaving the gun muzzle. Therefore, it is necessary to only adjust the initial angular velocity and moment coefficient of the pitch and yaw directions to obtain simulated magnetic azimuth curve that matches the measured magnetic azimuth curve, as shown in Figure 16.

5.2.3. Estimation of Pitch and Yaw Using RK4-UKF

The pitch and yaw angles of the projectile were estimated using the designed RK4-UKF algorithm, and the estimated values were compared with the flight-path angle and trajectory deflection angle obtained from the radar data, as shown in Figure 17a. The flight time of the projectile is about 32.4 s. The estimated value of the pitch angle is consistent with the flight-path angle and decreases with the decrease of the ballistic trajectory. The estimated value of the yaw angle is consistent with the trajectory deflection angle. Under the external action, the deflection direction is right. The amplitude of the projectile axis is large at the initial phase of ballistic because of initial disturbance, then decreases continuously, and increases again and at the terminal phase. Figure 17b,c show the motion of the projectile attitude clearly. The slow circular motion components of the projectile’s lateral attitude parameters, i.e., the pitch and yaw angles, oscillate around flight-path angle and trajectory deflection angle, respectively with continuously diminishing amplitudes, which fits with the pattern of flight stability of the projectile.
For a rotating projectile flying steadily in the air, the following approximations can be made: The discrepancy between the calculated pitch θ and flight-path angle θ a is taken as the pitch component δ 1 of the attack angle, and the discrepancy between the yaw ψ and trajectory deflection angle ψ 2 is taken as the yaw component δ 2 of the attack angle, as shown in Figure 18. These approximations are also helpful for measuring the attack angle of the projectile.

5.2.4. Discussion on Experimental Results

Based on the analysis of experimental data and filtering results, the following issues and relevant conclusions were obtained:
  • In a real-world scenario, the actual static and equatorial damping moment coefficients of the flying projectile often deviated from the theoretical values that were determined based on the projectile design. Therefore, it is necessary to adjust the theoretical moment coefficient when performing the initial alignment based on trajectory simulation.
  • Oscillation is bound to happen during the descending section of the actual trajectory. The larger the shooting angle, the larger the oscillation amplitude, which is the well-known Mayevsky problem. By contrast, the end section of the theoretical trajectory is free of oscillation. This is because the motion of projectile axis is obtained based on the pure kinematics theory, which assumes that the moment of momentum vector coincides with the projectile axis. Therefore, there is no projectile axis swing problem during the end section of the theoretical trajectory. The oscillation phenomenon that occurs during the descending section of the actual trajectory can be attributed to two factors: (1) A reduction of the gyro stability factor due to decreased rotational speed; (2) The dramatic change of the aerodynamic load that causes the projectile to oscillate when the projectile flight speed is in the transonic region.
  • The method for estimating pitch and yaw angels proposed in this paper is based on the constraints of dynamics equations of projectile. Through proper approximations, the relationship between the attitude and velocity angles can be determined, i.e., the slow-motion terms of the lateral attitude of the projectile are consistent with the velocity direction. This is also the basis for determining the rationality of the filtering results.
  • Limited by the current attitude measurement technology and experimental conditions, the true value of the projectile attitude cannot be obtained in the field experiment, and the accuracy of the estimation cannot be quantified. The experiment is mainly to verify the feasibility of the method in practical engineering applications. The method is proven to be feasible and effective through the analysis of the flight stability of the projectile. The quantification of estimation error by designing the verification experiment is the focus of the next step in the future.

6. Conclusions

In this paper, a novel method for estimating the pitch and yaw angles of a flying projectile was developed. Based on the analysis of the flight characteristics and external moment of the rotating projectile in steady flight, the dynamics constraint equations of the lateral attitude Euler angles and the velocity angles were derived using the projectile dynamics equation without relying on conversions between the relevant coordinate systems. The relationships between the pitch, yaw, and magnetic azimuth were established based on the spatial vector relationship. Finally, the pitch and yaw angles were estimated using the RK4-UKF algorithm. The feasibility and effectiveness of the proposed method were verified using simulation and experimental results, and different issues arising in the simulations and experiments were analyzed and discussed. The following points are worth noting:
  • Although the geomagnetic azimuth used in the proposed method was calculated using the zero-crossing method, any other method can also be used.
  • As the proposed method deals with high-spinning projectiles in steady flight, the magnetic declination and dip during the projectile flight can be obtained in two ways: (1) Calculation based on the geomagnetic model; (2) Calculation using the measured launch location based on the assumption that the magnetic declination and dip are constant at each location.
  • The object studied in this paper is the idealized projectile, which only considers the static moment and the equatorial damping moment, and assumes that there is no wind. The influence of the wind field model, Magnus moment and the moment caused by the shape asymmetry on the attitude of projectile will be considered in the future works, which makes the simulation model more accurate and improves the accuracy of the estimation.
The proposed method seeks to break through the dynamic characteristics of projectile and opens up new directions for developing attitude estimation methods. Catering to the need of developing intelligent bombs, it is expected to play an important role in the navigation and guidance of artillery shells and high-spinning rockets, and precision control of flying objects.

Author Contributions

L.A. proposed the idea, developed the method, and written the first version of this paper. L.W. guided the paper writing, reviewed the paper and offered the funding. N.L. revised the paper, discussed the idea, guided the paper writing and reviewed the paper before submission. J.F. offered the funding. Y.Z. discussed the idea of this paper.

Funding

This research was supported by the National Natural Science Foundation of China [61603191, 61903241], China Postdoctoral Science Foundation [2018M630424], Science Foundation of State Key [614260403010117], Jiangsu Province Post Doctoral Fund [1401023C].

Acknowledgments

The author gratefully acknowledges the helpful comments and suggestions of reviewers, which have improved the presentation. If you want the source code, please e-mail me, and I will share the all the materials except the experiment data.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Hepner, D.J.; Harkins, T.E. Determining inertial orientation of a spinning body with body-fixed sensors. Proc. SPIE Int. Soc. Opt. Eng. 2000, 4025, 68–78. [Google Scholar]
  2. Springmann, J.C.; Sloboda, A.J.; Klesh, A.T.; Bennett, M.W.; Cutler, J.W. The attitude determination system of the RAX satellite. Acta Astronaut. 2012, 75, 120–135. [Google Scholar] [CrossRef]
  3. Lefferts, E.J.; Markley, F.L.; Shuster, M.D. Kalman Filtering for Spacecraft Attitude Estimation. J. Guid. Control. Dyn. 1982, 5, 417–429. [Google Scholar]
  4. Pittelkau, M.E. Kalman Filtering for Spacecraft System Alignment Calibration. J. Guid. Control. Dyn. 2001, 24, 1187–1195. [Google Scholar] [CrossRef]
  5. Crassidis, J.L.; Markley, F.L. Three-Axis Attitude Estimation Using Rate-Integrating Gyroscopes. J. Guid. Control. Dyn. 2016, 39, 1513–1526. [Google Scholar] [CrossRef]
  6. Park, C.G.; Kim, K.; Kang, W.Y. UKF Based In-Flight Alignment Using Low Cost IMU. In Proceedings of the Aiaa Guidance, Navigation, & Control Conference & Exhibit, Keystone, CO, USA, 21–24 August 2006. [Google Scholar]
  7. Rhudy, M.; Gross, J.; Gu, Y.; Napolitano, M.R. Fusion of GPS and Redundant IMU Data for Attitude Estimation. In Proceedings of the Aiaa Guidance Navigation & Control Conference, Minneapolls, MN, USA, 13–16 August 2012. [Google Scholar]
  8. Ettl, J.; Schmidt, A.P.; Turner, J.; Kim, D.G. Using data fusion of DMARS-R-IMU and GPS data for improving attitude determination accuracy. In Proceedings of the International Conference on Space Operations, Daejeon, Korea, 16–20 May 2016. [Google Scholar]
  9. Changey, S.; Beauvois, D.; Fleck, V. A Mixed Extended-Unscented Filter for attitude estimation with magnetometer sensor. In Proceedings of the American Control Conference, Minneapolis, MN, USA, 14–16 June 2006; Institute of Electrical and Electronics Engineers: New York, NY, USA, 2006; pp. 2892–2897. [Google Scholar]
  10. Abdelrahman, M.; Park, S.Y. Integrated Attitude Determination and Control System via Magnetic Measurements and Actuation. Acta Astronaut. 2011, 69, 168–185. [Google Scholar] [CrossRef]
  11. Ivanov, D.S.; Ovchinnikov, M.Y.; Penkov, V.I.; Roldugin, D.S.; Doronin, D.M.; Ovchinnikov, A.V. Advanced Numerical Study of the Three-Axis Magnetic Attitude Control and Determination with Uncertainties. Acta Astronaut. 2017, 132, 103–110. [Google Scholar] [CrossRef]
  12. Psiaki, M.L.; Oshman, Y. Spacecraft Attitude Rate Estimation from Geomagnetic Field Measurements. J. Guid. Control. Dyn. 2003, 26, 244–252. [Google Scholar] [CrossRef]
  13. Christin, S.H. Satellite Attitude Determination Using Magnetometer Data Only. In Proceedings of the 47th AIAA Aerospace Sciences Meeting including The New Horizons Forum and Aerospace Exposition, Orlando, FL, USA, 5–8 January 2009. [Google Scholar]
  14. Roberts, T.M.; Lynch, K.A.; Clayton, R.E.; Disbrow, M.E. Magnetometer-Based Attitude Determination for Deployed Spin-Stabilized Spacecraft. J. Guid. Control. Dyn. 2017, 40, 2941–2947. [Google Scholar] [CrossRef]
  15. Halil Ersin Söken, S.S. Real-Time Attitude-Independent Magnetometer Bias Estimation for Spinning Spacecraft. J. Guid. Control Dyn. 2017, 41, 276–279. [Google Scholar] [CrossRef]
  16. Celani, F. Spacecraft Attitude Stabilization using Magnetorquers with Separation between Measurement and Actuation. J. Guid. Control Dyn. 2015, 39, 2184–2191. [Google Scholar] [CrossRef]
  17. Robinson, J.B.; Richie, D.J. Stabilization and Attitude Determination Methods for FalconSAT-3. J. Spacecr. Rocket. 2016, 507–519. [Google Scholar] [CrossRef]
  18. Rogers, J.; Costello, M.; Harkins, T.; Hamaoui, M. Effective Use of Magnetometer Feedback for Smart Projectile Applications. Navigation 2011, 58, 203–219. [Google Scholar] [CrossRef]
  19. Sedlak, J.E. Iterative magnetometer calibration. In Proceedings of the AIAA/AAS Astrodynamics Specialist Conference, Keystone, CO, USA, 21–24 August 2006. [Google Scholar]
  20. Hu, X.; Pang, H.; Fu, L.; Pan, M. Magnetometer calibration improvement using wavelet and genetic algorithm. IEEJ Trans. Electr. Electron. Eng. 2016, 11, S130–S137. [Google Scholar] [CrossRef]
  21. John, L.C.; Kok-Lam Lai Richard, R.H. Real-Time Attitude-Independent Three-Axis Magnetometer Calibration. J. Guid. Control Dyn. 2005, 28, 115–120. [Google Scholar]
  22. Wilson, M.J. Onboard Attitude Determination for Gun-Launched Projectiles. In Proceedings of the 43rd AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, USA, 10–13 January 2005. [Google Scholar]
  23. Wilson, M.J. Projectile Navigation and the Application to Magnetometers; Dissertations & Theses—Gradworks; University of Delaware: Newark, DE, USA, 2007. [Google Scholar]
  24. Changey, S.; Pecheur, E.; Wey, P. Real time Estimation of Supersonic Projectile Roll Angle using Magnetometers: In-lab Experimental Validation. IFAC Proc. Vol. 2009, 42, 123–127. [Google Scholar] [CrossRef]
  25. Changey, S.; Pecheur, E.; Bernard, L.; Sommer, E. Real time estimation of projectile roll angle using magnetometers: In-flight experimental validation. In Proceedings of the Position Location & Navigation Symposium, Myrtle Beach, SC, USA, 23–26 April 2012. [Google Scholar]
  26. Maley, J.M. Efficient Attitude Estimation for a Spin-Stabilized Projectile. J. Guid. Control. Dyn. 2015, 39, 339–350. [Google Scholar] [CrossRef]
  27. Rogers, J.; Costello, M. A Low-Cost Orientation Estimator for Smart Projectiles Using Magnetometers and Thermopiles. Navigation 2012, 59, 9–24. [Google Scholar] [CrossRef]
  28. Bar-Itzhack, I.Y.; Idan, M. Recursive attitude determination from vector observations Euler angle estimation. J. Guid. Control. Dyn. 1987, 10, 152–157. [Google Scholar] [CrossRef]
  29. Crassidis, J.L.; Markley, F.L.; Cheng, Y. Survey of Nonlinear Attitude Estimation Methods. J. Guid. Control Dyn. 2007, 30, 12–28. [Google Scholar] [CrossRef]
  30. Theil, S.; Appel, P.; Schleicher, A. Low Cost, Good Accuracy—Attitude Determination using Magnetometer and Simple Sun Sensor. In Proceedings of the 17th Annual AIAA/USU Conference on Small Satellites, Logan, UT, USA, 11–14 August 2003. [Google Scholar]
  31. Teshnizi, M.M.; Shirazi, A. Attitude estimation and sensor identification utilizing nonlinear filters based on a low-cost MEMS magnetometer and sun sensor. Aerosp. Electron. Syst. Mag. IEEE 2015, 30, 20–33. [Google Scholar] [CrossRef]
  32. Bekkeng, J.K.; Psiaki, M. Attitude Estimation for Sounding Rockets Using Microelectromechanical System Gyros. J. Guid. Control. Dyn. 2008, 31, 533–542. [Google Scholar] [CrossRef]
  33. Thompson, A.A. A Procedure for Calibrating Magnetic Sensors; Army Research Laboratory (ARL): Adelphi, MD, USA, 2002; Volume 524, pp. 1–2. [Google Scholar]
  34. Harkins, T.E.; Hepner, D.J. MAGSONDE (patent pending): A device for making angular measurements on spinning projectiles using magnetic sensors. Proc. SPIE Int. Soc. Opt. Eng. 2000, 4025, 60–67. [Google Scholar]
  35. Li, D.; BU, X. Attitude Measurement on High-spinning Projectile Using Magnetic Sensors and Acceletometers. Trans. Nanjing Univ. Aeronaut. Astronaut. 2008, 25, 106–112. [Google Scholar]
  36. Zhu, J.; Wu, P.; Bo, Y. A Novel Attitude Estimation Algorithm Based on the Non-Orthogonal Magnetic Sensors. Sensors 2016, 16, 730. [Google Scholar] [CrossRef]
  37. Zhao, H.; Su, Z.; Liu, F.; Li, C.; Li, Q. Magnetometer-Based Phase Shifting Ratio Method for High Spinning Projectile’s Attitude Measurement. IEEE Access 2019, 7, 22509–22522. [Google Scholar] [CrossRef]
  38. Boutayeb, M.; Sébastien, C.; Bara, J. A strong tracking extended Kalman observer for projectile attitude and position estimation. Proc SPIE 2007, 6538, 65381G-11. [Google Scholar]
  39. Allik, B.; Ilg, M.; Zurakowski, R. Ballistic Roll Estimation using EKF Frequency Tracking and Adaptive Noise Cancellation. IEEE Trans. Aerosp. Electron. Syst 2013, 49, 2546–2553. [Google Scholar] [CrossRef]
  40. Changey, S.; Fleck, V.; Beauvois, D. Projectile attitude and position determination using magnetometer sensor only. Proc. SPIE Int. Soc. Opt. Eng. 2005, 13, 802–806. [Google Scholar]
  41. Brink, K.M. Unscented Partial-Update Schmidt–Kalman Filter. J. Guid. Control. Dyn. 2018, 41, 929–935. [Google Scholar] [CrossRef]
  42. Julier, S.J.; Uhlmann, J.K.; Durrant-Whyte, H.F. A New Method for the Nonlinear Transformation of Means and Covariances in Filters and Estimators. IEEE Trans. Autom. Control 2000, 45, 477–482. [Google Scholar] [CrossRef]
  43. Julier, S.J.; Uhlmann, J.K. Unscented Filtering and Nonlinear Estimation. Proc. IEEE 2004, 93, 401–422. [Google Scholar] [CrossRef]
  44. Lisano, M.E. Nonlinear Consider Covariance Analysis Using a Sigma-Point Filter Formulation. In Proceedings of the Guidance and Control 2006: Proceedings of the 29th Annual AAS Rocky Mountain Guidance and Control Conference, Univelt, San Diego, CA, USA, 4–8 February 2006. [Google Scholar]
  45. Lisano, M.E. Comparing Consider-Covariance Analysis with Sigma-Point Consider Filter and Linear-Theory Consider Filter Formulas. In Proceedings of the 20th International Symposium on Space Flight Dynamics, Annapolis, MD, USA, 24–28 September 2007. [Google Scholar]
  46. Stauch, J.; Jah, M. Unscented Schmidt–Kalman Filter Algorithm. J. Guid. Control. Dyn. 2015, 38, 117–123. [Google Scholar] [CrossRef]
  47. Zingg, D.; Chisholm, T. Runge-Kutta methods for linear problems. In Proceedings of the 12th Computational Fluid Dynamics Conference, San Diego, CA, USA, 19–22 June 1995. [Google Scholar]
  48. Hempel, P.R. Application of implicit Runge-Kutta methods. In Proceedings of the AIAA & AAS, Astrodynamics Conference, Minneapolis, MN, USA, 13–16 August 2012. [Google Scholar]
  49. Aristoff, J.M.; Poore, A. Implicit Runge-Kutta Methods for Orbit Propagation. In Proceedings of the AIAA & AAS, Astrodynamics Specialist Conference, Minneapolis, MN, USA, 13–16 August 2012; pp. 13–16. [Google Scholar]
  50. Han, Z. Exterior Ballistics for Projectiles and Rockets; Beijing Institute of Technology Press: Beijing, China, 2008; pp. 127–157. [Google Scholar]
  51. Pu, F.; Rui, X. Exterior Ballistics; National Defense Industry Press: Beijing, China, 1980; pp. 131–176. [Google Scholar]
Figure 1. The angular relationships between these coordinate systems.
Figure 1. The angular relationships between these coordinate systems.
Sensors 19 05096 g001
Figure 2. Coordinate system O X 2 Y 2 Z 2 turns to Coordinate system O ξ η 2 ζ 2 .
Figure 2. Coordinate system O X 2 Y 2 Z 2 turns to Coordinate system O ξ η 2 ζ 2 .
Sensors 19 05096 g002
Figure 3. Installation diagram of two uniaxial magnetometers.
Figure 3. Installation diagram of two uniaxial magnetometers.
Sensors 19 05096 g003
Figure 4. Orientation map of the reference coordinate system and geomagnetic vector.
Figure 4. Orientation map of the reference coordinate system and geomagnetic vector.
Sensors 19 05096 g004
Figure 5. diagram of attitude estimation.
Figure 5. diagram of attitude estimation.
Sensors 19 05096 g005
Figure 6. Simulation of magnetic azimuth measurement. (a) Magnetic azimuth measurement, (b) Measuring noise.
Figure 6. Simulation of magnetic azimuth measurement. (a) Magnetic azimuth measurement, (b) Measuring noise.
Sensors 19 05096 g006
Figure 7. Estimated and true values of pitch and yaws. (a) Estimation of the pitch and yaw angles in the entire ballistic, (b) Estimation error of the pitch and yaw angles, (c) Estimation of the pitch angle within 1 s of the initial phase of ballistic, (d) Estimation of the yaw angle within 1 s of the initial phase of ballistic.
Figure 7. Estimated and true values of pitch and yaws. (a) Estimation of the pitch and yaw angles in the entire ballistic, (b) Estimation error of the pitch and yaw angles, (c) Estimation of the pitch angle within 1 s of the initial phase of ballistic, (d) Estimation of the yaw angle within 1 s of the initial phase of ballistic.
Sensors 19 05096 g007aSensors 19 05096 g007b
Figure 8. Oscillation trajectory of the projectile.
Figure 8. Oscillation trajectory of the projectile.
Sensors 19 05096 g008
Figure 9. Results of Monte Carlo simulation. (a) Mean of estimation errors,(b) Mean square error of estimation errors,(c) Maximum of the absolute value of estimation errors.
Figure 9. Results of Monte Carlo simulation. (a) Mean of estimation errors,(b) Mean square error of estimation errors,(c) Maximum of the absolute value of estimation errors.
Sensors 19 05096 g009
Figure 10. The launch angle and direction of the experiment.
Figure 10. The launch angle and direction of the experiment.
Sensors 19 05096 g010
Figure 11. The block diagram of the measurement system.
Figure 11. The block diagram of the measurement system.
Sensors 19 05096 g011
Figure 12. The photos of the measurement system. (a) Top view, (b) 45° view, (c) Side view
Figure 12. The photos of the measurement system. (a) Top view, (b) 45° view, (c) Side view
Sensors 19 05096 g012
Figure 13. The photos of assembly. (a) Assembly before experiment, (b) Recycled assembly.
Figure 13. The photos of assembly. (a) Assembly before experiment, (b) Recycled assembly.
Sensors 19 05096 g013
Figure 14. The magnetic azimuth calculated using the zero-crossing method. (a) Magnetic azimuth, (b) Magnetic azimuth in 1s.
Figure 14. The magnetic azimuth calculated using the zero-crossing method. (a) Magnetic azimuth, (b) Magnetic azimuth in 1s.
Sensors 19 05096 g014
Figure 15. Rotational speed calculated using the zero-crossing method.
Figure 15. Rotational speed calculated using the zero-crossing method.
Sensors 19 05096 g015
Figure 16. Initial alignment.
Figure 16. Initial alignment.
Sensors 19 05096 g016
Figure 17. Pitch and yaw angles estimated by RK4-UKF. (a) Estimation of pitch and yaw of the entire ballistic, (b) Estimation of pitch angle in 2 s, (c) Calculation of yaw angle in 2 s.
Figure 17. Pitch and yaw angles estimated by RK4-UKF. (a) Estimation of pitch and yaw of the entire ballistic, (b) Estimation of pitch angle in 2 s, (c) Calculation of yaw angle in 2 s.
Sensors 19 05096 g017aSensors 19 05096 g017b
Figure 18. Pitch and yaw components of the attack angle. (a) Pitch component of the attack angle, (b) Yaw component of the attack angle.
Figure 18. Pitch and yaw components of the attack angle. (a) Pitch component of the attack angle, (b) Yaw component of the attack angle.
Sensors 19 05096 g018

Share and Cite

MDPI and ACS Style

An, L.; Wang, L.; Liu, N.; Fu, J.; Zhong, Y. A Novel Method for Estimating Pitch and Yaw of Rotating Projectiles Based on Dynamic Constraints. Sensors 2019, 19, 5096. https://doi.org/10.3390/s19235096

AMA Style

An L, Wang L, Liu N, Fu J, Zhong Y. A Novel Method for Estimating Pitch and Yaw of Rotating Projectiles Based on Dynamic Constraints. Sensors. 2019; 19(23):5096. https://doi.org/10.3390/s19235096

Chicago/Turabian Style

An, Liangliang, Liangming Wang, Ning Liu, Jian Fu, and Yang Zhong. 2019. "A Novel Method for Estimating Pitch and Yaw of Rotating Projectiles Based on Dynamic Constraints" Sensors 19, no. 23: 5096. https://doi.org/10.3390/s19235096

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop