A Novel Attitude Estimation Algorithm Based on the Non-Orthogonal Magnetic Sensors.

Because the existing extremum ratio method for projectile attitude measurement is vulnerable to random disturbance, a novel integral ratio method is proposed to calculate the projectile attitude. First, the non-orthogonal measurement theory of the magnetic sensors is analyzed. It is found that the projectile rotating velocity is constant in one spinning circle and the attitude error is actually the pitch error. Next, by investigating the model of the extremum ratio method, an integral ratio mathematical model is established to improve the anti-disturbance performance. Finally, by combining the preprocessed magnetic sensor data based on the least-square method and the rotating extremum features in one cycle, the analytical expression of the proposed integral ratio algorithm is derived with respect to the pitch angle. The simulation results show that the proposed integral ratio method gives more accurate attitude calculations than does the extremum ratio method, and that the attitude error variance can decrease by more than 90%. Compared to the extremum ratio method (which collects only a single data point in one rotation cycle), the proposed integral ratio method can utilize all of the data collected in the high spin environment, which is a clearly superior calculation approach, and can be applied to the actual projectile environment disturbance.

1. Introduction

The geomagnetic field is the earth’s natural resource, which can provide a natural coordinate system for aeronautics, astronautics and marine applications. The vehicle’s attitude information can be achieved by using magnetic detection technology to measure the components of the geomagnetic field. A projectile body rotates about its longitudinal axis at very high speed. Due to the limitation of the measuring range of the triaxial gyroscope, it is difficult to use the gyroscope in an actual high-speed, high-spin projectile environment. In contrast, a magnetic sensor has advantages such as a rapid response speed, small size, high resistance to overload, and being cumulatively error free [1,2,3]. A magnetic sensor is applicable to the high-speed, high-spin general projectile attitude measurement [4,5]. Projectile attitude information includes heading angle, pitch angle and rolling angle, and its characteristic is an approximately constant heading angle. Therefore, the projectile attitude information can be estimated using roll angle and pitch angle.

There are several methods for measuring the projectile attitude by unitizing magnetic sensors, including the two-axis non-orthogonal magnetic sensor, three-axis orthogonal magnetic sensor [6,7], and four-axis magnetic sensor. The non-orthogonal magnetic sensors [8] utilize two mounted non-orthogonal magnetic sensors to collect geomagnetic field data during one projectile spin cycle. Because the heading angle and pitch angle are nearly invariable when a projectile rotates one cycle, the roll angle of the projectile can be calculated directly. There are two non-orthogonal magnetic sensor methods for calculating the pitch angle of spinning projectiles: the zero crossing method [9,10] and the extremum ratio method [11]. Both methods utilize the relationship between the eigenvalue ratios of the sensor output and the pitch angle to calculate the projectile attitude. Both of these methods can obtain the estimated value of projectile attitude at a special point. The difference between these methods is as follows: the zero crossing method utilizes phase information from two sensors to calculate the projectile attitude [12], whereas the extremum ratio method utilizes extremum value information from two sensors to calculate the projectile attitude. However, because those two methods use the special point of the sampled magnetic data to calculate the attitude value as the projectile rotates one cycle, the two methods are vulnerable to random disturbance. As the error of sampling value becomes smaller, the error of the estimated projectile attitude becomes smaller, and vice versa [13]. The twin-channel tracking differentiator [14] method utilizes a differential filter to process data from the triaxial magnetic sensors, thereby reducing the random disturbance encountered by the projectile attitude.

To reduce the random disturbance that exists in the extremum ratio method, an integral ratio method is proposed in this paper. This method can reduce the random disturbance by establishing the integral model and can thus achieve greater attitude calculation accuracy than the extremum ratio method.

The remainder of this paper is organized as follows. In Section 2, the projectile attitude measurement principle of magnetic sensor is introduced. In Section 3, the mathematical model of integral ratio method is established. The analytical expression of integral ratio method is derived in Section 4. In Section 5, the simulation results demonstrate that the accuracy of the proposed algorithm is significantly higher than that of the extremum ratio method. Conclusions are drawn in Section 6.

2. Attitude Measurement Principle of the Magnetic Sensors

The geographical coordinate is , the x axis points to the magnetic north, the y axis points to the sky, and the Z axis points to the magnetic east. As shown in Figure 1, vector H represents the magnitude and direction of the magnetic field in the geographical coordinate, Ψ is the heading angle, θ is the pitch angle, and γ is the roll angle.

Figure 1

Coordinate transformation schematic diagram.

The intersection angle between vector H and its projection in horizontal plane Px is geomagnetic inclination I. The carrier coordinate Ox can be achieved by rotating geographical coordinate Ox I degrees around the z axis. The x axis in Ox coincides with the magnetic field vector H, and H represents the geomagnetic magnitude and direction of Ox.

The instantaneous centroid of projectile is chosen as the origin of the projectile coordinate Ox. The Ox axis coincides with the lengthwise axis of projectile, pointing to the head, which is positive. Oy is perpendicular to Ox, and upward is positive. Oz is perpendicular to Ox, and the direction is determined by using the right hand rule. Vector H represents the magnitude and direction of the magnetic field in the projectile coordinate Ox.

The coordinates H, H, H of vector H in the projectile coordinate Ox can be represented as the transformation from carrier coordinates to projectile coordinates . where is the pitch angle, which includes magnetic dip, ; is the relationship that can be described by the direction cosine matrix.

Vector H in carrier coordinates Ox is where h is the magnitude of the vector H.

The coordinates H, H, H of vector H in the projectile coordinates Ox can be represented as [11]

Two non-orthogonal single-axis magnetic sensors, S1 and S2, are assembled separately on A and B along with axis Ox of the projectile. As shown in Figure 2, the axis Ox coincides with the longitudinal axis of the projectile. Two sensitive axes are in the plane Ox. The sensitive axis of S1 is assembled along with axis Oz, and the sensitive axis of S2 is assembled in the angle of with axis Ox.

Figure 2

Schematic diagram of mounted magnetic sensors.

According to Equation (4), the measured data of magnetic sensor S1 form H, and the measured data of magnetic sensor S2 are the combination of H and H As a result, the measured data H1 and H2, of two magnetic sensors, can be expressed by , and h:

The projectile rotates at nearly constant speed during one spinning circle, and the heading angle Ψ and pitch angle θ are nearly invariable. The roll angles during one spinning cycle can be calculated by obtaining one roll angle at a particular time. Assuming that the particular time is the zero point time of H1 and H2, the unknown magnetic field intensity scalar h in Equations (5) and (6) can be eliminated when the measured data of H1 or H2 are zero. The effects of environmental disturbance can thus be reduced.

If the measured data H1 = 0, then Equation (5) can be rewritten as

Furthermore, we can obtain where j = 1,2. The two components of Equation (8) cannot be zero simultaneously.

If the measured data H2 = 0, then Equation (6) can be rewritten as

If then

When , because , then . Thus, Equation (11) can be described as

Substituting Equation (12) into gives

When , , Equation (13) has real solution where , , and .

When , H2 has two zero points.

When , H2 has one zero point.

When , H2 does not have a zero point.

From the above, we can see that the value of roll angle γ can be calculated from Equations (8) and (14) when the outputs of two magnetic sensors, S1 and S2, are equal to zero. There are four zero-crossing points that can be used individually to calculate the corresponding roll angle γ. Consequently, the roll angles at any time during one cycle can be calculated. Assuming that the heading angle Ψ does not change [15], when the included angle λ between magnetic sensor S2 and axis Ox is given, the value of γ is related to the pitch angle . As a result, to further ensure the value of the roll angle γ, it is important to utilize the magnetic data obtained as the projectile rotates one cycle to calculate . Thus, assuming that the rotation speed remains constant, the calculation error of rolling angle is determined by the calculation error of the pitch angle as the projectile rotates one cycle.

3. Mathematical Model

3.1. Mathematical Model of the Extremum Ratio Method

Assume that the heading angle and the pitch angle are invariable in one cycle. The extremum ratio method utilizes the characteristics that there exist one maximum value and one minimum value in every cycle. The pitch angle after the projectile rotates one cycle can be calculated. During one projectile rotation cycle, the roll angle changes between 0 and . When the output curves H1 and H2 are at their extremum values, Equations (5) and (6) should meet the following conditions [6]:

The mathematic model of the integral ratio method can be described as where H1 and H2 are the maximum or minimum values achieved by H1 and H2, respectively. represents the ratio of H1 and H2. If the heading angle Ψ and the included angle λ between S2 and axis are known, then the estimated pitch angle can be calculated by Equation (17).

3.2. Mathematical Model of the Integral Ratio Method

The mathematical model of the extremum ratio method utilizes the ratio between two particular points on magnetic sensor curves to calculate the pitch angle . Its computational accuracy depends on the numerical precision of the sampling values H1 and H2 at the extreme point. Under practical conditions, due to various types of disturbance, some sampling values may have high error rates. In particular, if impulsive interference exists at any non-extreme point, a rather large sampling value will be generated. As a result, that value will be treated as the largest sampling value to calculate, thereby generating large errors.

A novel integral model is derived in this paper. During the rotation process, all of the samples of H1 and H2 are used to perform the integral calculation; subsequently, the ratio calculation is performed. The integral model is described as:

The estimated value of can be expressed as where N is the total sampling number during one projectile rotation cycle and n is the sampling point. , , and are true value of magnetic sensors S1 and S2, respectively. and are sampled value of S1 and S2, respectively. and are both Gaussian white noise with a mean of zero and a variance of . In Equation (19), one value can be calculated in every cycle using the quadratic sum and mean calculation with the sampling values from two magnetic sensors.

The mathematical expectation of can be described as where ; according to the properties of normal distribution, we can get

Therefore, , and

Substituting Equation (22) into Equation (20) gives

The mathematical expectation of can be given by

Similarly, the mathematical expectation of can be given by

Therefore, and are not unbiased estimators. Compared the mathematical expectation with unbiased estimation of and , the magnitude of single sampling error variance is amplified times. So the effect on parameter estimation is small and acceptable. and can be regarded as powers of two magnetic sensors’ sampled signal, and the power of the magnetic sensor signal is actually sampled as a characteristic value in the integral model to calculate the projectile attitude.

Thus, Equations (5) and (6) can be rewritten as where , , , , and .

The integral expression of can be described as

Applying integral calculation on each items in Equation (29) separately gives

Substituting Equations (30)–(32) into Equation (29) gives

The integral expression of can be described as

Applying the integral calculation on each items in Equation (34) separately gives

Substituting Equations (30)–(32), (35) and (36) into Equation (34) gives

Using Equations (33) and (37), the integral model can be described as

Substituting into Equation (38) gives

We define and . Then, Equation (39) can be rewritten as

The pitch angle is therefore given by

In Equation (41), the two components cannot be zero simultaneously.

The period of pitch angle is . The period of function with respect to pitch angle is . As the projectile rotates one cycle, one value and the corresponding four values can be obtained. Assuming that , Figure 3 gives the comparison between and the two integral outputs of and . It can be seen that when the pitch angle of projectile takes any value, feature values ( and) can be got from integral models and the estimated value of can also be obtained throughout missile flight. However, one may correspond to two or four possible values of , and the value of the pitch angle can be chosen based on the quadrant in which the pitch angle is located during the actual flight.

Figure 3

Comparison between the period of function () and the two integral outputs of feature values and .

4. Novel Integral Ratio Method

A novel integral ratio method is proposed in this paper, and the flow chart of this method is shown in Figure 4. Because this method utilizes all of the data collected during the projectile rotation process, it is convenient to preprocess those data.

Figure 4

The flow chart of integral ratio method.

4.1. Data Preprocessing

Noise error exists in the measured data from the magnetic sensor. To reduce the noise disturbance, the filtering algorithm is usually adopted to preprocess the data smoothly. A maximum least squares filtering algorithm is used.

When the samples number is less than the designated value p, the increasing memory filter is successively used to perform the filtering. p is chosen as 10 in this paper.

When the sample number is greater than p, the fixed memory filter can be used. p is the filter order of moving average part. where is the sampling value at time and is the least squares filtering value at time .

4.2. The Principle of the Integral Ratio Method

Extreme ratio method only using a single extremum of sample data as characteristic value for processing, which is susceptible to random interference. But the improved integration ratio method utilizes the sum of all data squares of projectile spin as a characteristic value for processing, which can reduce the random noise on the estimated values.

The implementation procedure of the proposed integral ratio method is given as follows:

Step 1: Preprocessing the data from the magnetic sensor using Equations (42) and (43).

The sampling values from the magnetic sensor are filtered by utilizing the maximum least squares filtering algorithm. After filtering the noise error, the value can be input into the integral model to calculate the model function .

Step 2: Calculating the integral value of and .

Determine the integral expression .

Determine the integral expression .

Step 3: Calculating the integral ratio ..

One value can be calculated according to Equation (19) by using the quadratic sum and the integral of the samples of two magnetic sensors.

Step 4: Calculating the pitch angle .

During one projectile rotation cycle, assuming that both the included angle and the heading angle are invariable, the pitch angle can be determined by Equation (41).

Step 5: Calculating the roll angle using Equations (8) and (14).

5. Simulation and Results

5.1. Comparison of the Algorithms

In this section, the noise restraining abilities of the extremum ratio method and the proposed integral ratio method are compared to each other. Assuming that the projectile’s heading angle is , the included angle between magnetic sensor S2 and axis is , and the data range of pitch angle is . The noise of the magnetic sensor is zero mean white noise. Figure 5 and Figure 6 show a comparison of the pitch angle error for different noise intensities ( and , respectively) after applying the extremum ratio and the proposed integral ratio methods.

Figure 5

Comparison of the pitch angle error with .

Figure 6

Comparison of the pitch angle error with .

Figure 5 and Figure 6 show that both algorithms have rather large errors around . The reason for the error is that the numerical difference of the magnetic component is too small near the zero point and the error of extremum is too large. The calculation error of the pitch angle also depends on the noise intensity. The error during the variance 0.01 is one order of magnitude greater than the error observed with a variance of 0.001. In the range of the pitch angle, the calculation error of the integral ratio method is much smaller than that of the extremum ratio method.

Table 1 shows the mean and variance of the pitch angle errors, which were calculated by applying the extremum ratio method and the integral ratio method, as well as the error comparison between these methods under different noise intensities. When the noise variance is 0.001, the mean and variance of the errors by applying integral ratio method can be reduced by 84.6% and 89.2%, respectively, compared with that of the errors obtained by applying the extremum ratio method. When the variance of noise is 0.01, the mean and variance of the errors by applying integral ratio method can be reduced by 90.2% and 96.1%, respectively, compared with that of the errors obtained by applying the extremum ratio method.

Table 1

The comparison of the pitch angle error.

Method	Error Mean		Error Variance
	(unit: 0.001 rad)		(unit: 0.001)
	σ2 = 0.001	σ2 = 0.01	σ2 = 0.001	σ2 = 0.01
Extremum ratio method	0.45389	11.2	0.017722	0.87621
Integral ratio method	0.070004	1.1	0.0019182	0.034168
Error reduction	84.6%	90.2%	89.2%	96.1%

Table 1 indicates that the mean and variance of the errors obtained by applying the integral ratio method are significantly smaller than those of the errors obtained by applying the extremum ratio method. In addition, as the random errors from the magnetic sensor become larger, the integral ratio method becomes clearly superior in terms of calculation accuracy. Thus, the integral ratio method is clearly superior when the magnetic sensor operates under an environment with a high level of disturbance.

The extremum ratio method utilizes the extremum to calculate the pitch angle as the projectile rotates one cycle. If the extremum encounters any disturbance, the calculation will produce errors. However, the integral ratio method uses the result of the square of the magnetic sensor output to perform the integral calculation. Because the integral calculation performs the function of a filter when the data encounters any white noise, the integral ratio method has strong noise suppression ability.

From the above-mentioned information, we can see that the proposed integral ratio method is clearly superior in both the mean value and the variance of the calculation error compared to the extremum ratio method. When the variance of random noise is 0.01, the error of the integral ratio method is only 10% of the error of the extremum ratio method.

5.2. Ballistic Simulation

The dynamic model of the projectile with magnetic sensor assembled is described as where x, y and z are the position components and , and are the velocity components.

The ballistic coefficient , D is the diameter of projectile, i is the elastic coefficient, and is the mass of projectile.

, is the standard virtual temperature. is the virtual temperature, is the function of air density, kg/m3, , and is the flight height.

is the resistance function, and .

is the gravitational acceleration, m/s2, and m.

The equations of the projectile attitude is denoted as

The initial conditions of the flight trajectory are set as follows: D = 152 mm, L = 1300 mm, m = 52.8 kg, rad/s, and J = 7.5 kg·m2. The muzzle velocity is 550m/s, the initial heading angle is , the pitch angle is , and the rolling angle is . The included angle between the pitch angle of magnetic sensor S2 and axis is . The noise of the magnetic sensor is white noise, with a mean of 0 and a variance of 0.01 or 0.001. Both the magnetic sensor sampling period and the updating period of the projectile attitude are 1 ms.

Figure 7 shows the curves of the heading angle, pitch angle and rolling angle. The simulation curve of the rolling angle runs only 1 s, from 50 s to 51 s.

Figure 7

The curves of the attitude angle during projectile flight.

Figure 8 shows the curve of the projectile attitude obtained by utilizing the proposed integral ratio method. The simulation curve of the rolling angle runs only 1 s. The result shows that the proposed integral ratio method can produce an estimated value of the projectile attitude during the whole trajectory. Figure 9 shows the error curve of the projectile attitude by utilizing the integral ratio method.

Figure 8

The attitude angle calculated by the integral ratio method.

Figure 9

Attitude angle calculating error of the integral ratio method.

Figure 10 indicates that the calculation error of the pitch angle obtained by utilizing the integral ratio method is much smaller than that obtained by utilizing the extremum ratio method.

Figure 10

The pitch angle calculated by the extremum ratio method and the integral ratio method.

Table 2 lists the variance of the calculation error of the pitch angle and the rolling angle obtained by utilizing two different methods. The data show that the variance of calculation error of pitch angle obtained by utilizing the integral ratio method is one order of magnitude smaller than that obtained by utilizing the extremum ratio method. The variance of the calculation error of the rolling angle obtained by utilizing the integral ratio method is only half that obtained by utilizing the extremum ratio method.

Table 2

The error variance comparison between the extremum ratio method and the proposed integral ratio method.

Method	Pitch Angle	Roll Angle
	(unit: 0.001)	(unit: 0.001)
Extremum ratio method	1.4526	0.20321
Integral ratio method	0.18652	0.093807
Error reduction	87.2%	53.8%

6. Conclusions

To solve the computational problem of a projectile’s pitch angle, the extremum ratio method can be used to select the extreme point to perform the calculation. According to the problem relating to the fact that the extremum ratio method will result in a large error under a noise disturbance, a novel integral ratio method was proposed, and the corresponding mathematical model was derived. In this paper, the computational expression of pitch angle was derived by utilizing the integral ratio method. The simulation results provided comparison diagrams between the extreme ratio method and the integral ratio method under different noise intensities. Compared to the extremum ratio method, the integral ratio method has several advantages. When the variance of noise becomes larger, the integral ratio method is clearly superior in terms of the calculation error. In addition, the extremum ratio method only selects one data point in every spin cycle, whereas the integral ratio method selects multi-group data in every spin cycle; as a result, multiple data processing algorithms can be performed to get more accurate attitude estimation. Furthermore, the integral model allows for digital processing. An existing digital device can be utilized to perform high-speed, high-precision AD (analog digital conversion) sampling.