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

Abstract

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

1. 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 micro-motions. Specifically, the true warheads have precession and nutation movements, while decoys wobble after being released from the true warheads. When observing the mid-course 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 micro-motion 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=(\lambda /4\pi )\Delta \varphi$, where $\lambda$ and $\Delta \varphi$ 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 $\Delta \varphi$ is unambiguous only with $\Delta \varphi <2\pi$. In other words, the maximum unambiguous range is $\lambda /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 micro-curve 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.

2. 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^{{}^{\prime }}-XYZ$), while Figure 1b presents the details of a warhead. During the observation, the radar keeps stationary.

2.1. 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 $\overrightarrow{Oz}$ with the angular velocity ${\omega }_s$, and the other is the warhead rotating around the precession direction $\overrightarrow{OC}$ with angular velocity ${\omega }_c$. The angle $\phi$ between $\overrightarrow{OC}$ and the cone symmetric axis $\overrightarrow{Oz}$ is referred as to the precession angle. In Figure 1b, H, $\sigma$, 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, $\parallel \overrightarrow{O{P}_1}\parallel$), 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 $\alpha$ and $\beta$ denote the initial azimuth and elevation angle of the warhead, respectively; $\overrightarrow{R_0}={(U_0,V_0,W_o)}^T$ denote the initial position vector of the origin of the reference coordinate system O; and $\parallel \overrightarrow{R_0}\parallel =R_0$, where $\parallel ·\parallel$ represents the Euclidean norm, then the unit vector of radar $LOS$ is

$$\overrightarrow{n}=\frac{\overrightarrow{R_0}}{R_0}=\left(cos\alpha cos\beta ,sin\alpha cos\beta ,sin\beta \right).$$

(1)

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 $\overrightarrow{r_0}={\left(X_0,Y_0,Z_o\right)}^T$ at the initial time $t=0$. Let $\left(\varphi ,\theta ,\psi \right)$ denote the initial Euler angles, where $\varphi$, $\theta$, and $\psi$ rotate about z-axis, x-axis, and z-axis, respectively, and the corresponding initial rotation matrix is

$$R_{Init}=\left[\begin{array}{ccc}cos\varphi & -sin\varphi & 0\\ sin\varphi & cos\varphi & 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& cos\theta & -sin\theta \\ 0& sin\theta & cos\theta \end{array}\right]\left[\begin{array}{ccc}cos\psi & -sin\psi & 0\\ sin\psi & cos\psi & 0\\ 0& 0& 1\end{array}\right].$$

(2)

Then, the new location of P in the reference coordinate system can be represented by $R_{Init}·\overrightarrow{r_0}$. Supposing that the angular vector is $\overrightarrow{w_c}={(w_x,w_y,w_z)}^T$, then the unit vector of the rotation is ${\overrightarrow{w_c}}^{{}^{\prime }}={({w_x}^{{}^{\prime }},{w_y}^{{}^{\prime }},{w_z}^{{}^{\prime }})}^T=R_{Init}·\frac{\overrightarrow{w_c}}{\parallel \overrightarrow{w_c}\parallel }$. According to the Rodrigues formula [8], the rotation matrix at time t can be written as

$$R\left(t\right)=I+Ksin\mathsf{\Omega }t+K^2\left(1-cos\mathsf{\Omega }t\right),$$

(3)

where I is the unit matrix, $\mathsf{\Omega }=\parallel \overrightarrow{{\omega }_c}\parallel$, and $K={\widehat{w_c}}^{{}^{\prime }}$ is a skew symmetric matrix of the unit direction vector:

$$K=\left[\begin{array}{ccc}0& -{w_{cz}}^{{}^{\prime }}& {\mathrm{w}}_{\mathrm{cy}}^{{}^{\prime }}\\ {w_{cz}}^{{}^{\prime }}& 0& {-w_{cx}}^{{}^{\prime }}\\ {-w_{cy}}^{{}^{\prime }}& {w_{cx}}^{{}^{\prime }}& 0\end{array}\right]=\left[\begin{array}{ccc}0& -cos\phi & sin\phi \\ cos\phi & 0& 0\\ -sin\phi & 0& 0\end{array}\right].$$

(4)

Substituting (4) into (3), the rotation matrix becomes

$$R\left(t\right)=\left[\begin{array}{ccc}cos\mathsf{\Omega }t& -cos\phi sin\mathsf{\Omega }t& sin\phi sin\mathsf{\Omega }t\\ cos\phi sin\mathsf{\Omega }t& 1-{cos}^2\phi +{cos}^2\phi cos\mathsf{\Omega }t& sin\phi cos\phi -sin\phi cos\phi cos\mathsf{\Omega }t\\ -sin\phi sin\mathsf{\Omega }t& sin\phi cos\phi -sin\phi cos\phi cos\mathsf{\Omega }t& 1-{sin}^2\phi +{sin}^2\phi cos\mathsf{\Omega }t\end{array}\right].$$

(5)

The distance from the radar to the point $P^{{}^{\prime }}$ can be derived as

$$r\left(t\right)=\parallel \overrightarrow{Q{P}^{{}^{\prime }}}\parallel =\parallel \overrightarrow{QO}+\overrightarrow{O{O}^{{}^{\prime }}}+\overrightarrow{O^{{}^{\prime }}{P}^{{}^{\prime }}}\parallel =\parallel \overrightarrow{R_0}+\overrightarrow{V}t+R_t\overrightarrow{r_0}\parallel ,$$

(6)

where $\overrightarrow{V}$ is the velocity vector, i.e., $\overrightarrow{O{O}^{{}^{\prime }}}=\overrightarrow{V}t$, 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,

$$r_i\left(t\right)=〈\overrightarrow{LOS},\overrightarrow{P_i}〉+R_0=Y_{i}sin\gamma +Z_{i}cos\gamma +R_0,$$

(7)

where $〈·〉$ denotes the inner product of two vectors, and $\gamma$ denotes the angle between target symmetry axis and radar line of sight at the time t, $cos\gamma =〈\overrightarrow{LOS},\overrightarrow{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

$$r_1\left(t\right)=-dcos\phi cos\left(\phi +\alpha \right)-Hsin\phi sin\left(\phi +\alpha \right)cos\mathsf{\Omega }t+R_h,$$

(8)

$$r_2\left(t\right)=r\sqrt{1-{cos}^2\phi {cos}^2\left(\phi +\gamma \right)}+\left(\begin{array}{c}H-d-r\frac{cos\phi cos\left(\phi +\gamma \right)}{\sqrt{1-{cos}^2\phi {cos}^2\left(\phi +\gamma \right)}}\end{array}\right)·\rho +R_h,$$

(9)

$$r_3\left(t\right)=-r\sqrt{1-{cos}^2\phi {cos}^2\left(\phi +\gamma \right)}+\left(\begin{array}{c}H-d+r\frac{cos\phi cos\left(\phi +\gamma \right)}{\sqrt{1-{cos}^2\phi {cos}^2\left(\phi +\gamma \right)}}\end{array}\right)·\rho +R_h,$$

(10)

where $\rho =sin\phi sin\left(\phi +\gamma \right)cos{w}_{c}t$, and $R_h$ represents the radial distance change caused by the position of the local coordinate center O

$$R_h=\left(H-d\right)sin\phi sin\left(\phi +\alpha \right)cos\left(\mathsf{\Omega }t\right).$$

(11)

2.2. Radar Echoes

The linear frequency modulation (LFM) signals are widely used in wideband radar systems [10] and can be expressed as

$$s\left(t\right)=\mathrm{rect}\left(\frac{t}{T_p}\right)exp\left[\begin{array}{c}j2\pi \left(f_{c}t+\frac{1}{2}k{t}^2\right)\end{array}\right],$$

(12)

where $\mathrm{rect}\left(\zeta \right)=\left\{\begin{array}{c}1,\left|\begin{array}{c}\zeta \end{array}\right|\le 0.5\hfill \\ 0,\left|\begin{array}{c}\zeta \end{array}\right|>0.5\hfill \end{array}\right.$, 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]:

$$s\left(t,t_m\right)=\sum _{i=1}^N{A}_i\mathrm{rect}\left(\frac{t-2r_i\left(t_m\right)/c}{T_p}\right)exp\left\{-j\frac{4\pi r_i\left(t_m\right)}{\lambda }\right\}exp\left\{j\pi k{\left(t-\frac{2r_i\left(t_m\right)}{c}\right)}^2\right\}+n\left(t_m\right),$$

(13)

where $t_m=m{T}_r$ is the slow time, $A_i$ is the $P_i$ reflectivity, $\lambda =c/f_c$ represents the carrier wavelength, c is the wave propagation velocity, $n_r\left(t_m\right)$ 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

$$\begin{array}{cc}\hfill S\left(t,t_m\right)& =\sum _{i=1}^N{A}_i\sin \mathrm{c}\left[B\left(t-\frac{2r_i\left(t_m\right)}{c}\right)\right]exp\left\{-j\frac{4\pi r_i\left(t_m\right)}{\lambda }\right\}+n\left(t_m\right)\hfill \\ \hfill & =\sum _{i=1}^N{A}_i^{{}^{\prime }}\left(t_m\right)\exp \left\{-j\frac{4\pi r_i\left(t_m\right)}{\lambda }\right\}+n\left(t_m\right),\hfill \end{array}$$

(14)

where $A_i^{{}^{\prime }}$ denotes the amplitude after pulse compression, and $A_i^{{}^{\prime }}\left(t_m\right)$ represents the amplitude containing the sinc function.

3. Micro-Motion Parameter Estimation

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

$$\begin{array}{cc}\hfill S_i\left(t_m\right)& =A_i^{{}^{\prime }}\left(t_m\right)\exp \left\{-j\frac{4\pi r_i\left(t_m\right)}{\lambda }\right\}+n\left(t_m\right).\hfill \end{array}$$

(15)

Then, the signal of the m-th pulse of the i-th scattering point in the slow-time dimension can be expressed as

$$\begin{array}{cc}\hfill S_i\left(m\right)& =A_i^{{}^{\prime }}\left(m\right)\exp \left\{-j\frac{4\pi r_i\left(m\right)}{\lambda }\right\}+n\left(m\right).\hfill \end{array}$$

(16)

3.2. 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\left(m\right)$,

$$s\left(m\right)=\left[\begin{array}{c}{s}_1,s_2,\cdots ,s_m\end{array}\right].$$

(17)

Step 2: Add the pairs of positive and negative white noises ${\left(-1\right)}^q{a}_0{n}^j\left(m\right)$ to the sequence $s\left(m\right)$, and the mixed signal $s\left(m\right)+{\left(-1\right)}^q{a}_0{n}^j\left(m\right)$ is achieved. Perform EMD decomposition on the mixed signals, and stop the EMD decomposition immediately after filtering out the first IMFs,

$$s\left(m\right)+{(-1)}^q{a}_0{n}^j\left(m\right)=im{f}_1^j\left(m\right)+l_1^j\left(m\right),$$

(18)

where $q=1,2$, $j=1,2,\cdots ,M/2$, M represents the number of the added positive and negative white noises, and $n^j\left(m\right)$ denotes the jth white Gaussian noise. The ensemble mean $\overline{im{f}_1\left(m\right)}$ of all the $im{f}_1^j\left(m\right)$ can be calculated as

$$\overline{im{f}_1\left(m\right)}=\frac{1}{M}\sum _{j=1}^{M}im{f}_1^j\left(m\right)=s\left(m\right)-\frac{1}{M}\sum _{j=1}^M{l}_1^j\left(m\right).$$

(19)

Step 3: Subtracting $\overline{im{f}_1\left(m\right)}$ from $s\left(m\right)$ and the achieved residual signal $l_1\left(m\right)$ is

$$l_1\left(m\right)=s\left(m\right)-\overline{im{f}_1\left(m\right)}.$$

(20)

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 ${\left(-1\right)}^q{a}_1{E}_1\left(n^j\left(m\right)\right)$ to the residual signal $l_1\left(m\right)$ and then obtain the new signal $l_1\left(m\right)+{\left(-1\right)}^q{a}_1{E}_1\left(n^j\left(m\right)\right)$. Decomposing the new signal by EMD will achieve the second ensemble $im{f}_2^j\left(m\right)$ and the corresponding residual signal $l_2\left(m\right)$, we have

$$l_1\left(m\right)+{\left(-1\right)}^q{a}_1{E}_1\left(n^j\left(m\right)\right)=im{f}_2^j\left(m\right)+l_2^j\left(m\right),$$

(21)

$$l_2\left(m\right)=l_1\left(m\right)-\overline{im{f}_2\left(m\right)}.$$

(22)

Step 4: The pairs of positive and negative adaptive noises ${\left(-1\right)}^q{a}_k{E}_k\left(n^j\left(m\right)\right)$ are added to the residual signal to achieve the new signal $l_k\left(m\right)+{\left(-1\right)}^q{a}_k{E}_k\left(n^j\left(m\right)\right)$, where $a_k$ is usually chosen as $0.1\sim 0.2$ times of the standard deviation of the new residual signal $l_k\left(m\right)$ [28]. Then, the $\left(k+1\right)$-th ensemble $\overline{im{f}_{k+1}\left(m\right)}$ is

$$\overline{im{f}_{k+1}\left(m\right)}=\frac{1}{M}\sum _{j=1}^{M}im{f}_{k+1}^j\left(m\right)=l_k\left(m\right)-\frac{1}{M}\sum _{j=1}^M{l}_{k+1}^j\left(m\right).$$

(23)

Step 5: Repeat step 4 until all modal components are extracted. With all the ensemble mean $\overline{im{f}_k\left(m\right)}$, the final residual component $L\left(m\right)$ is obtained as follows:

$$L\left(m\right)=s\left(m\right)-\sum _{k=1}^G\overline{im{f}_k\left(m\right)},$$

(24)

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\left(m\right)$ is decomposed into a series of modal components and residual components, merely IMF with ZCR greater than 10 is used to reconstruct the envelop.

$$\widehat{s}\left(m\right)=\left[\begin{array}{c}{\widehat{s}}_1,{\widehat{s}}_2,\cdots ,{\widehat{s}}_m\end{array}\right]=\sum _{k=1}^D\overline{im{f}_k\left(m\right)},$$

(25)

where D represents the maximum order of IMF that meets the threshold, and the corresponding ambiguous phase is

$$\begin{array}{cc}\hfill \widehat{\Phi }\left(m\right)& =\left[\begin{array}{c}\Phi \left({\widehat{s}}_1\right),\Phi \left({\widehat{s}}_2\right),\cdots ,\Phi \left({\widehat{s}}_m\right)\end{array}\right]\hfill \\ \hfill & =\left[\begin{array}{c}arctan\frac{Im\left({\widehat{s}}_1\right)}{Re\left({\widehat{s}}_1\right)},arctan\frac{Im\left({\widehat{s}}_2\right)}{Re\left({\widehat{s}}_2\right)},\cdots ,arctan\frac{Im\left({\widehat{s}}_m\right)}{Re\left({\widehat{s}}_m\right)}\end{array}\right].\hfill \end{array}$$

(26)

Next, the least-squares method [32] is employed to find the real phase $\tilde{\Phi }\left(m\right)$. Then, the micro-curve can be estimated by the PDR as

$$R\left(m\right)=\frac{\lambda }{4\pi }\tilde{\Phi }\left(m\right).$$

(27)

3.3. Micro-Motion Parameter Estimation

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

$$F=\mathrm{FFT}\left\{\widehat{s}\left(m\right)\right\},$$

(28)

and the precession frequency is the frequency corresponding to the peak point [33].

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

$$\widehat{\phi }=\underset{\phi }{error}\left\{\begin{array}{c}\left|\begin{array}{c}F-\mathrm{FFT}\frac{\widehat{s}\left(m\right)}{\left|\begin{array}{c}\widehat{s}\left(m\right)\end{array}\right|}\end{array}\right|\end{array}\right\},$$

(29)

where $\widehat{\phi }$ denotes the precession angle, which is usually less than ${15}^{°}$ to ensure the stability of the target’s attitude.

3.3.3. 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=\sqrt{r^2+H^2}$, $\sigma =\arctan \left(r/H\right)$, and $\mathsf{\xi }=R_1-R_3=Hcos\gamma \left(t\right)-rsin\gamma \left(t\right)$, and then

$${\mathsf{\xi }}_{max}=Scos\left(\alpha -\phi +\sigma \right),$$

(30)

$${\mathsf{\xi }}_{min}=Scos\left(\alpha +\phi +\sigma \right),$$

(31)

$$S=\frac{{\mathsf{\xi }}_{max}-{\mathsf{\xi }}_{min}}{2cos\left(\alpha +\sigma \right)cos\phi },$$

(32)

$$\mathsf{\xi }\left(t\right)=\frac{{\mathsf{\xi }}_{max}-{\mathsf{\xi }}_{min}}{2cos\left(\alpha +\sigma \right)cos\phi }cos\left(\gamma \left(t\right)+\sigma \right),$$

(33)

$$\sigma =\underset{\sigma }{\arg \min }{∥\begin{array}{c}\mathsf{\Gamma }\left(\sigma \right)\end{array}∥}_1,$$

(34)

$$\mathsf{\Gamma }\left(\sigma \right)=R\left(\sigma \right)-\tilde{R},$$

(35)

$$R\left(\sigma \right)=\left[\begin{array}{c}\mathsf{\xi }\left(\sigma |t_0\right),\mathsf{\xi }\left(\sigma |t_1\right),\cdots ,\mathsf{\xi }\left(\sigma |t_N\right)\end{array}\right],$$

(36)

where (36) can be used to search the half cone angle, and $\tilde{R}$ is the distance difference between two scattering points in the estimation micro-curve.

3.3.4. Cone Height and Base Radius Estimation

The cone height and the base radius can be obtained according to the properties of the geometric parameters,

$$\left\{\begin{array}{c}H=Scos\left(\sigma \right),\hfill \\ r=Ssin\left(\sigma \right).\hfill \end{array}\right.$$

(37)

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

	Parameters	Values
warhead	cone height H	0.96 m
	distance between the mass center and the cone top h	0.64 m
	base radius of cone r	0.25 m
	spin frequency $f_s$	4 Hz
	conic rotation frequency $f_z$	2 Hz
	precession angle $\phi$	${10}^{°}$
radar	carrier Frequency $f_c$	10 GHz
	bandwidth B	2 GHz
	pulse width $t_p$	10 us
	pulse repetition period $prf$	1 kHz
	dwell time T	1 s

.

4.1. Scattering Point Separation

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.

4.2. 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 cone-top and cone-bottom scattering point exhibit nonlinear and non-stationary characteristics. We adopt the original EMD method and the EEMDCAN method to decompose the non-stationary 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 $-\pi$ or $+\pi$ 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 $\pi$ or $-\pi$ 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 EEMDCAN 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.

4.3. 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. The error rate (%) comparison of the EEMDCAN with EMD methods under different SNR conditions.

SNR (dB)	Methods	P	$\mathit{\theta }$	$\mathit{\sigma }$	s	H	r
30	EEMDCAN	0	0.88	0.85	1.69	1.34	0.79
	EMD	0	2.63	3.46	2.67	2.36	2.80
25	EEMDCAN	0	0.90	0.86	1.69	1.75	0.84
	EMD	0	2.65	3.38	2.80	2.80	3.15
20	EEMDCAN	0	0.89	0.83	1.67	1.72	0.83
	EMD	0	2.64	3.71	3.74	3.18	3.23
15	EEMDCAN	0	0.90	0.88	0.17	1.74	0.81
	EMD	0	5.90	5.32	3.91	4.04	3.68
10	EEMDCAN	0	0.91	0.79	1.64	1.34	0.84
	EMD	0	7.50	5.98	3.09	4.31	3.76
5	EEMDCAN	0	1.15	3.78	0.10	0.16	3.12
	EMD	0	8.23	5.87	5.45	6.45	4.84
0	EEMDCAN	0	1.41	5.40	2.55	2.94	1.54
	EMD	0	10.65	12.34	9.91	9.41	11.34
−5	EEMDCAN	0	2.81	30.76	10.62	12.76	15.39
	EMD	0	11.20	40.35	13.35	16.13	16.19
−10	EEMDCAN	0	2.21	33.56	8.01	10.29	21.15
	EMD	0	11.34	42.46	23.28	21.11	24.32

, 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

$$RMSE=\sqrt{\frac{1}{M}\sum _{m=1}^M{(r\left(m\right)-\widehat{s}\left(m\right))}^2},$$

(38)

where $r\left(m\right)$ and $\widehat{s}\left(m\right)$ 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.

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