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

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.


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.

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 − XYZ, 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]. 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 − , the ballistic coordinate system − 2 2 2 , the projectile axis coordinate system − and the second projectile axis coordinate system − 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 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 − and − 2 2 are non-rolling coordinate systems that do not roll with the projectile. The axis of each coordinate system is the vertical axis of the projectile and the only difference between the coordinate planes and 2 2 is a turning angle [50].

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]: where ⃗⃗ is the geomagnetic field vector, 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 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 − , the ballistic coordinate system − 2 2 2 , the projectile axis coordinate system − and the second projectile axis coordinate system − 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 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 − and − 2 2 are non-rolling coordinate systems that do not roll with the projectile. The axis of each coordinate system is the vertical axis of the projectile and the only difference between the coordinate planes and 2 2 is a turning angle [50].

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]: where ⃗⃗ is the geomagnetic field vector, 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

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]: 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 1A ,ϕ S 1B ) and (ϕ S 2A ,ϕ S 2B ). The ratio Φ can be calculated as 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.
Sensors 2019, 19, x FOR PEER REVIEW 4 of 21 output signal of the magnetometer is zero, and the roll phase of the projectile represents the zerocrossing. There are two zero-crossings in a single cycle. The zero-crossing method uses two uniaxial magnetometers ( 1 and 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 ( 1 , 1 ) and ( 2 , 2 ). The ratio can be calculated as 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 twoelement 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.

Method for Estimating Pitch and Yaw
The movement of the projectile in air consists of two parts: the centroid motion and the aroundcentroid 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].

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: 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.

Method for Estimating Pitch and Yaw
The movement of the projectile in air consists of two parts: the centroid motion and the aroundcentroid 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].

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: 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 zz are derivatives of the static and equatorial damping moment coefficients.
The external moments include static and equatorial damping moments, and their vector forms are 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, 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: The component form of the static moment in the projectile axis coordinate system O − ξηζ is given as (7) 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: 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]: 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 If the influence of wind is ignored, the static moment components of the Oη and Oζ axes in Equation (7) can be written as Substituting Equation (9) into the above equation, we obtain Defining A m = ρSl 2A m z , the static moment components M zη and M zζ can be rewritten as In the same way, the equatorial damping moment in Equation (5) can be written in component form in the coordinate system O − ξηζ as follows: Similarly, under the condition of no wind, defining C m = − ρSld 2A m zz , the components M zzη and M zzζ of the equator damping moment can be rewritten as Considering both Equations (13) and (15), the components M η and M ζ of the total external moment can be rewritten as Substituting Equation (16) into (3), the dynamics constraint equations including flight-path angle, trajectory deflection angle and two Euler angles are obtained as follows: 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 ω ξ .

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 − XYZ are denoted as 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  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 ⃗⃗ and ⃗⃗ 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 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 is known before the experiment.
The projectile axis unit vector can be obtained by projecting the projectile axis vector onto the reference coordinate system − as shown in Figure 1.
The magnetic azimuth 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:

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 Thus, the geomagnetic unit vector can be expressed as 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 − XYZ as shown in Figure 1.
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:

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.

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 The measurement equation is The workflow of the RK4-UKF algorithm is as follows: • Calculation of the sigma point set • Prediction phase • Correction phase

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 The measurement equation is The workflow of the RK4-UKF algorithm is as follows: • Calculation of the sigma point set • Prediction phase Sensors 2019, 19, 5096

State Equation
Given the continuous nonlinear state equations in Equation (17), the state variables are written as Then, Equation (17) can be written as 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.

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): where measurement noise V is the Gaussian white noise, given as V ∼ N(0, R), and R = σ 2 σ M .

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. 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.
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.  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. 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.

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.

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 • .

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:

1.
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.

2.
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.

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. 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 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. Figure 10. The launch angle and direction of the experiment.
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 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. 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 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. Figure 10. The launch angle and direction of the experiment.
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 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 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 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. 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.

Figure 15.
Rotational speed calculated using the zero-crossing method. 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.

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 0 can be calculated either using the zero-crossing method, or based on the initial velocity as follows [50,51]: Figure 15. Rotational speed calculated using the zero-crossing method.

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]: 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. where the twist pitch of is rifling, is the diameter of the projectile and 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.

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.

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 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. 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.
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 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.

Discussion on Experimental Results
Based on the analysis of experimental data and filtering results, the following issues and relevant conclusions were obtained: 1. 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. 2. 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

Discussion on Experimental Results
Based on the analysis of experimental data and filtering results, the following issues and relevant conclusions were obtained: 1.
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.

2.
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. 3.
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.

4.
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.

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:

1.
Although the geomagnetic azimuth used in the proposed method was calculated using the zero-crossing method, any other method can also be used.

2.
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. 3.
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.