Separation and feature extraction of micro‐motion signal of ballistic target

The Journal of EngineeringVolume 2021, Issue 12 p. 828-837 ORIGINAL RESEARCH PAPEROpen Access Separation and feature extraction of micro-motion signal of ballistic target Yuxi Li, Corresponding Author Yuxi Li 2821098869@qq.com orcid.org/0000-0002-7358-441X Graduate College, Air Force Engineering University, Xi'an, Shaanxi, China Correspondence Li Yuxi, Postal: No.1 Jiazi, Changle East Road, Baqiao District, Xi'an City, Shaanxi Province. Email: 2821098869@qq.comSearch for more papers by this authorCunqian Feng, Cunqian Feng Air and Missile Defense College, Air Force Engineering University, Xi'an, Shaanxi, ChinaSearch for more papers by this authorXuguang Xu, Xuguang Xu Graduate College, Air Force Engineering University, Xi'an, Shaanxi, ChinaSearch for more papers by this authorLixun Han, Lixun Han orcid.org/0000-0002-0449-1957 Graduate College, Air Force Engineering University, Xi'an, Shaanxi, ChinaSearch for more papers by this authorDayan Wang, Dayan Wang Graduate College, Air Force Engineering University, Xi'an, Shaanxi, ChinaSearch for more papers by this author Yuxi Li, Corresponding Author Yuxi Li 2821098869@qq.com orcid.org/0000-0002-7358-441X Graduate College, Air Force Engineering University, Xi'an, Shaanxi, China Correspondence Li Yuxi, Postal: No.1 Jiazi, Changle East Road, Baqiao District, Xi'an City, Shaanxi Province. Email: 2821098869@qq.comSearch for more papers by this authorCunqian Feng, Cunqian Feng Air and Missile Defense College, Air Force Engineering University, Xi'an, Shaanxi, ChinaSearch for more papers by this authorXuguang Xu, Xuguang Xu Graduate College, Air Force Engineering University, Xi'an, Shaanxi, ChinaSearch for more papers by this authorLixun Han, Lixun Han orcid.org/0000-0002-0449-1957 Graduate College, Air Force Engineering University, Xi'an, Shaanxi, ChinaSearch for more papers by this authorDayan Wang, Dayan Wang Graduate College, Air Force Engineering University, Xi'an, Shaanxi, ChinaSearch for more papers by this author First published: 22 September 2021 https://doi.org/10.1049/tje2.12083AboutSectionsPDF ToolsRequest permissionExport citationAdd to favoritesTrack citation ShareShare Give accessShare full text accessShare full-text accessPlease review our Terms and Conditions of Use and check box below to share full-text version of article.I have read and accept the Wiley Online Library Terms and Conditions of UseShareable LinkUse the link below to share a full-text version of this article with your friends and colleagues. Learn more.Copy URL Share a linkShare onFacebookTwitterLinked InRedditWechat Abstract During the process of the ballistic target's mid-flight, it is very important to accurately identify the target. The separation of target micro-Doppler curves and the extraction of characteristic parameters are the key to accurate identification. Aiming at the slow convergence speed of the traditional Fast-ICA algorithm, this paper proposes an improved Fast-ICA algorithm, which realizes the separation of the micro-Doppler curve of the scattering points. On this basis, the trajectory of scattering points in space is analyzed, the expression of the target feature value is established, and the target feature parameters are effectively extracted. According to the results of the simulation experiment, the algorithm can achieve better signal separation and the characteristic parameters of the target are effectively extracted. 1 INTRODUCTION Nowadays, due to the increasing threat of ballistic missiles, all countries are stepping up research and building ballistic missile defence systems [1]. Ballistic target recognition technology is one of the core technical problems that need to be solved in ballistic missile defence. With the rapid development of multiple independently re-entry vehicles and missile decoy technology [2], during the flight, ballistic missiles will release fake warheads and decoys, and at the same time, they will also produce many fragments, which will interfere with the radar and improve the survivability of warheads. Therefore, the ballistic target recognition technology based on traditional feature quantities can no longer meet the needs of modern high-tech warfare. As the inherent properties of ballistic targets, micro motion characteristics are difficult to imitate, and there are obvious differences in the movement forms of warheads, decoys, and fragments. Therefore, they are used to distinguish and identify real warheads [3]. In recent years, they have received extensive attention from domestic and foreign research scholars. The extraction of relevant parameters of ballistic missiles and the identification of true and false warheads are the key to the establishment of ballistic missile defence systems. Among them, the length information and precession characteristics of the target are important for distinguishing warheads and decoys [4]. Since precession will produce micro-Doppler modulation of radar echo, the micro-Doppler of radar echo can be analyzed, and then the characteristic parameters of the target can be obtained, which are used to identify [5-8]. The literature [9] uses the generalized Radon transform method to transform the range profile information to obtain the precession angle. On this basis, the micro-Doppler joint equations are constructed, and the related parameters are solved. The literature [10] separates and extracts the micro-Doppler information according to the change trend of the micro-Doppler information change rate of the scattering centre. Then, the non-linear least squares fitting method is used to obtain the micro-Doppler information to achieve the target reconstruction. The literature [11] can determine the echo signal form of the scattering point according to the target micro-motion type, and at the same time, construct the atom library, and use the ROMP method to extract the parameters. The literature [12] uses the non-linear least squares estimation method to calculate the amplitude-equal correlation parameters of the scattering centre, and achieves signal separation through correlation processing. On this basis, the micro-motion characteristics of the ballistic target are used to estimate the target's precession characteristics and structural parameters. The literature [13] determined the evaluation criteria for the optimal design of the curve aperture according to the principle of irrelevance, and realized the target three-dimensional feature extraction by using the basis tracking method based on global optimization. In order to detect the micro-Doppler signal induced by the helicopter blades, a real-time passive coherent location (PCL) system based on digital video broadcasting (DVB) is proposed in the literature [14]. This method is more likely to adapt to urban environments. With the time-frequency characteristics of the micro-Doppler component signal, the literature [15] proposes a multi-order STFRFT time-frequency algorithm, which effectively reduces the amount of calculation and improves the probability of the target identification. Based on the above research status, this article first determines the type of scattering points. On this basis, the Fast-ICA algorithm is used to separate the obtained echo signals. Aiming at the slow convergence speed of the traditional Fast-ICA algorithm, the iterative formula is improved to increase the convergence speed. Under the condition of a monostatic radar, only the structure and movement characteristics of the target's micro-movement feature in the direction of the radar's line of sight can be extracted. For the target recognition of the warhead, although it is easy to obtain the rotation period information of the warhead target, it is difficult to obtain the true bottom diameter. Under different radar perspectives, due to the complexity of target attitude changes, its micro-Doppler characteristics will show significant differences, which will affect the accuracy of target recognition. In order to overcome the attitude sensitivity of the target's micro-Doppler feature, this paper uses the networked radar to extract the target feature parameters. By deriving the relevant expressions of the characteristic parameters, the radar micro-Doppler information joint equation group is constructed to solve the relevant parameters and obtain the relevant characteristic parameters of the target. 2 PRECESSION MODEL OF BALLISTIC TARGET In order to simplify the analysis, this paper takes the tailless cone ballistic target as the research object. The centre of rotation of the ballistic target is O. As shown in Figure 1 the reference coordinate system O-XYZ is established with O as the origin. The ballistic target takes the axis Z of the coordinate system as the precession axis of the cone motion, and its cone rotation angular velocity is ω c . It spins on its own axis of symmetry, and its spin angular velocity is ω s . The angle between the cone axis and the spin axis, that is, the precession angle is θ. At the beginning, the azimuth and elevation angles of the radar line of sight LOS in the reference coordinate system O-XYZ are α and β respectively. The body coordinate system of the ballistic target is O-xyz, rotates ϕ e around the z-axis, rotates θ e around the x-axis, and rotates φ e around the z-axis to form a reference coordinate system. The distance from the centre of rotation to the apex A of the cone is l 1 , and the distance from the centre of the bottom of the cone is l 2 . FIGURE 1Open in figure viewerPowerPoint Scattering model of precession cone ballistic target The initial position of the scattering point on the target in the body coordinate system is r 0 = ( r x 0 , r y 0 , r z 0 ) T , then the position vector of the scattering point in the reference coordinate system at time t is r = R c R s R i n i t r 0 (1)where R c and R s are the spin and cone spin matrices of the ballistic target, respectively, R i n i t is the Euler rotation matrix, which is determined by the initial Euler angle ( ϕ e , θ e , φ e ) . According to the literature [16] and Rodrigues formula [17], the expression of the correlation matrix can be obtained as follows: R c = I + ⌢ ω c ′ sin ( Ω c t ) + ⌢ ω c ′ 2 ( 1 − cos ( Ω c t ) ) (2)where I is the identity matrix, ω c = ( ω cX , ω cY , ω cZ ) T is the cone angular velocity vector, Ω c = ∥ ω c ∥ , ω c ′ = ω c / Ω c , ⌢ ω c ′ are the antisymmetric matrix of ω c ′ : R s = I + ⌢ ω s ′ sin ( Ω s t ) + ⌢ ω s ′ 2 ( 1 − cos ( Ω s t ) ) (3)where ω s = ( ω sX , ω sY , ω sZ ) T is the spin angular velocity vector, Ω s = ∥ ω s ∥ , ω s ′ = ω s / Ω s , ⌢ ω s ′ are the antisymmetric matrix of ω s ′ . Construct the following oblique symmetric matrix [18] ⌢ ω c ′ = 0 − ω cZ ′ ω cY ′ ω cZ ′ 0 − ω cX ′ − ω cY ′ ω cX ′ 0 ⌢ ω s ′ = 0 − ω sZ ′ ω sY ′ ω sZ ′ 0 − ω sX ′ − ω sY ′ ω sX ′ 0 (4) R init = cos ϕ e − sin ϕ e 0 sin ϕ e cos ϕ e 0 0 0 1 1 0 0 0 cos θ e − sin θ e 0 sin θ e cos θ e × cos φ e − sin φ e 0 sin φ e cos φ e 0 0 0 1 (5) Let P be any scattering point on the cone warhead, and the distance between the scattering point P and the radar at time t is R ( t ) = R 0 + v t + r = n LOS T • ( R 0 + v t + R c R s R init r 0 ) (6)where R 0 is the position vector of the warhead's rotation centre from the radar, and n L O S is the direction of the radar's line of sight. According to formula (6), the echo signal of the target can be obtained. First, the signal is processed by the fast Fourier transform, then the envelope is removed, and finally a high-resolution range image can be obtained [16] S d ( f , t m ) = ∑ l = 1 L σ l T p sin c T p f + 2 μ c R Δ ( t m ) • exp − j 4 π f c c R Δ ( t m ) (7)where T p is the pulse width, μ is the modulation frequency, f c is the radar carrier frequency, and R Δ ( t m ) is the radial distance between the scattering centre and the rotation centre, L is the total number of scattering points of the target, and sinc ( x ) is the symplectic function. 3 OBTAINING MICRO-MOTION INFORMATION 3.1 Determination of the type of scattering point According to the literature [19], the specific expression of the distance between scattering point P and radar is as follows: R ( t ) = ∑ i = 1 4 A i sin ( ω ( i ) t + φ i ) + R 0 (8)where A i is the modulation coefficient, ω = ( ω c + ω s , ω c − ω s , ω c , ω s ) is the frequency information, φ i is the phase of the micro-Doppler information, and R 0 is a constant. The specific expression of A i and φ i can be obtained according to the literature [19]. According to formula (8), the distance between the scattering point and the radar can be obtained, and the distance formula is derived n times to obtain the nth order micro-motion distance of the scattering point, as shown in the following formula: R ( n ) ( t ) = ∑ i = 1 4 A i ω n ( i ) cos ω ( i ) t + φ i + π 2 ( n mod 4 ) (9)where mod is the remainder operator. The ratio of the nth order micro-motion distance of different scattering points is shown in the following formula: R 1 ( n ) ( t + Δ t ) R 2 ( n ) ( t ) = ∑ i = 1 4 A 1 i ω n ( i ) cos ω ( i ) ( t + Δ t ) + φ 1 i + π 2 ( n mod 4 ) ∑ i = 1 4 A 2 i ω n ( i ) cos ω ( i ) t + φ 2 i + π 2 ( n mod 4 ) (10) when Δ t → 0 , the above formula can be simplified to R 1 ( n ) ( t ) R 2 ( n ) ( t ) = ∑ i = 1 4 A 1 i ω n ( i ) cos ω ( i ) t + φ 1 i + π 2 ( n mod 4 ) ∑ i = 1 4 A 2 i ω n ( i ) cos ω ( i ) t + φ 2 i + π 2 ( n mod 4 ) (11) Since the micro-motion distance change trend of different types of scattering points is related to different parameters, and the micro-motion distance change trend of the same type of scattering points is related to the same parameter. Therefore, the changing trends of the nth order micro-motion distance of different types of scattering points over time are quite different, while the same type of scattering points show relatively similar change characteristics. Furthermore, by comparing the changing trends of the micro-Doppler curves of different scattering centres, the classification and marking of the micro-Doppler curves of different types of scattering centres can be realized [10]. 3.2 Separation of micro-Doppler curves According to the previous section, different types of scattering points can be determined. Next, the curves need to be separated in order to extract the relevant parameters of the scattering points. This article introduces the Fast independent component analysis algorithm to separate the signals. As shown in Figure 2 the calculation process of independent component analysis: Assuming that the components in the source s ( t ) are independent of each other, from observation x ( t ) , the components are separated by the decomposition system B, so that output y ( t ) approaches s ( t ) [20]. FIGURE 2Open in figure viewerPowerPoint The calculation process of independent component analysis The Fast-ICA algorithm is a fast optimization iterative algorithm. The difference from other common neural network algorithms is that this algorithm uses a batch processing method, and a large amount of sample data participates in the operation in each iteration. This algorithm combines fixed-point iteration and non-Gaussian maximization algorithm, and has the advantages of good stability, fast convergence speed, and small amount of calculation. It uses the largest negative entropy as a search direction, and can extract independent sources sequentially, which fully embodies the traditional linear transformation idea of Projection Pursuit [21]. Because the Fast-ICA algorithm uses the largest negative entropy as a search direction, first discuss the negative entropy decision criterion. From the information theory, we know that among all random variables with equal variances, the Gaussian variable has the largest entropy, so we can use entropy to measure non-Gaussian variables. The commonly used modified form of entropy is negative entropy. According to the central limit theorem, if the random variable X consists of the sum of many independent random variables S i ( i = 1 , 2 , 3 , … , N ) , as long as S i has a finite mean and variance, no matter what distribution it is, the random variable X is closer to the Gaussian distribution than S i . In other words, S i is more non-Gaussian than X. Therefore, in the separation process, the mutual independence between the separation results can be expressed by the non-Gaussian measurement of the separation results. When the non-Gaussian measurement reaches the maximum, it indicates that the separation of the independent components has been completed. The negative entropy expression of a random variable is defined as N ( x ) = H ( x g ) − H ( x ) (12)where x g is a Gaussian random variable with the same variance as x, and H ( • ) is the differential entropy of the random variable: H ( x ) = − ∫ p x ( ξ ) lg p x ( ξ ) d ξ (13) According to information theory, among random variables with the same variance, Gaussian random variables have the largest differential entropy. When x has a Gaussian distribution, N ( x ) = 0 . The stronger the non-Gaussian of x, the smaller its differential entropy and the larger the value of N ( x ) , so N ( x ) can be used as a measure of the non-Gaussian of the random variable x. According to formula (13), the probability density distribution function x needs to be known to calculate the differential entropy, which is obviously impractical, so the following approximate formula is used N ( x ) = E g x − E g x g 2 (14)where E [ • ] is the mean value operation; g ( • ) is a non-linear function, usually g ( x ) = x exp ( − x 2 2 ) can be used. The Fast-ICA learning rule is to find a direction so that N ( W T x ) has the greatest non-Gaussian. According to formula (14), the approximate value of negative entropy can be calculated to measure non-Gaussianity. The variance of W T x is constrained to 1. For whitened data, this is equivalent to the norm of the constraint W being 1. The derivation of the Fast-ICA algorithm is as follows: First, the maximum approximation of the negative entropy of W T x can be obtained by optimizing E [ G ( W T x ) ] . According to the Kuhn-Tucker condition, under the constraint of E { ( W T x ) 2 } = ∥ W ∥ 2 = 1 , the optimal value of E [ G ( W T x ) ] can be obtained at the point that satisfies the following formula. Use Newton's method to calculate and simplify, the iterative formula can be obtained W k + 1 = E x g W k T x − E g ′ W k T x W k (15) Normalized W k + 1 = W k + 1 / W k + 1 (16)where g ( • ) is the first derivative of G ( • ) ; g ′ ( • ) is the second derivative of G ( • ) . The above iterative process is second order convergent. In order to improve the efficiency of the method for solving non-linear equations, the convergence order can be increased. The classic Newton iterative method is modified, and a two-step method of Newton's third order convergence with simple operation and stable convergence is proposed, which further simplifies the update iteration method of the improved Fast-ICA algorithm: W k + 1 ∗ = E x g W k T x − E g ′ W k T x W k (17) W k + 1 = E x g W k T x + E x g W k ∗ T x − E g ′ W k T x W k − β W k + 1 ∗ (18)where β = E { W k T x g ( W k T x ) } . Compared with the implementation steps of the traditional Fast-ICA algorithm, the improved Fast-ICA algorithm only differs in the update iteration steps, and the other steps are basically the same. The specific algorithm steps are as follows: Centralize the observation data x so that its mean value is 0; Whiten the data, x → z ; Choose the number m of components to be estimated, and set the number of iterations p ← 1 ; Choose an initial weight vector (random) W p ; Let { W k + 1 ∗ = E { x g ( W k T x ) } − E { g ′ ( W k T x ) } W k W k + 1 = E { x g ( W k T x ) } + E { x g ( W k + 1 ∗ T x ) } − E { g ′ ( W k T x ) } W k − E { W k T x g ( W k T x ) } W k + 1 ∗ , the selection of the non-linear function g, see the previous article; W p = W p − ∑ j = 1 p − 1 ( W p T W j ) W j .; Order W p = W p / ∥ W p ∥ ; If W p does not converge, go back to Step 5; Let p = p + 1 , if p ≤ m , go back to Step 4. 4 TARGET PRECESSION PARAMETER ESTIMATION 4.1 Three-dimensional cone rotation vector extraction The angle between the radar line of sight and the cone axis ε = cos − 1 ( ( 0 , 0 , 1 ) • n L O S ) (19) The angular frequencies of the micro-Doppler curves observed by radars with different observation angles are the same, but have different amplitudes. When the amplitude M c i of the micro-Doppler curve of the scattering point of the cone is extracted (the subscript i represents the serial number of the radar), the following formula can be obtained M c i = r sin ε i (20)where r is the radius of the circle of rotation at the apex of the cone. 4.2 Precession period, spin period extraction According to formula (8), it can be known that the radar echo mainly includes signals with four frequency components, namely ω c + ω s , ω c − ω s , ω c and ω s . The component with angular frequency ω c + ω s makes the amplitude of the echo oscillate at angular frequency ω c + ω s , and the other three angular frequency components together form the envelope of the echo curve. According to the literature [19], the specific expressions of A i and φ i can be known. In a precession period T pr , the number of maximum points of the micro-Doppler curve is N p = T pr ( Ω c + Ω s ) / ( 2 π ) , so the spin period can be obtained as T s = T pr T c N p T c − T pr (21)where T c = 2 π / Ω c . The precession period T pr can be obtained by the correlation method of the curve, and N p can be obtained by the extreme value search method. After extracting the cone spin vector, T c is a known quantity, so the spin period T s of the target can be obtained according to Equation (21) 4.3 Feature extraction of precession angle and bottom radius The micro-Doppler effect of the sliding scattering point should be modulated by the target spin and cone spin. In addition to extracting the precession period and the spin period, characteristic parameters such as the precession angle and the radius of the bottom surface can also be extracted from the radar echo. According to the literature [18], for the ith radar, there is a correlation between the maximum and minimum difference M p i of the micro-Doppler curve of the sliding scattering points B and C and the target precession parameters and structural parameters. Since the precession angle β is usually small, according to the related derivation in the literature [18], the following formula can be obtained: M p i = 2 r 1 cos β sin ε i + 2 r 1 sin β cos ε i (22)where r 1 is the radius of the bottom surface of the cone target. 4.4 Target length estimation The target length is equal to r 0 + r / r sin β sin β . It is assumed that the maximum length of the micro-Doppler curve of all the scattering points of the target during the precession period extracted on the radar line-of-sight range image is L max . According to the related derivation in the literature [18], the approximate estimated value of the target length is A O 1 ¯ ≈ ( ( L max − r sin ε i ) / cos ε i − r 1 sin β ) / cos β (23) Calculate the estimation error produced by this formula below, and record l x = 1 sin β r 0 cos β + r 1 sin β + cos ε i r 0 2 + r 1 2 (24)when l x < 0 , the error e l 1 = | r 1 − r 0 tan β | | tan ε i | (25)when 0 < l x < 2 r 1 , the error e l 2 = 1 − cos ( α − β ) cos ε i cos β cos ε i r 0 2 + r 1 2 (26)when l x > 2 r 1 , the error e l 3 = − 2 r 1 tan β + ( r 1 + r 0 tan β ) tan ε i (27) e l 1 is proportional to | tan ε | , and e l 2 is inversely proportional to | cos ε | . Since the precession angle is relatively small, − 2 r 1 tan β is usually relatively small, so e l 3 and | tan ε | are also approximately proportional. Therefore, when estimating the target length, the radar echo data with the largest | cos ε | should be selected for processing to minimize the target length estimation error. 5 SIMULATION The distance from the centre of mass of the cone to the vertex is l 1 = 2 m , the distance from the centre of the bottom surface is l 2 = 1 m , and the radius of the bottom surface is r = 1 m . The cone spin frequency of the target is ω c = 4 π , and the spin frequency is ω s = 5 π . The precession angle is θ = π / π 18 18 , and the initial Euler angle is ( ϕ e , θ e , φ e ) = ( π , π / π 18 18 , π ) . In the distributed networked radar, it is assumed that there are three radars. The radar bandwidth is 3 G H z , and the pulse width is T p = 40 μ s . The transmitted signals are all linear frequency modulation signals. The signal carrier frequency is f = 10 GHz , and the pulse repetition frequency is P R F = 1000 Hz . Radar 1 is located at the reference coordinate system ( 100 , 300 , − 500 ) k m , and the line of sight of the radar is n L O S 1 = ( − 1 / 1 35 35 , − 3 / − 3 35 35 , 5 / 5 35 35 ) . Radar 2 is located at the reference coordinate system ( − 200 , − 200 , − 500 ) k m , and the line of sight of the radar is n L O S 2 = ( 2 / 2 33 33 , 2 / 2 33 33 , 5 / 5 33 33 ) . Radar 3 is located at the reference coordinate system ( − 100 , − 400 , − 500 ) k m , and the line of sight of the radar is n L O S 3 = ( 1 / 1 42 42 , 4 / 4 42 42 , 5 / 5 42 42 ) . Suppose the signal-to-noise ratio is S N R = 0 dB . 5.1 Signal separation Before the signal is separated, the flat component is compensated by utilizing a corner detection algorithm [22]. Therefore, the factors of the flat component are not considered in the experiment. Take Radar 1 as an example to carry out a simulation experiment of signal separation. According to the distance formula and related parameters obtained above, the cone target movement is simulated, and the radar echo with a duration of 2 s is received together, and the target micro-Doppler curve as shown in Figure 3 is obtained. FIGURE 3Open in figure viewerPowerPoint Target micro-Doppler curve According to the principle described in Section 3.2, after the target echo is obtained, the type of scattering point can be determined. Using formula (10), we can analyze how fast the n-order micro-distance changes at different scattering points. Figure 4 shows how fast the micro-Doppler curves of different scattering centres change over time at n = 3. It can be seen from Figure 4 that at n = 3, the third-order micro-motion distance between the scattering centre of the cone top and the scattering centre of the bottom edge sliding is very different. Therefore, according to the principle described above, the type of the scatter centre of the cone target can be more accurately identified through the image. The model assumed here is relatively simple. If the model is relatively complex, such as with a tail wing, it can also be determined using the above-mentioned principles and methods. FIGURE 4Open in figure viewerPowerPoint Order micro-movement distance of different scattering centres Figure 5 simulates the separation of the obtained cone target micro-Doppler signal according to the principle of the Fast-ICA algorithm FIGURE 5Open in figure viewerPowerPoint Signal after separation First, set any matrix parameter to process the observed signal, and then centralize and whiten the processed signal. Finally, realize the separation of the signal. The separated signal is shown as in Figure 5 According to Figure 5 it can be seen that the micro-Doppler curves of the scattering points