Abstract
The accurate measurement of roll angular rate for high spinning projectile has long been a challenging problem. Aiming to obtain the accurate roll angular rate of high spinning projectile, a novel extraction and filter algorithm, BSCZT-KF, is proposed in this paper. Firstly, a compound angular motion model of high spinning projectile is established. According to the model, we translate the roll angular rate measurement problem into a frequency estimation problem. Then the improved CZT algorithm, BSCZT, was employed to realize an accurate estimation of the narrowband signal frequency. Combined with the peak detection method, the BSCZT-KF algorithm is presented to further enhance the frequency estimation accuracy and the real-time performance. Finally, two sets of actual flight tests were conducted to verify the effectiveness and accuracy of the algorithm. The test results show that the average error of estimated roll angular rate is about 0.095% of the maximum of roll angular rate. Compared with the existing methods, the BSCZT-KF has the highest frequency estimation accuracy for narrowband signal.
1. Introduction
In recent years, with the development of science and technology, conventional artillery shells are developing in the direction of guidance and intelligence. The impact accuracy and overall combat effectiveness of the projectiles have been significantly improved, and conventional artillery shells have entered the era of low-cost, high-precision strikes. The precise guidance of the projectile cannot be separated from the accurate measurement of the real-time attitude of the projectile. However, the environment in which the projectile is located is relatively harsh and usually has the following characteristics. High overload. The launching overload of the projectile is generally above 10,000g. For example, the maximum axial launch overload of the US “copper snake” terminal guided projectile can reach up to 15000g [1]. High spinning. The projectile rotates at high speed around its longitudinal axis during flight, which can reduce the interference of thrust eccentricity, aerodynamic asymmetry, and uneven mass distribution on the stability of the projectile. Therefore, most of the current guided projectiles adopt the rolling mechanism. For instance, the rotation speed of a certain type of 155mm projectiles is as high as 250r/s. Such high rotation speed is a test for the range and resolution of the gyroscope. Limited space. The internal structure of the guided projectile is more complex than conventional shells, and the installation space of each component is limited.
Under the condition of high overload, high rotation speed, and limited space, the accurate measurement of the attitude of high spinning projectiles has not been effectively solved. According to the published literature, the currently used attitude measuring sensors mainly have gyroscope [2], accelerometer [3], magnetometer [4], GNSS [5], continuous wave radar [6], infrared sensor [7], solar azimuth sensor [8], and so on. Among them, the gyroscope is a typical angular rate sensor. However, the gyroscope cannot measure the excessive rotation speed of the projectiles with a limited range. Accelerometers can be sensitive to gravity vectors and are widely used in attitude estimation of unmanned aerial vehicles. However, because the motion of the projectiles during the flight is complicated, it is impossible to find its attitude through the gravity vector. The magnetometer is extremely suitable for measuring the attitude of high spinning projectiles with small size, large range, and strong overload resistance. But the magnetometer is easily interfered by the external environment and eddy current effects, and the calibration of magnetometer is cumbersome. The GNSS can provide information on the position and velocity of the projectile, so that the attitude information of the projectile can be calculated, but the GNSS data update rate is low and the signal may be lost. HD cameras, infrared sensors, and solar azimuth sensors are susceptible to environmental and weather conditions, so their application scenarios are limited.
In order to measure the attitude of high spinning projectiles, various scholars have proposed many solutions using the above sensors. Yu et al. [9] measured the attitude of high spinning projectile with a non-orthogonally mounted magnetometer. They found a one-to-one correspondence between the extreme value ratio of the sensor output and the pitch angle. Thus the pitch angle can be obtained by the sensor output and the calibration curve. Then the roll angle and yaw angle can be calculated using three-axis magnetic strength data. Inspired by Yu's work, Xiang et al. [10] used the three orthogonally mounted magnetometers to measure the attitude of high spinning projectile. Both of the above methods can calculate the attitude of the projectile using the magnetometer. Unfortunately, the attitude estimation accuracy was easily interfered by the external environment, and the calibration work was cumbersome. A roll cycle can only update the attitude information twice. Raul et al. [11] measured the attitude of the high spinning projectiles using a combination of accelerometer, GNSS, and semiactive image detector. In this solution, the acquisition accuracy of the gravity vector was closely related to the established dynamic model of the projectile, while the signal of GNSS might be lost, and the data update rate was low. Christophe et al. [12] estimated the angular velocity of the projectile using low-cost sensors (two-axis accelerometers and three-axis magnetometers), combined with a 6-degree-of-freedom ballistic model. In this method, the parameter setting of the 6-DOF ballistic model directly determined the estimation accuracy of the angular rate of the projectile, and the system was poor in robustness. The solutions proposed in [13–15] are similar to the above solutions. All of them suffer from poor autonomy, high complexity, and vulnerability to external environmental interference.
In order to solve the problem that the low-range gyroscope cannot measure the angular rate of high spinning projectiles, a compound angular motion model of high spinning projectiles is established by analyzing the projectile characteristics, and a roll angular rate estimation method (BSCZT-KF) for high spinning projectile is proposed. We first use the peak detection method to obtain the period and frequency of the transverse magnetometer signal, and adopt the least squares method to perform polynomial fitting on the rolling frequency. Then the frequency of transverse gyroscope signal is extracted using the improved CZT algorithm, BSCZT. Finally, the frequency fitting function and the BSCZT are fused by KF to estimate the roll angular rate of the projectile, which further improve the accuracy and the real-time performance. The hardware system of the scheme is simple in structure, strong in autonomy, low in algorithm time complexity, and high in data update rate, so it is extremely suitable for attitude measurement of high spinning projectiles.
2. The Angular Motion Model of High Spinning Projectile
The coordinate system of this paper mainly includes three coordinate systems, which are local geographic coordinate system, projectile coordinate system, and quasi-elastic system. Among them, the local geographic coordinate system (East-North-Sky) is selected as the navigation reference coordinate system. The axis points to east. The points to north. The points vertically to sky. The projectile coordinate system takes the center of mass of the projectile as the origin. The points to the head of the projectile. The is the radial axis. The is vertical to the plane consisting of and . The quasi-elastic system coincides with the projectile coordinate system at the initial moment. The difference between the quasi-elastic system and the projectile coordinate system is that the quasi-elastic system does not rotate with the projectile. In this paper, two-axis MEMS gyroscope and single-axis MEMS magnetometer are used to estimate the rolling angular rate. The sensitive axes of gyroscope coincide with and . The sensitive axis of magnetometer coincides with . It can be seen from Figure 1.
According to the literature [16], the angular motion of a spinning projectile can be decomposed into three types of motion: the self-rotation motion around the center of mass of the projectile, the coning motion around velocity vector, and the swing motion around radial axis.
2.1. The Self-Rotation Motion around the Center of Mass of the Projectile
According to the literature [16], the kinematic equation of the self-rotation motion around the center of mass of the projectile iswhere are the components in the direction of three axes of the projectile coordinate system. are pitch angle, yaw angle, and roll angle, respectively.
By solving (1), we can get the following.
Based on the auxiliary angle formula of the trigonometric function, the expressions of are modified as follows.
It can be seen from (3) that the Y and Z axes gyroscope signals are a set of sine and cosine curves. The amplitude is , and the frequency is . The phase difference between the two axes is 90°.
2.2. The Coning Motion around Velocity Vector
The rotating projectile can reduce the interference of the eccentricity, aerodynamic asymmetry, uneven mass distribution, and other factors on the stability of the projectile through its own rotation. When the rotation speed of the projectile is higher than a certain threshold, the projectile will not be reversed due to the turning moment caused by gravity, but will be rotated around the velocity vector according to a certain periodicity. The angular velocity vector of the classical cone motion in the projectile coordinate system iswhere is the angular velocity vector of classical cone motion. is the precession frequency. is the half-cone angle.
However, the coning motion around velocity vector of projectile is different from the classical cone motion. It is based on the classical cone motion and adds the high speed rotation along the axial direction of the projectile. Let the self-rotation speed be ; the attitude quaternion of the rotation process is as follows.
The attitude transformation matrix of rotation around is as follows.
Combining (4) and (6), we can obtain the angular motion model of rotation around velocity vector.
2.3. The Swing Motion around Radial Axis
Subjected to the external moment, the high spinning projectile will conduct the coning motion around the velocity vector and oscillate around the radial axis of the projectile. When the rotation speed of the projectile is large, the amplitude of the swing is small. When the rotation speed decreases, the amplitude of the swing becomes larger. In the dynamic stability analysis of high spinning projectiles, the amplitude of the oscillation is one of the important parameters. Let the angular rate of oscillation be and the amplitude of the oscillation be ; then the angle of oscillation is as follows.
Based on (8), we can obtain the angular velocity vector of oscillation in the quasi-elastic coordinate system.
By (6), we can get the angular motion model of the oscillation around the radial axis in the projectile coordinate system.
2.4. The Compound Angular Motion Model of High Spinning Projectile
Combining (3), (7), and (10), we can get the compound angular motion model of high spinning projectile.
Taking the Y-axis angular motion model as an example, its compound angular rate is as follows.
It can be seen from (12) that the Y-axis compound angular rate contains three sinusoidal components. Their frequencies are , , and , respectively. Similarly, the angular rate of the Z-axis is the same as the angular rate of the Y-axis, and it is also composed of several sinusoidal components. Therefore, the compound angular motion model of the radial axis of the high spinning projectile can be defined as follows:where denotes the number of sinusoidal components. are the amplitudes of angel oscillation. are the phases. is the angular motion frequency. In addition to the frequency of rotation, precession, and nutation, there are various frequencies of high-frequency microvibration in the actual high spinning projectile. Therefore, it is undoubtedly reasonable to use (13) to describe the angular motion of a high spinning projectile.
3. The Time-Frequency Analysis Method for Roll Angular Rate Estimation
According to the high spinning projectile compound angular motion model described in Section 2, the angular motion of the radial axis (Y-axis/Z-axis) contains multiple frequency components, including the frequency of rotation, precession, nutation, and other high-frequency vibrations. Therefore, the time-frequency analysis method can be used to analyze the gyro signal of the radial axis to obtain the respective frequency components, thereby obtaining the roll angular rate of the projectile. At this point, the roll angular rate estimation problem of the high spinning projectile is transformed into the frequency estimation problem of the radial axis gyro signal.
For the frequency estimation of the signal, the commonly used methods include peak detection, Short Time Fourier Transform (STFT) [17], Wavelet Transform (WT) [18], Wigner-Ville Distribution (WVD) [19], Hilbert-Huang Transform (HHT) [20], Chirp Transform (CT) [21], Polynomial Chirp Transform (PCT) [22], Zoom FFT [23], and Chirp-Z Transform (CZT) [24]. Among them, the peak detection method calculates the signal period by detecting the peak value of the signal, thereby obtaining the frequency of the signal. This method can only output the signal frequency once in a rolling cycle. It is also susceptible to noise interference, so the frequency estimation accuracy is low. If the estimated angular rate of the high spinning projectile is required to be within the error of 3°, the accuracy of the frequency estimation should reach 0.0083 Hz. Therefore, only the Zoom Fourier Transform (ZFFT) and the Chirp-Z Transform (CZT) can be used to realize high-precision estimation of signal frequency.
If the signal frequency is a constant value, both the Zoom Fourier Transform (ZFFT) and the Chirp-Z Transform (CZT) can be used to extract the signal frequency. However, the roll angular rate of the spinning projectile changes with time, the change speed is slow, and the change amplitude is small. Thence, the radial gyro output signal of the spinning projectile can be defined as a narrowband nonstationary signal. According to the description in [25], the CZT method has the best frequency estimation accuracy for narrowband nonstationary signals. Unfortunately, the CZT method also has obvious disadvantages. Frequency estimation accuracy is closely related to data length. For less collected data, the estimation accuracy will be reduced. Frequency resolution is closely related to the degree of frequency band subdivision. The denser the frequency band is, the higher the frequency resolution is, but the algorithm time complexity will also increase.
In order to overcome the above problems, this paper proposes a BSCZT-KF algorithm for estimating the angular rate of high spinning projectile. By combining the roll frequency fitting function with the BSCZT algorithm by KF, it is possible to use a shorter time series and coarser frequency band subdivision to achieve high-precision frequency estimation of signals at all times, and have a lower time complexity degree.
3.1. The BSCZT Algorithm
3.1.1. The Segmented CZT Algorithm
CZT transform can achieve frequency subdivision for part frequency bands in the entire signal band. With the same input data, a much higher frequency resolution than the FFT transform can be obtained for the frequency band of interest. However, the time complexity of CZT transform is higher. In particular, when the data sequence is long, it takes a lot of time to execute the CZT algorithm, which is not suitable for real-time analysis of the signal spectrum. In order to improve the efficiency of CZT transform and to process data flexibly, it is necessary to improve the CZT algorithm.
Assuming that the collected data sequence is , its Z Transform is as follows.
If we divide into L segments, each length is N, and then the i-th segment is as follows.
Equation (14) can be converted into the following.
The CZT transform of is denoted by . L segments have L transformations. Equation (16) can be rewritten as follows.
Let , , then we can obtain
The execution results of SCZT algorithm and CZT algorithm are the same, but the execution efficiency is very different. If the data length is and the length of is M, the algorithm time complexity of FFT, CZT and SCZT is, respectively, as follows.
It can be seen from (19) that the time complexity of the SCZT algorithm is smaller than that of the CZT algorithm, and the SCZT algorithm performs more efficiently.
3.1.2. The Blackman Segmented CZT Algorithm
In order to suppress the influence of spectrum leakage and fence effect on signal frequency estimation accuracy, it is necessary to select a suitable window function to process the gyroscope measurement signals. Commonly used window functions include Gaussian window, Hanning window, Chebyshev window, Hamming window, and Kaiser window. According to the literature [26], this paper chooses the Blackman window, which is often used to detect two signals with similar frequencies but different amplitudes. The frequency of rotation, precession, and nutation components in the radial gyro signal are relatively close, so the Blackman window can be used to better distinguish different frequency components.
Assuming that a measurement window containing M consecutive sampling points is
and Blackman window of length M iswhere
The windowing process is realized by multiplying a by b. Then we can use SCZT algorithm to extract the frequency of , and make this frequency as the instantaneous frequency in the middle of the sequence. The Blackman Segmented CZT algorithm is as follows:where denotes the Blackman window sequence of length N.
3.2. The BSCZT-KF Algorithm
3.2.1. The Selection of Subdivision Band
From (12), it can be concluded that the angular rate signal measured by the Y-axis gyroscope contains three main frequency components, , and . When using the BSCZT algorithm to perform time-frequency analysis on the Y-axis gyro signal, it is necessary to set the subdivision frequency band. The setting of the starting frequency, the cutoff frequency, and the degree of refinement of the refinement band directly determines the accuracy of the frequency estimation. At the same time, in order to accurately obtain the projectile rolling frequency, the selected subdivision band should contain , but does not contain and . Therefore, it is necessary to estimate the approximate range of the rolling frequency in advance.
First, as can be seen from the three frequency expressions, the body roll frequency is in the middle of the three frequencies. Secondly, we can use the peak detection method to make a rough estimate of the projectile rolling frequency, and then set the appropriate subdivision band. The peak detection method usually uses the magnetic field data collected by the magnetometer installed radially on the projectile.
The geomagnetic field, like the gravity field, is an inherent feature of the earth. Generally, geomagnetic seven elements are used to represent geomagnetic information at a certain point on the earth. The seven elements of geomagnetism are magnetic declination (D), magnetic inclination (I), horizontal intensity (), north component (), east component (), vertical component (), and total field (M). The relationship between them is shown in (24).
The three attitude angles that define the attitude of the missile are yaw (), pitch (), and roll (). According to the three attitude angles, the rotation matrix of the projectile can be obtained as shown in (25).
According to (25), the relationship between the magnetic field strength in the geomagnetic coordinate system and the one in the projectile coordinate system is as shown in (26).
According to (24)–(26), we can get the following.
By sorting out (27), we can obtain
where the magnetic inclination and total field are constant within the range of the projectile.
The following transformation is performed on the theoretical output of the Y-axis.
According to the auxiliary angle formula of trigonometric function, we can get the following.
According to (30), the output of the Y-axis magnetometer is a sinusoid with amplitude and frequency . Therefore, the peak detection method can be used to obtain the period of the output of the Y-axis magnetometer, thereby obtaining the rolling frequency of the projectile and thus achieving a rough estimation of the rolling frequency of the projectile. Then, we can set the range of the subdivision band according to the rough estimate of the projectile roll frequency.
3.2.2. Roll Frequency Fitting
The subdivision band is selected by using the peak detection method and the radial magnetic strength data to estimate the rolling frequency of the projectile. For subsequent sampled data, we can also use the peak detection method to estimate the projectile roll frequency. However, the peak detection method can only output the frequency at one point, and cannot get the frequency of the intermediate moment of the roll. According to experience, the rolling frequency of the projectile changes slowly. We can use the least squares method to fit the known rolling frequency to get the function of the rolling frequency with respect to time, so that the frequency value of each moment can be obtained. The roll frequency fitting function is as follows:
where denotes the roll frequency. denote polynomial coefficients, which are obtained by least squares fitting.
It should be pointed out that for the case where the rolling frequency is severe, the rolling frequency should adopt the method of segment fitting. For the case where the rolling frequency changes gently, the global fitting method can be adopted. The fitting coefficients of the roll frequency need to be updated after every roll cycle.
3.2.3. The Design of Kalman Filter
According to the above content, the radial axis magnetometer signal can be used to estimate the angular rate of the rolling angle of the projectile, then the subdivision band and subdivision accuracy of the BSCZT method are set, and the accurate projectile roll frequency is obtained by the BSCZT algorithm. However, there still exist the following problems. Firstly, the algorithm time complexity is still high. Although the segmentation CZT algorithm (SCZT) is adopted, the time complexity of the algorithm can be appropriately reduced, but the time complexity of the SCZT algorithm is also related to the value of the frequency subdivision M. In order to obtain a higher frequency resolution, the value of the subdivision degree M is generally large, so that the time complexity of the SCZT algorithm becomes higher. Secondly, frequency estimation has a delay. The frequency obtained by BSCZT is the frequency at the moment in the middle of the window, so the frequency estimation has a delay of N/2 (N is the window length). For the last N/2 data of the sequence, it is not even possible to estimate the frequency.
In order to solve the above problem, that is, the algorithm has high time complexity and frequency estimation delay, this paper combines the BSCZT algorithm with the obtained frequency fitting function by KF, which not only improves the real-time performance of frequency estimation, but also can improve frequency estimation accuracy based on sparsely subdivision band. The BSCZT-KF algorithm mainly includes three parts: BSCZT module, roll frequency fitting module, and KF fusion module. The KF process is mainly divided into a status update process and a measurement update process, which are described as follows.
(1) The State Equation. The selected state vectors are , where is the roll frequency at the moment k. are the 5th-order fitting coefficients of the changing rate of rolling frequency. The 5th-order fitting coefficients are obtained by peak detection method. The state equation is as follows:
where F denotes the state transition matrix. G denotes the noise gain matrix. W denotes the noise matrix. The values of F and G are as follows.
(2) The Measurement Equation. The roll frequency of the projectile is selected as the measurement, which can be obtained by the BSCZT algorithm. The measurement equation is as follows:
where V is the measurement noise matrix. H is the measurement matrix, which is expressed as follows.
3.2.4. The Flow of BSCZT-KF Algorithm
The flowchart of BSCZT-KF algorithm is illustrated in Figure 2.
The BSCZT-KF algorithm has four steps.
Step 1 (projectile roll frequency fitting). The peak detection method is used to process the first 2 seconds of collected data, and several groups of rolling frequency of projectile and the moment of occurrence are obtained. The data length selected here can be changed according to the actual sampling frequency and the projectile roll frequency. Then we can get the function of the rolling frequency of the projectile relative to time using the least squares polynomial curve fitting. After the above operation, on the one hand, a rough estimate of the roll frequency can be obtained for the subdivision band selection of the BSCZT algorithm; on the other hand, the initial value of the state variable in the KF state equation can be obtained.
Step 2 (fine estimation of the projectile rolling frequency). The frequency of radial gyro signal can be extracted by BSCZT algorithm, then we can obtain the roll frequency of projectile at the moment . denotes the current sampling moment. denotes the sampling period. N is the length of Blackman window. It should be noted that before using the BSCZT algorithm, it is necessary to set the range of the subdivision band and the subdivision degree according to the rough estimation of the projectile roll frequency.
Step 3 (KF process). At the moment , we first update the state equation of KF. Then we make the frequency estimated by BSCZT as the measurement to update the measurement process. Finally, we can get the best estimation of roll frequency of projectile after the complete KF process.
Step 4. If the new roll cycle arrives, the function of roll frequency should be refitted using the new frequency obtained by the peak detection method. Then, skip to Step 2.
According to the above steps, an accurate estimation of the rolling frequency of the projectile can be achieved. In terms of the shortcomings of the above BSCZT algorithm, the BSCZT-KF algorithm can be combined with the rolling frequency fitting on the basis of the rough frequency band division to obtain a more precise estimation of frequency, and the algorithm time complexity is reduced. If it is necessary to obtain a real-time frequency estimation, we can directly calculate the frequency from the roll frequency fitting function. The accuracy of the frequency fitting function coefficients is guaranteed by the KF process.
4. Experiment and Analysis
In order to verify the validity and frequency estimation accuracy of the proposed BSCZT-KF algorithm, the two sets of real flight data are processed according to the above algorithm steps, and the roll frequency and roll angular rate of the projectile are obtained. It should be pointed out that the first set data is from a certain type of 120mm projectile, which spins with a speed of 2r/s, and the second set data is from a certain type of 130mm projectile, which spins with a speed of 5r/s. For the purpose of verifying the accuracy of the BSCZT-KF algorithm, the angular rate of the X-axis (rolling axis) measured by the MEMS gyroscope is selected as a reference. The relationship between the angular rate of the X-axis measurement of the on-load gyroscope and the roll angular rate is as follows.
Equation (37) is derived from the kinematic equation of the rotation around the center of mass, where denotes the angular rate of the X-axis measurement of the on-load gyroscope. denotes the roll angular rate of projectile. denotes the changing rate of yaw angle. denotes the pitch angle. Because the changing rate of yaw angle is relatively small during the actual flight, (37) can be written as follows.
Therefore, it can be considered that the X-axis angular velocity measured by the projectile-mounted gyroscope is the rolling angular rate of the projectile. Obviously, the gyroscope X-axis measurement data can be used to verify the frequency estimation accuracy of the BSCZT-KF algorithm.
4.1. 120mm Projectile Rolling Angular Rate Estimation
The actual flight data of the 120mm projectile is obtained by the three-axis gyroscope and three-axis magnetometer. The sample frequency is 1000Hz. The mounting position of the two sensors coincides with the projectile coordinate system. Since only the data on the outer ballistics of the projectile is of interest, the data in the internal ballistic process is truncated. External ballistic data of three-axis gyroscopes and three-axis magnetometers are shown in Figures 3 and 4.
4.1.1. Projectile Roll Frequency Fitting
As can be seen from Figure 4, the magnetometer can sense the rolling motion of the projectile, and the period of the radial axis magnetic strength data is the same as the projectile rolling period. Therefore, the peak detection method can be used to estimate the rolling frequency of the projectile. According to the obtained frequency sequence, the function of the rolling frequency is obtained by the least squares polynomial fitting method.
Firstly, the Y-axis magnetic intensity data is filtered by a low-pass filter, and then the peak value of the Y-axis magnetic intensity signal and its time are obtained by the peak detection method. The roll period of the projectile can be obtained by the difference between adjacent peak moments. According to the rolling period, the rolling frequency of the projectile can be obtained. The roll frequencies are shown in Table 1.
The rolling frequency is converted into the rolling angular rate, and the function of the angular rate of projectiles is obtained by the least squares polynomial curve fitting method. The curves can be seen in Figure 5.
In Figure 5, the black line represents the X-axis angular rate actually measured by the gyroscope, and the red line represents the roll angular rate fitting data. It can be seen that the roll angular rate fitting function is close to the actual measured X-axis angular rate, but there is still a large error in the details.
4.1.2. The Extraction of Roll Frequency by BSCZT
The Y-axis gyroscope signal contains the rolling frequency component, so the time-frequency analysis of the Y-axis gyroscope signal can be used to obtain the rolling frequency. Here, the Y-axis gyro data is processed by the BSCZT algorithm proposed in this paper. First, we use the FFT algorithm to analyze the spectrum of the Y-axis gyroscope signal, as shown in Figure 6.
It can be seen from Figure 6 that there are three frequency components in the Y-axis gyro signal, which are consistent with the established compound angular motion model of the rolling projectile. Among them, the rolling frequency of the projectile is about 2Hz, which is in the middle of the three frequency components. This point once again confirms the correctness of the compound angular motion model of the rolling projectile established in this paper.
In the actual flight process, the projectile rolling frequency changes around 2 Hz, so only the frequency band around 2 Hz needs to be subdivided to obtain accurate projectile rolling frequency. The subdivision band range is set to 1.5~2.5Hz, and the subdivision accuracy M is set to 1000. The frequency resolution can reach 0.001Hz, and the ideal roll angular rate accuracy can reach 0.36°.
In order to improve the frequency estimation accuracy of the BSCZT algorithm as much as possible, the Blackman window length is set to 2000, and the time delay is 1s. In the CZT transformation of a window data sequence, in order to reduce the computational time complexity, the SCZT algorithm is used for CZT transformation, the number of segments is 5, and the length of each segment is 400. The frequency estimation result obtained by the BSCZT algorithm is shown in Figure 7.
It can be seen from Figure 7 that the projectile rolling frequency estimated by the BSCZT algorithm is very close to the actually measured rolling frequency. It is calculated that the average error between them is 0.0053 Hz, and the corresponding roll angular rate error is 1.908°, which can fully meet the required roll angular rate estimation accuracy.
4.1.3. KF Process
After the high-order fitting of the rolling frequency of the projectile and the rolling frequency estimation based on the BSCZT algorithm, the KF state update process and the measurement update process can be performed. Since the frequency estimation of the BSCZT algorithm has a time delay, the KF process is a guarantee of real-time frequency estimation.
Firstly, based on the high-order fitting function of the rolling frequency of the projectile, the state equation is established according to the content of Section 3.2.3, and the state equation is continuously updated with the passage of time. When the BSCZT algorithm is used to obtain the fine estimation of the roll frequency, the measurement update of the KF process is performed, and the coefficients of the high-order fitting function are corrected. In addition, each time a roll cycle is passed, the projectile roll frequency function is refitted.
Finally, after the above steps, real-time high-precision estimation results of the rolling frequency of the projectile can be obtained.
We used STFT, Zoom FFT, and BSCZT-KF algorithms to estimate the roll angular rates, respectively. Their results are as shown in Figure 8.
After being calculated, the average angular rate error between the roll angular rate measured by the X-axis gyroscope and the roll angular rate estimated by BSCZT-KF algorithm is 0.8782°, which is 0.121% of the maximum roll angular rate. The average angular rate error by STFT is 13.5765°, which is 1.89% of the maximum roll angular rate. The average angular rate error by Zoom FFT is 2.7615°, which is 0.38% of the maximum roll angular rate. Therefore, the BSCZT-KF algorithm has the highest estimation accuracy of roll angular rate.
4.2. 130mm Projectile Rolling Angular Rate Estimation
The above experiment is performed with 2r/s projectile flight data. In order to further verify the validity and accuracy of the proposed algorithm and its applicability to high spinning projectile, this section will use the above method to analyze the flight data of projectile with the rotation speed of about 5r/s.
The installation position of the three-axis gyroscope is the same as the above content, and the collected raw data of the three-axis gyroscope is as shown in Figure 9.
As can be seen from Figure 9, the actual roll angular rate of the projectile varies from 1400° to 2300°. First, the spectrum analysis of the Y-axis gyroscope output is performed by FFT. Since the range of roll angular rate changes is relatively large, a piece of data with a slow frequency change is selected for spectrum analysis, such as data between 30 seconds and 40 seconds. The spectrum obtained by FFT is shown in Figure 10.
As can be seen from Figure 10, the Y-axis gyroscope signal mainly contains three frequency components, and the frequency in the middle position is about 4 Hz. Combined with the roll angular rate (1500°r/s) at 3 seconds to 4 seconds in Figure 9, the projectile roll frequency is equal to the intermediate frequency of the three frequency components of the Y-axis gyro signal. In addition, Figure 10 once again proves the correctness of the compound angular motion model of the rolling projectile established in this paper.
The Y-axis gyroscope data is processed by BSCZT-KF to obtain the rolling frequency and roll angular rate of the projectile. The detailed processing will not be described again. Please refer to the above content. The rolling frequency estimated by the BSCZT algorithm is shown in Figure 11, wherein the selected subdivision band is 3 Hz to 7 Hz, and the band subdivision accuracy is selected to be 5000. Since this set of gyro data is long, the Blackman window length is set to 4000 and the frequency estimation delay is 2 seconds.
As can be seen from Figure 11, the BSCZT algorithm can better estimate the rolling angular frequency of the projectile in the part where the rolling frequency of the projectile changes slowly, but the frequency estimation error is larger in the part where the frequency changes faster. At the same time, the method has a time delay of 2 seconds, and the frequency cannot be estimated for the last 2 seconds of the gyro data.
According to the above algorithm steps, the proposed BSCZT-KF algorithm is used to analyze the Y-axis gyro, and the obtained body rolling angular rate is as shown in Figure 12.
As can be seen from Figure 12, the error of the roll angular rate estimated by the BSCZT-KF algorithm is larger than that of the BSCZT algorithm. The reason is that when the Y-axis magnetometer data is used for high-order fitting of the roll angular rate, all the data before the current time are used for fitting. For the case where the overall angular rate of roll is slow, the use of all historical data fitting has little effect on the fitting effect, but for the case where the roll angular rate changes drastically (as in the 20 seconds in Figure 9), using all the historical data to fit, and the order of the fitted function remains the same, the fitting accuracy will be seriously degraded. Therefore, this paper uses a piecewise high-order fitting method to obtain the function of the rolling frequency of the projectile. The resulting projectile rolling angular rate is shown in Figure 13.
The results show that after changing the overall fitting in the BSCZT-KF algorithm to the piecewise fitting, the estimation accuracy of the rolling angular rate of the projectile is significantly improved, wherein the estimated average error is 1.9417°, which is about 0.095% of actual roll angular rate.
4.3. Algorithm Performance Analysis
The above two experiments can prove that the BSCZT algorithm proposed in this paper can use the output of the radial sensor to estimate the rolling frequency and the rolling angular rate of the spinning projectile with high precision, which provides an efficient solution for the measurement of the rolling angular rate of the high spinning projectile. It can be seen from the algorithm steps and parameter settings that the accuracy and complexity of the algorithm are affected by many factors. The following three main factors are discussed.
4.3.1. The Length of Window
In the signal frequency extraction using the BSCZT algorithm, in order to suppress the effects of spectrum leakage and fence effect, window processing is usually adopted. However, the choice of window length will affect the frequency extraction accuracy and real-time performance of the algorithm.
The longer the window length, the higher the frequency estimation accuracy, but the higher the algorithm time complexity.
The window length must not be less than four roll cycles; otherwise the algorithm will fail.
The existence of the sliding window causes the frequency estimation to have a delay of N/2. The longer the window, the larger the delay of the frequency estimation, and the worse the real-time performance of the algorithm.
Therefore, it is necessary to set a reasonable window length according to the rolling frequency of the projectile, the data sampling rate, and the real-time requirements of the algorithm.
4.3.2. Frequency Band Subdivision
In the BSCZT algorithm, we need to set the subdivision of the frequency band. In theory, the denser the frequency band subdivision, the higher the frequency resolution, and the higher the frequency estimation accuracy. However, according to the specific experimental results, it is found that when the frequency band subdivision reaches a certain intensity, the frequency subdivision is further improved, and the frequency estimation accuracy is no longer improved. On the contrary, it brings a huge computational burden and improves the time complexity of the algorithm. Therefore, in practical applications, a reasonable degree of frequency band subdivision should be set.
The sampling frequency of the data in this paper is 1000Hz, and the selected frequency band is around 5Hz. Therefore, the frequency subdivision parameter M is set to 2000 to meet the accuracy requirement of frequency estimation. If the frequency subdivision is continued, the accuracy of the frequency estimation will not change, and the time complexity of the algorithm will become higher.
4.3.3. The Roll Frequency of Projectile
In the above actual flight data analysis experiments of 120mm and 130mm projectiles, the roll angle rates of the projectile are not more than 7 r / s. For the case of high spinning projectile, such as 100~200 r/s, because it is difficult to obtain actual flight data, the numerical simulation method will be adopted to verify the feasibility of the BSCZT-KF algorithm.
The simulation condition is as follows. The simulation time is 6 seconds and the data sampling frequency is 1000Hz. The initial yaw angle of the projectile is 60° and keeps still during the flight. The initial pitch angle is 30°, and the rate of change during flight is -10°. The initial roll angle is 0 and the roll frequency is as follows.
Under the above simulation conditions, the local curve of the rolling angle of the projectile is shown in Figure 14.
According to the change of the attitude angle of the projectile and the kinematic equation of the rotation around the center of mass, the generated roll axis gyroscope data is shown in Figure 15.
According to the IGRF11 geomagnetic model, the geomagnetic intensity is 56864nT at the place whose longitude is 122.9°E, latitude is 47.47°N, altitude is 230 meters, and the magnetic dip is 65.89°. For the above-mentioned change of the attitude of the projectile, the three-axis magnetometer output is shown in Figure 16.
By using the Y-axis magnetometer data, a rough estimation of the roll angular rate can be obtained by the peak detection method. Since the rolling frequency of the projectile fluctuates around 150 Hz, the average rolling period is 6.67 ms, and there are only 7 sampling points in one rolling cycle. Therefore, the detected peak moment error is large, and the final calculated rolling cycle and rolling frequency error are also correspondingly larger. The relationship between the error of the peak detection method and the rolling frequency of the projectile is shown in Table 2.
As can be seen from Table 2, for the high spinning projectiles, due to the limitation of the data sampling rate, the error of the rolling frequency of the projectile estimated by the peak detection method is large. The rolling frequency fitting error in the BSCZT-KF algorithm proposed in this paper also increases with the increase of the rolling frequency. Therefore, it is necessary to adjust the noise matrix of the state equation according to the rolling frequency.
It is worth noting that the BSCZT algorithm can be used to estimate the rolling frequency of the projectile, which can achieve higher estimation accuracy. For the rolling frequency in Table 2, the error of the rolling frequency estimated by the BSCZT algorithm is about 0.0154 Hz. Under the above numerical simulation conditions, the roll frequency estimation accuracy of the BSCZT algorithm is independent of the rolling frequency of the projectile, which proves that the BSCZT algorithm is extremely suitable for the roll angular rate estimation of high spinning projectile.
For the case of 250 Hz roll frequency, the simulation time is 6 seconds and the data sampling rate is 1000 Hz. The initial yaw angle of the projectile is 60° and keeps still during the flight. The initial pitch angle is 30°, and the rate of change during flight is -10°. The initial roll angle is 0 and the roll frequency is as follows.
Considering that the data sampling rate 1000 Hz is a little slow for the case of 250 Hz roll frequency, we set the data sampling rate as 2k Hz, 10k Hz, 20k Hz, and 50k Hz, separately. Then we estimate the roll frequency by the proposed BSCZT-KF method. The results of frequency estimation are shown in Table 3 and Figure 17.
As can be seen from Table 1 and Figure 1, the BSCZT-KF algorithm can accurately estimate the roll frequency at 250 Hz case. The estimation error is about 0.012% of real roll frequency. What is more, with the increasing of sampling rate, the accuracy of estimation is enhanced. However, the improvement in accuracy is not very obvious. From another perspective, the BSCZT-KF method has no special requirements for sampling rate. As higher sampling rates place higher demands on the hardware, it is important to choose proper sampling rate and determine satisfied estimation accuracy.
5. Conclusions
Accurate measurement of the attitude of the projectile is the key technology of navigation and guidance. For high spinning projectiles, the measurement of roll angular rate has always been a challenge. Although various techniques can be used to replace the gyroscope for high spinning projectiles roll rate measurement, the measurement accuracy, real-time performance, and system stability are difficult to guarantee. Aiming at the defects of the prior art, this paper establishes a compound angular motion model of high spinning projectiles, by studying its motion characteristics, and proposes a time-frequency analysis method BSCZT-KF algorithm for accurate estimation of the rolling frequency of high spinning projectiles in real time. In order to verify the validity and accuracy of the proposed algorithm, two sets of actual flight data are processed. The experimental results show that the compound angular motion model of the rolling project proposed in this paper is corrected and prove that the BSCZT-KF algorithm can adopt the radial axis gyro data to efficiently and accurately estimate the rolling frequency of the spinning projectile. Finally, the important factors of the algorithm are discussed. The discussion shows that the parameters in the algorithm should be set according to the rolling frequency and data sampling rate. It is worth noting that the accuracy of the BSCZT algorithm for the frequency estimation of the rolling projectile is independent of the rolling frequency of the projectile, indicating that the algorithm is extremely suitable for attitude measurement of high spinning projectile. For the subsequent research work, it will expand on the real-time performance of the algorithm.
Data Availability
The 120mm projectile and 130mm projectile measurement data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare no conflicts of interest.
Acknowledgments
This work was supported by National Natural Science Foundation (NNSF) of China under Grant 61771059 and 61471046 and Beijing Natural Science Foundation under Grant 4172022.