Micro-Motion Parameter Extraction for Ballistic Missile with Wideband Radar Using Improved Ensemble EMD Method

: Micro-motion parameters extraction is crucial in recognizing ballistic missiles with a wideband radar. It is known that the phase-derived range (PDR) method can provide a sub-wavelength level accuracy. However, it is sensitive and unstable when the signal-to-noise ratio (SNR) is low. In this paper, an improved PDR method is proposed to reduce the impacts of low SNRs. First, the high range resolution proﬁle (HRRP) is divided into a series of segments so that each segment contains a single scattering point. Then, the peak values of each segment are viewed as non-stationary signals, which are further decomposed into a series of intrinsic mode functions (IMFs) with different energy, using the ensemble empirical mode decomposition with the complementary adaptive noise (EEMDCAN) method. In the EEMDCAN decomposition, positive and negative adaptive noise pairs are added to each IMF layer to effectively eliminate the mode-mixing phenomenon that exists in the original empirical mode decomposition (EMD) method. An energy threshold is designed to select proper IMFs to reconstruct the envelop for high estimation accuracy and low noise effects. Finally, the least-square algorithm is used to do the ambiguous phases unwrapping to obtain the micro-curve, which can be further used to estimate the micro-motion parameters of the warhead. Simulation results show that the proposed method performs well with SNR at − 5 dB with an accuracy level of sub-wavelength.


Introduction
In the missile defense system, it is crucial to distinguish the true warheads from decoys in the mid-course phase, where warheads and decoys exhibit different micromotions. Specifically, the true warheads have precession and nutation movements, while decoys wobble after being released from the true warheads. When observing the midcourse ballistic missiles using radar, the precession and nutation movements will modulate the echoes. Consequently, it is possible to extract micro-motion parameters from the radar echoes to recognize the true warheads [1][2][3][4][5][6].
In the past few decades, many methods have been developed to exploit the micromotion features to recognize the mid-course ballistic missiles [7][8][9][10][11][12]. It is difficult for the narrow-band radar to discriminate the different scattering points lying on the true warhead due to its low range resolution. Thus, the warhead is usually considered as a whole, and many methods based on the whole radar cross section (RCS) are developed [13][14][15]. RCS is affected by various factors such as target shape, target material, illuminating angle, etc. As a result, the estimated micro-motion parameters based on the RCS sequence may not be fine, accurate, and reliable enough [16]. By contrast, the wideband radar is much more appealing. The high range resolution profile (HRRP) benefiting from the large available bandwidth makes it possible to separate the scattering points from each other. Therefore, many more advanced methods for the wideband radar have been proposed to extract finer and more accurate micro-motion features [17][18][19][20][21].
It is reported that the phase-derived range (PDR) method based on HRRP sequences can achieve the half-wavelength accuracy, and thus has attracted increasing attentions [22,23]. It measures the target range using the phase information based on R = (λ/4π)∆φ, where λ and ∆φ represent the wavelength and the unambiguous phase, respectively. However, there are two challenges for the the PDR method in the wideband radar. One is how to extract the ambiguous phase, corresponding to the envelop across the HRRP sequence. As the scattering point will occupy several range bins, it is particularly significant to determine the correct range bins to fit the envelop. The other is how to unwrap the ambiguous phase, as ∆φ is unambiguous only with ∆φ < 2π. In other words, the maximum unambiguous range is λ/2, which is impractical. Therefore, the phase unwrapping is indispensable.
Note that several improved PDR methods have been applied for the micro-motion parameters extraction. In [7], a micro-Doppler phase delta method was proposed to tackle the phase ambiguity. It calculates the phase difference between two adjacent pulses, resolves the number of phase ambiguities, and further uses the PDR method to properly measure the distance. However, as mentioned in [7], it requires several related measurement errors to be small enough to correctly resolve the phase ambiguity. In other words, it requires high signal-to-noise ratio (SNR). In [17], the authors proposed a new PDR method to extract Doppler phases from HRRP sequences through a matched filter. This algorithm was developed in [7], where they applied the same method to the direct sampling linear frequency modulated signal to avoid the influences of the recording error. However, the high SNR requirement and the limitation to stationary targets remained unresolved. In [10], the authors presented a high-precision PDR method with no constraints on SNRs for phase ambiguity. The ambiguous phase is estimated by applying the conjugate multiplication between two adjacent pulses, and relaxes the requirement of SNR. However, the variation range between adjacent frames is required to be less than one-quarter of the wavelength, which limits its application.
As is well known, EMD is a powerful tool for the non-stationary signal analysis [24,25]. Its main advantage is that it can decompose the non-stationary signal into a series of intrinsic mode functions (IMFs). The non-stationary signal can be reconstructed by the filtered IMFs to envelop association for micro-curve extraction [26]. However, the original EMD algorithm heavily relies on the extreme points of the non-stationary signal. Besides, EMD is greatly affected by the high-frequency interference in the signal.
In this paper, we develop an improved PDR method and attempt to deal with the micro-motion parameters extraction in the wideband radar. To overcome the low-SNR issue and high-frequency interference, instead of directly applying EMD, we adopt an ensemble EMD with complementary adaptive noise (EEMDCAN) [27] to extract the scattering point envelop across the HRRPs. First, we divide the HRRP into a series of segments such that each of them contains a scattering point. We take the peak value of each segment along the HRRP as the coarse estimation of the corresponding scattering point, and then regard all coarse estimations of the scattering point across the HRRPs as a non-stationary signal. Then, the EEMDCAN method is applied to each non-stationary signal to obtain a series of IMFs for each scattering point. The decomposition terminates when the conditions of modal amplitude and estimation function proposed by Rilling G. are satisfied [28]. Only the IMFs with a zero-crossing rate (ZCR) larger than a threshold will be employed to reconstruct the accurate envelop, while the remaining IMFs will be considered as interference and discarded [28]. Based on the accurate envelop, the ambiguous phases can be directly obtained, and the least-squares algorithm is employed to unwrap the ambiguous phases. Finally, parameters estimation is performed. For the target with precession, the microcurve of scattering points can be regarded as sinusoidal curves, of which the amplitudes are affected by the precession angle, cone height, and radius. Therefore, the precession angle, precession frequency, cone height, and radius can be obtained through the estimated micro-curve.
The most important contribution of this paper is that an improved PDR method is proposed to overcome the difficulties in recognizing the true warheads with low SNRs. Specifically, the proposed method can adaptively estimate the envelop according to the characteristics of the signal. Even when the SNR is lower than −5 dB, this method can still estimate the envelop with high accuracy at a sub-wavelength level.
The remaining of this paper is organized as follows. Section 2 introduces the geometric model and related backgrounds. Section 3 presents the parameter estimation method, which contains scattering point separation, micro-curve extraction, and parameter extraction methods. In Section 4, simulation results are provided to verify the proposed method. Section 5 concludes the paper.

Geometry and Signal Model
The cone model of a mid-course warhead observed by a wideband radar is illustrated in Figure 1, where Figure 1a presents the radar coordinate system (Q − UVW), the warhead local coordinate system (O − xyz), and the reference coordinate system (O − XYZ), while Figure 1b presents the details of a warhead. During the observation, the radar keeps stationary.

Warhead Motion
In general, the warhead is modeled as a rigid body, where the distance between any two particles does not vary during any motion [1]. For cone targets, precession is a common form of micro-motion, i.e., when the target is moving, it is accompanied by periodic cone rotation around the cone rotation axis [8]. As Figure 1 shows, the warhead is modeled as a cone, and in a mid-course warhead, there usually exist two motions [8,29]: one is the warhead spinning around the symmetric axis − → Oz with the angular velocity ω s , and the other is the warhead rotating around the precession direction − → OC with angular velocity ω c . The angle ϕ between − → OC and the cone symmetric axis − → Oz is referred as to the precession angle. In Figure 1b, H, σ, r, and d stand for the height, half opening angle, the base radius of the cone, and the distance from the top to O (namely, − − → OP 1 ), respectively. According to the scattering theory, there are usually three strong scattering points in the HRRP sequences of the cone target, namely P 1 , P 2 , and P 3 , where P 1 is fixed, while P 2 and P 3 change their positions as the radar line of sight (LOS) changes. In the radar coordinate system, let α and β denote the initial azimuth and elevation angle of the warhead, respectively; − → R 0 = (U 0 , V 0 , W o ) T denote the initial position vector of the origin of the reference coordinate system O; and − → R 0 = R 0 , where · represents the Euclidean norm, then the unit vector of radar LOS is Rotation will change the target posture, and consequently change the position of the scattering point relative to the radar. In the target coordinate system, assume the scattering point P is with a position vector − → r 0 = (X 0 , Y 0 , Z o ) T at the initial time t = 0. Let (φ, θ, ψ) denote the initial Euler angles, where φ, θ, and ψ rotate about z-axis, x-axis, and z-axis, respectively, and the corresponding initial rotation matrix is Then, the new location of P in the reference coordinate system can be represented by According to the Rodrigues formula [8], the rotation matrix at time t can be written as where I is the unit matrix, Ω = − → ω c , and K =ŵ c is a skew symmetric matrix of the unit direction vector: Substituting (4) into (3), the rotation matrix becomes cos Ωt − cos ϕ sin Ωt sin ϕ sin Ωt cos ϕ sin Ωt 1 − cos 2 ϕ + cos 2 ϕ cos Ωt sin ϕ cos ϕ − sin ϕ cos ϕ cos Ωt − sin ϕ sin Ωt sin ϕ cos ϕ − sin ϕ cos ϕ cos Ωt The distance from the radar to the point P can be derived as where − → V is the velocity vector, i.e., and R 0 is the initial range value. Ballistic targets are usually far away from the radar, and the far-field condition can be satisfied. Thus, the instantaneous radial distance from the scattering point to the radar can be calculated by the projection of the scattering point to the radar's LOS, namely, where · denotes the inner product of two vectors, and γ denotes the angle between target symmetry axis and radar line of sight at the time t, cos γ = − − → LOS, − → Z , i = 1, 2, 3. We assume that the radar bandwidth is large enough so that it can discriminate each scattering point of the cone target, and the cone target will span several range bins. Therefore, the micro-motions of the scattering points could be distinguished from the bulk of the target in the radial range dimension. In this paper, we ignore the influence of bulk motion, and assume R 0 to be 0. According to the works in [30], the distance between scattering points and radar can be expressed as where ρ = sin ϕ sin(ϕ + γ) cos w c t, and R h represents the radial distance change caused by the position of the local coordinate center O

Radar Echoes
The linear frequency modulation (LFM) signals are widely used in wideband radar systems [10] and can be expressed as , t denotes the fast time, T p denotes the pulse duration, f c denotes the carrier frequency, and k denotes the slope frequency modulation. Then, the signal returned from the cone can be modeled as the summation of the returns from all scattering points [31]: where t m = mT r is the slow time, A i is the P i reflectivity, λ = c/ f c represents the carrier wavelength, c is the wave propagation velocity, n r (t m ) denotes the noise term of the received echo signal, and N denotes the number of scattering points. After pulse compression, the obtained HRRPs can be represented as where A i denotes the amplitude after pulse compression, and A i (t m ) represents the amplitude containing the sinc function.

Scattering Point Separation
The attitude of the warhead is relatively stable in the mid-course phase, and the envelops of the different scattering points do not cross each other. Thus, HRRP sequences can be divided into corresponding segments according to the number of scattering points. The HRRP sequences of each scattering point can be presented as Then, the signal of the m-th pulse of the i-th scattering point in the slow-time dimension can be expressed as

Micro-Curve Extraction
HRRP sequences can be used to distinguish different scattering points, but it is challenging to obtain the peak value for the envelop from numerous range-bins. As shown in (16), the energy distribution of the scattering points is close to sinc functions, and the peaks are sensitive to noise. Thus, it is required to determine the accurate scattering point position in each HRRP. The EEMDCAN has good time-frequency characteristics, so it is introduced to deal with the nonlinear and non-stationary signals. To reduce the interference of the residual noise remaining in the IMF component, positive and negative noises are added in the decomposition process of each layer of the EEMDCAN approach.
As shown in Figure 2, the proposed procedure is described in detail as follows.
Step 1: For each scattering point segment, select the peak as the estimated scattering point. As a result, we can obtain a rough envelop, namely, s(m), s(m) = s 1 , s 2 , · · · , s m .
Step 2: Add the pairs of positive and negative white noises (−1) q a 0 n j (m) to the sequence s(m), and the mixed signal s(m) + (−1) q a 0 n j (m) is achieved. Perform EMD decomposition on the mixed signals, and stop the EMD decomposition immediately after filtering out the first IMFs, Step 3: Subtracting im f 1 (m) from s(m) and the achieved residual signal l 1 (m) is Define a k and E k () as the standard deviation of the kth residual signal and IMF of the signal by the EMD, respectively. Add the pairs of positive and negative adaptive noises (−1) q a 1 E 1 n j (m) to the residual signal l 1 (m) and then obtain the new signal l 1 (m) + (−1) q a 1 E 1 n j (m) . Decomposing the new signal by EMD will achieve the second ensemble im f j 2 (m) and the corresponding residual signal l 2 (m), we have Step 4: The pairs of positive and negative adaptive noises (−1) q a k E k n j (m) are added to the residual signal to achieve the new signal l k (m) + (−1) q a k E k n j (m) , where a k is usually chosen as 0.1 ∼ 0.2 times of the standard deviation of the new residual signal l k (m) [28]. Then, the (k + 1)-th ensemble im f k+1 (m) is Step 5: Repeat step 4 until all modal components are extracted. With all the ensemble mean im f k (m), the final residual component L(m) is obtained as follows: where G is the total number of the ensemble mean of IMFs. We employ the modal amplitude and estimation function proposed by Rilling G. as the condition to terminate the decomposition [28]. After EEMDCAN decomposition, s(m) is decomposed into a series of modal components and residual components, merely IMF with ZCR greater than 10 is used to reconstruct the envelop.ŝ (m) = ŝ 1 ,ŝ 2 , · · · ,ŝ m = where D represents the maximum order of IMF that meets the threshold, and the corresponding ambiguous phase iŝ Re(ŝ 2 ) , · · · , arctan Im(ŝ m ) Re(ŝ m ) .
Next, the least-squares method [32] is employed to find the real phase Φ(m). Then, the micro-curve can be estimated by the PDR as

Precession Frequency Estimation
According to the works in [1], for the precession target, the micro-curve of predicted scattering points can be regarded as sinusoidal curves. The Fourier transform is performed to obtain the spectrum of the micro-curve: and the precession frequency is the frequency corresponding to the peak point [33].

Precession Angle Estimation
As the amplitude of the micro-curve is coupled with the precession angle, the height, and radius of the cone, the micro-curve is normalized to eliminate the influence of the amplitude factor. After excluding the influence of the amplitude factor, the spectrum amplitude of the micro-curve is only affected by the precession angle. If the estimated precession frequency is regarded as a known parameter, the precession angle estimation problem can be transformed into an optimization problem aŝ whereφ denotes the precession angle, which is usually less than 15 • to ensure the stability of the target's attitude.

Half-Cone Angle and Length of Busbar Estimation
As shown in Figure 1, the length of bus-bar S, half-cone angle, and the distance between the top and the bottom of the cone can be expressed as S = √ r 2 + H 2 , σ = arctan(r/H), and ξ = R 1 − R 3 = H cos γ(t) − r sin γ(t), and then where (36) can be used to search the half cone angle, and R is the distance difference between two scattering points in the estimation micro-curve.

Cone Height and Base Radius Estimation
The cone height and the base radius can be obtained according to the properties of the geometric parameters, H = S cos(σ), r = S sin(σ).

Simulations and Discussion
This section shows the numeric simulations carried out to evaluate the proposed method. The radar parameters and warhead parameters are shown in Table 1.  Figure 3 shows the HRRP sequences with SNR equal to 15 dB and it can be seen that the wideband radar can distinguish the scattering points. In addition, the cone target occupies some range bins due to its volume. Note that it is also possible to divide HRRP sequences into two parts: one contains the cone-top scattering point, and the other contains the cone-bottom scattering point.

Micro-Curve Extraction
Micro-curve extraction is divided into two steps: one is envelop extraction, and the other is phase unwrapping. Envelop extraction includes peak signal decomposition and signal reconstruction. It can be observed from Figure 4 that the extracted peaks of the conetop and cone-bottom scattering point exhibit nonlinear and non-stationary characteristics. We adopt the original EMD method and the EEMDCAN method to decompose the nonstationary peak envelops, and obtain the results shown in Figure 5. We can see that the first-order IMF contains the highest frequency energy, and there is no mode mixing in the decomposition process. However, as the red arrows indicate, the obtained IMF components by original EMD are obviously mixed with each other. With the obtained IMFs by EEMDCAN, we reconstruct the peak envelops, which are presented in Figure 4a,b as red dashed lines.  Based on the reconstructed envelops, we calculate the phase distributions in Figure 6a,d. When the range changing is greater than one range bin (i.e., c/2B) in the radial range dimension, namely, across-range unit migration, phase increases by −π or +π and consequently phase hopping occurs. Thus, phase compensation is required before the phase unwrapping to guarantee the accuracy. We calculate the phase difference between two adjacent HRRPs, and take the points where the difference is greater than 0.5 as phase hopping, as depicted in Figure 6b,e. Then, additional compensation of π or −π is added and a smooth sine-like curve will be achieved, as shown in Figure 6c,f. Besides, Figure 7 depicts that the micro-curve extraction using the EMD and EEMDCAN methods with different SNR values. It can be observed that the micro-curves extraction using EEMD-CAN method is much closer to the theoretical curve, compared to the one using the EMD method. Especially, when the SNRs is low, say −10 dB, the EEMDCAN method obviously outperforms the EMD method.

Parameters Estimation
With the extracted and unwrapped phase distribution curves via the EEMDCAN and EMD methods, we compare the estimated warhead cone parameters in Table 2, where 100 times Monte Carlo simulations are carried out for each SNR value. The results also validate that EEMDCAN outperforms EMD with higher accuracy.
We compare the proposed method with the modified Kalman filter (MKF) method proposed in [4] with the same parameters listed in Table 1. The RMSE with different SNRs are given in Figure 8, where the RMSE is obtained by where r(m) andŝ(m) denote the theoretical and estimated values of micro-curve, respectively.
In the simulation, the MKF method is used to fit envelops and calculate the ambiguous phases, and the least-square method is adopted to phase unwrapping. Finally, the PDR method is used to obtain the micro-curves. Both the MKF and the proposed methods achieved the sub-wavelength level accuracy.
It is shown in Figure 8 that the proposed method can achieve lower RMSE values, especially when SNR is low. Note that in [4], when using the MKF method to fit the envelop, it is assumed the Kalman filter is a constant velocity model, which may cause fitting distortion.

Conclusions
In this paper, we proposed an improved PDR method to conquer the difficulties in recognizing the true warheads using wideband radar with low SNRs. Specifically, the proposed method can adaptively estimate the envelop according to the characteristics of the signal. Even when the SNR is lower than −5 dB, this method can still estimate the envelop with high accuracy at a sub-wavelength level. Numerical simulations validate the anti-noise performance and robustness of the proposed method. The directions of future research could be on the scenario where scattering points are not distinguishable.