An Improved Phase-Derived Range Method Based on High-Order Multi-Frame Track-Before-Detect for Warhead Detection

: It is crucial for a ballistic missile defense system to discriminate the true warhead from decoys. Although a decoy has a similar shape to the warhead, it is believed that the true warhead can be separated by its micro-Doppler features introduced by the precession and nutation. As is well known, the accuracy of the phase-derived range method, to extract micro-Doppler curves, can reach sub-wavelength. However, it suffers from an inefﬁciency of energy integration and high computational costs. In this paper, a novel phase-derived range method, using high-order multi-frame track-before-detect is proposed for micro-Doppler curve extraction under a low signal-to-noise ratio (SNR). First, the sinusoidal micro-Doppler range sequence is treated as the state, and the dynamic model is described as a Markov chain to obtain the envelopes and then the ambiguous phases. Instead of processing the whole frames, the proposed method only processes the latest frame at an arbitrary given time, which reduces the computational costs. Then, the correlation of all pairs of adjacent pulses is calculated along the slow time dimension to ﬁnd the number of cells that the point scatterer crosses, which can be further used in phase unwrapping. Finally, the phase-derived range method is employed to get the micro-Doppler curves. Simulation results show that the proposed method is capable of extracting the micro-Doppler curves with sub-wavelength accuracy, even if SNR = − 15 dB, with a lower computational cost.


Introduction
Discriminating warheads from decoys is one of the most significant and challenging tasks in a ballistic missile defense system [1][2][3]. Generally, decoys have a similar shape to the true warhead and are released in the mid-course stage to deceive the missile defense system. Fortunately, the warhead is usually heavier than the decoys, and it has micro-motions, such as spin and precession to keep a stable attitude in the mid-course phase, while the micro-motion forms of the decoy are mostly tumbling and large-scale nutation [4,5]. When using radar to observe the warheads and decoys, the different micro-motions will introduce different micro-Doppler features in the received echoes. This makes it possible to separate them from each other [6][7][8].
In the past few decades, various recognition methods have been proposed based on the micro-Doppler extraction, which can be roughly divided into five categories: radar cross-section (RCS) [4,9,10], high-resolution range profile (HRRP) [2,11], inverse synthetic aperture radar images [12][13][14], backscatter coefficient [15,16], and orbit information [17,18]. The RCS sequence is easy to obtain, and contains rich information about the target, such as the size, structure, surface material, and posture of the target [4]. However, RCS is affected by various factors, such as target material, target shape, etc. As a result, the micro-Doppler features extraction based on the RCS sequence may not be reliable enough [7]. Compared with RCS sequences, HRRP and inverse synthetic aperture radar images can obtain refined structure information of the target, which provides favorable conditions for high-precision micro-Doppler estimation [19,20]. However, with the development of the invade technique, the shape, size, and coating of the decoys are almost the same as the warheads, and it is hard for the methods based on HRRP and inverse synthetic aperture radar images to distinguish them, only according to their imaging results [2,4]. Note that, in the mid-course phase, motions of warheads and decoys are almost not affected by air resistance, so their translational motions are essentially the same. Consequently, using the orbit information to separate the two targets is difficult [2].
Recently, the envelope association method from the HRRP sequence for micro-Doppler curve extraction via wideband radar has attracted increasing attention [20][21][22]. In [20], the authors proposed a phase-derived range method for micro-Doppler parameter extraction. A one-order keystone transform is employed to fit the micro-Doppler envelope, which can achieve the straightening of the range migration without parameters. The one-order keystone transform is linear, but the warhead micro-Doppler envelope is sine-like and nonlinear, and directly using the one-order Keystone transform cannot fit the envelope well. An improved second-order keystone transform was also developed in [23]. It can remove the first-order or second-order coupling, but suffers from complex phase residual terms.
The method based on extended Kalman filter (EKF) was introduced to deal with the non-linear motion in [24]. It first employs the Taylor expansion to obtain the local linearity of the non-linear motion and then uses the Kalman filter to estimate the motion state of the target. Note that when dealing with the motion with a high non-linearity degree, EKF will be degraded by the truncation error. Moreover, EKF is sensitive to the SNR, and low SNR will further deteriorate its performance. In [21], the authors proposed a modified Kalman filter for the envelope extraction of the cone-shaped targets, which is more robust and precise than the EKF. However, the performance of the modified Kalman filter is limited due to unknown noise covariance under low SNR conditions.
In recent years, methods based on track-before-detect (TBD) that directly process the raw data or the thresholded data with a small value, have been proposed for target detection under low SNR conditions [17,25]. Generally, such TBD methods can be divided into two categories: the single-frame and multiple-frames TBD [17]. Similar to the Kalman filter, the single-frame TBD methods estimate the target states by predicting and updating the posterior probability density, such as the random finite set algorithms [26,27], particle filter [28], and histogram probabilistic multi-hypothesis tracker [29]. Comparatively, multiple-frames TBD methods normally use a sliding window to jointly process multiple data frames at each measurement time, and are capable of achieving superior performance for target detection under low SNR conditions, at the expense of high computing complexity [17], such as dynamic programming TBD [30], Hough transform [31], velocity matched filtering [25], and maximum likelihood probabilistic data association [32].
In this paper, an improved phase-derived range method with a high-order multiframe TBD (HO-MF-TBD) is presented to extract the micro-Doppler curves under low SNR circumstances. First, extracting the micro-Doppler sequence from the HRRP, and the micro-Doppler ranges form a Markov chain. Then, the HO-MF-TBD is employed to estimate the states of the Markov chain, and only the latest frame is processed at an arbitrary given time. The micro-Doppler curve envelopes and the ambiguous phase can be obtained. At the same time, one can calculate the correlation of two adjacent pulses of each segment to get the moving distance, or range migration units, of the target, which is reflected by the number of range cells it crosses on the HRRP from one pulse to the other. Then, the unambiguous Doppler phase can be easily obtained from the ambiguous Doppler phase and the range migration units. Finally, the micro-Doppler curve extraction can be done by the commonly-used phase-derived range method. The main contributions of this paper include: (1) The HO-MF-TBD method is proposed for envelope extraction, with model robustness and lower computational costs under low SNR conditions. (2) A new phase unwrapping technique, which calculates the moving distance from two adjacent frames, is proposed to improve the ranging accuracy.
The remainder of this paper is organized as follows. Section 2 gives the geometric model and signal model. Section 3 explains the micro-Doppler envelope extraction using the HO-MF-TBD algorithm. A new phase unwrapping method is presented in Section 4. In Section 5, several simulation results are presented to verify the proposed method. Concluding remarks are provided in Section 6.

Geometric Model and Signal Model
In this section, the geometric and signal model is introduced. To ease the reading, the main notations are listed in Table 1.

Notations Meanings
number of total frames of patch processing s m measurement matrix in frame m z m state matrix in frame m S 1:K measurement set Z 1:K state set θ m ambiguous phase Ω m unambiguous phase

Geometric Model
The geometric model used in this paper is based on the model in [2,7,22]. As shown in Figure 1a, the radar is located at the origin of the coordinate system (U, W, V), and remains stationary during observation. The target coordinate system (x, y, z) takes the target centered as the coordinate origin, and the target is accompanied by cone rotation in flight. A smooth cone warhead contains three point scatterers, namely p1, p2, and p3, respectively. p1 is fixed, while p2 and p3 change their positions with the changes of radar line of sight. The theoretical expression of the micro-Doppler curve from p1, p2, and p3 to radar is where R 0 , H, d, r, ω c , and φ represent the distance from warhead centered to radar, the warhead height, the distance from warhead centered to warhead top, the warhead bottom radius, the angular velocity, and the initial phase, respectively. α and β denote the elevation angle of radar line of sight and precession angle, respectively. Generally, p3 is occluded because of the occlusion effect [21].
Radar R a n g

Signal Model
The linear frequency modulation signal is widely used in wideband radar systems [11]. The transmitted linear frequency modulation signal can be expressed as wheret, t m , T p , f c and γ denote the fast time, the slow time, the pulse duration, the carrier frequency, and the slope frequency modulation, respectively. Usually, the reference signal is used for de-chirp processing, which is a replica of the transmitted signal, and it can be expressed as where t re f = 2R re f /c denotes the time delay caused by the speed of light passing through the reference distance, R re f represents the reference distance, and T re f represents the length of the receiving window, which is larger than T p , so that the reference signal can cover the transmitted signal. Then, the echo signal from a point scatterer can be expressed as where ρ i (t m ) represents the amplitude containing the sinc function of point scatterer i, t i denotes the time caused by the motion change of point scatterer i, and n(t m ) ∼ CN(0, σ 2 u ) is the receive noise with Gaussian distribution, the mean value is 0 and the variance is σ 2 u . After de-chirp processing, the signal can be expressed as where ρ i (t m ) represents the amplitude containing the sinc function after de-chirp processing of point scatterer i. Generally, the pulse width T p of the rectangular signal emitted by the radar is very narrow, so the distance change caused by the micro-movement of the target within T p can be ignored, and the translation of the target can be approximated as a uniform motion [11]. Therefore, the radial range from the target to the radar can be expressed as where R o , v(t m ) and r i (t m ) represent the radial range of the target centered, radial velocity, and range variation induced by micro-motion, respectively. Substituting (7) into (6) and performing fast Fourier transform on (6) can get the target's HRRP. Reference (6) shows that the speed takes a great influence on the target's HRRP, especially for high-speed ballistic targets. Therefore, speed compensation must be performed to obtain a usable HRRP image. Fortunately, the speed measurement accuracy of modern radars has already met the requirements of speed compensation [33]. Therefore, this article will not discuss the impact of target speed on HRRP. The target model is shown in Figure 1a. The HRRP of the target can be expressed as where I,ρ i (t m ), λ and n H (t m ) represent the number of point scatterers, the amplitude of the target contain sinc function, the wavelength of the carrier, and the complex white Gaussian noise after fast Fourier transform, respectively. Generally, the attitude of the warhead is relatively stable in the mid-course phase, and the micro-Doppler of the wideband radar is easily extracted from the bulk motion [21]. Then the micro-Doppler model of the i-th point scatterer can be expressed as [11] r i (t m ) = a i cos(w i t m + ϕ i ), (9) where a i , w i and ϕ i represent the amplitude, angular velocity, and initial phase of the i-th micro-Doppler point scatterer, respectively. The micro-Doppler of the target is a sin function, as shown in the Formula (9), which is a high-order motion. The echo is digitized, which is arranged in a matrix of M × N grid of cells based on the resolution of radar, as shown in Figure 1b, where M and N represent the matrix length of the slow-time dimension and the fast-time dimension, respectively. Then, the pulse-echo signal of the slow-time dimension of the i-th point scatterer of Formula (8) can be rewritten as where T r denotes the pulse repetition interval, t m = mT r . Figure 1c respects that the point scatterer across several range cells. Therefore, the HO-MF-TBD is employed to associate the envelope of the point scatterer optimally.

Micro-Doppler Envelope Extraction Using HO-MF-TBD
This section focuses on describing the HO-MF-TBD method of micro-Doppler envelope extraction.

Measurement Model
In Section 2, the HRRP sequence of each point scatterer is divided into M × N cells based on the resolution of radar, namely ∆ t × ∆ r , where ∆ t and ∆ r denote the resolutions of M and N directions, respectively. Take M × N as measurement values, and the mea-surement intensity z m (n t , n r ) recorded in quantized cell (n t , n r ) in the mth frame can be formulated as z m (n t , n r ) = W(n t , n r , m), where n t = 1, 2, . . . , M, and n r = 1, 2, . . . , N, W(n t , n r , m) represents the additive white Gaussian noise with independently and identically distributed. The measurements of K frames in a multi-frame batch are denoted as In this paper, an effective energy integration is presented for weak target detection with low thresholds. First, calculate the amplitude of all K frames in a processing batch, and record the cells with amplitude greater than 0. The amplitudes are assumed to be independently distributed in the map depending on the distribution of targets and noise. Therefore, the measurement can be expressed as where A m denotes the amplitude of point scatterer, n

Target Model
The target state s m in the mth frame is assumed to obey the first-order Markov process, which is express as s m+1 |s m ∼ p(s m+1 |s m ), s m ∈ D where p represents the probability density function, and D denotes D-dimensional different state. According to the Bayesian estimation theory, the state transition equation is where I D represents D × D identity matrix, and ⊗ the Kronecker product.

HO-MF-TBD Method
According to the Bayesian estimation theory, the measurement model describing the relationship between the state s m and the measurements z m , where h m is a function of the state s m and the measured noise n m at the frame m. The track and measurement data of a target in a batch of K frames can be expressed as S 1:K = (s 1 , s 2 , · · · , s K ) and Z 1:K = (z 1 , z 2 , · · · , z K ), respectively. Based on the Bayesian estimation theory, the estimator of the state sequence can be obtained by the posterior probability density function of the state sequence P(S 1:K |Z 1:K ), and I(s 1 |Z 1:K ) = L(z 1 |s 1 ), where I(s m |Z 1:K ) represents the value function corresponding to the target state s m , G s m (·) is used to store the state transition relationship between frames, and L(z m |s m ) represents the measurement plane likelihood function where τ(s m−1 ) represents the state range of the transition of the state within one frame time. The computational complexity of (21) is O(| D |τ(K − 1), it increases linearly with the number of scans K, and D is a D-dimensional state space. The calculation amount is reduced by limiting the value range of τ(s m−1 ), take τ(s m−1 ) equal to 5 as an example, the search range and accumulation method are shown in Figure 2. P(s m |s m−1 ) represents the cost function of the state transition. The estimated envelope iŝ where V DT is a pre-set detection threshold, and if I(Ŝ K |Z 1:K ) > V DT , then for m = K − 1, · · · , 2, 1,Ŝ In Figure 2, the dotted arrow represents the search range, and the solid arrow represents the points selected from the search range to estimate envelopeŜ K . However,Ŝ K may deviate from the true path due to the influence of noise. Therefore, median filtering is introduced to reduce the error caused by noise.Ŝ K is recorded asŜ K after median filtering.
The ambiguous micro-Doppler phase can be obtained according to the envelope where Im(Ŝ K ) and Re(Ŝ K ) denote the imaginary and real parts ofŜ K . The unambiguous micro-Doppler phase can be obtained by the phase unwrapping method, which will be introduced in Section 4.
· · · · · · · · · · 1 · · · · · · · · M · · · · · · · N Slow time Fast time Figure 2. Take τ = 5 as an example to illustrate the envelope extraction using the HO-MF-TBD method. The dotted arrow represents the search range, and the solid arrow represents the points selected from the search range to estimate the envelope.

Phase Unwrapping
The mth ambiguous micro-Doppler phase obtained in Section 3 is defined as θ m , and the unambiguous phase Ω m can be expressed as where k m is an integer to ensure 0 ≤ Ω m − 2k m π + ε m ≤ 2π, and ε m is the measurement error. The D m and d m are defined as the difference between the (m + 1)th and mth unambiguous and ambiguous micro-Doppler phase, respectively. Then where N m = k m+1 − k m . The real difference between the mth and (m + 1)th echo signal pulse is D m = d m + 2πN m . So the key point to get the k m . It is known that the similarity between the two signals is measured by the correlation coefficient. Thus, the correlation of any adjacent pulse pairs is calculated along the slow time dimension to find the number of range cells is crossed by the point scatterer. The correlation between the mth and (m + 1)th pulse is expressed as where f ix(∆N step ) denotes range cells the point scatterer across, and where λ is the wavelength of the radar carrier, and ∆x = c/2B respects resolution in the fast time direction. From (27) to (30), it can yield the N m . The real radial range can be obtained where 1 ≤ m ≤ n − 1. The flowchart for the implementation of the proposed method is shown in Figure 3.

Simulations and Discussion
Some Monte Carlo simulation experiments are carried out to illustrate the performance of the proposed method. The height and radius of the metal cone are 0.96 and 0.25 m, respectively. The spin frequency, precession angle, pulse repetition period (PRF), and rotation frequency are 4 Hz, 10 degrees, 1 KHz, and 2 Hz, respectively. The radar works at the carrier frequency of 10 GHz and with a bandwidth of 2 GB. The parameters are shown in Table 2, and the warhead is shown in Figure 1a. Figure 4 shows the HRRP sequence according to the parameters in Table 2 with the SNR at −15 dB. Figure 4 shows that the target occupies some range bins due to the wideband radar and warhead volume.

Micro-Doppler Envelope Extraction Using HO-MF-TBD Method
For each point scatterer, HO-MF-TBD is used to associate the envelope. The envelope extraction results are shown in Figure 5a. Figure 5a shows that there are errors in some range cells. Therefore, a median filter is used for smoothing to eliminate the accumulated error. The smoothed envelopes are shown in Figure 5b.
Based on the envelopes, the ambiguous Doppler phase is obtained directly, and presented in Figure 6a,d. The phase hopping occurs where the range changing is greater than one range bin in Figure 6d. Thus, phase compensation is required to be performed before the phase unwrapping to guarantee the unwrapping accuracy. The phase difference between two adjacent phases is calculated, and take the point where the difference is greater than 0.5 as phase hopping, as depicted in Figure 6b,e. Then, additional compensation π or −π is carried out to obtain the smooth ambiguous Doppler phase, which is shown in Figure 6c,f.

Phase Unwrapping
Correlation is the degree of similarity between two sequences. The correlation of any adjacent pulse pairs is calculated along the slow time dimension to find the number of range bins crossed by the point scatterer. For example, Figure 7a-c show the relevant results of the 45th and 46th, 46th and 47th, and 47th and 48th frames, respectively. Figure 7a-c have the highest correlation at the 480th, 482nd, and 482nd range cells, respectively, which means that the point scatterer moves from the 480th range cell to the 482nd range cell from frame the 46th to frame 47th, crossing two range cells. However, the point scatterer does not cross the range cell from frame 47 to frame 48. The number of range cells that the point scatterer crosses can be further used in phase unwrapping. Then, the micro-Doppler curve can be estimated accurately by the phase-derived range method. The comparison between the proposed method and the Keystone-based method is shown in Figure 8. It can be seen that the proposed algorithm is better than the keystonebased method. The main reason is that the approximate fitting is required to fit the high-order micro-Doppler model when using the keystone-based algorithm.  To further evaluate the performance of the proposed method, the RMSEs and the time cost comparisons among the proposed method, EKF method, dynamic programming TBD (DP-TBD) method, and the empirical mode decomposition (EMD) method under different SNRs are simulated. The RMSE is given by where Y is the number of Monte Carlo trials, x m,y andx m,y represent the actual target position and the corresponding estimated position at time m in the yth (1 ≤ y ≤ Y) Monte Carlo trial. d(x m,y , x m,y ) denotes the Euclidean distance betweenx m,y and x m,y . The time cost is evaluated by the time running each of the four algorithms on the same computer.
Only the time for envelope extraction is calculated. In this experiment, the four algorithms mentioned above work in a unified radar system. The warhead parameters and the radar parameters are shown in Table 2. The SNRs are set as −15 to 10 dB, with 5 dB step lengths, and 100 Monte Carlo simulations are performed for each SNR. Only the envelope extraction are using different algorithms, while other steps are kept the same. The simulation results of RMSEs and time cost are shown in Tables 3 and 4, respectively. Note that the time cost listed in Table 4 may be affected by CPU memory occupation. However, it can generally reflect the complexity of an algorithm. From the wavelength calculation formula, the wavelength can calculated as λ = c/ f = 3 × 10 8 /10 × 10 9 = 0.03 m. Table 3 illustrates that the RMSE of the proposed method, the EKF method, the traditional DP-TBD method, and the EMD method almost reach the subwavelength level accuracy using the phase-derived range method. Generally, the RMSE of the proposed algorithm is better than the EKF method, the traditional DP-TBD method, and the EMD method under low SNR conditions. The main reason is that according to the motion law, the envelope of the current frame is usually distributed near the envelope of the previous frame. Therefore, the error caused by noise can be eliminated by limiting the range of dynamic programming.
The EKF method takes the least time, while it has the worst RMSE performance. This is because of the influence of truncation error and the unknown noise covariance under low SNR conditions. The RMSEs of the traditional DP-TBD method are similar to the proposed method, but among the four comparison algorithms, the traditional DP-TBD algorithm takes the longest time because the implementation process requires dynamic programming. Therefore, the classical DP-TBD algorithm needs to be optimized to meet the high real-time requirements of an anti-missile system.
In addition, the EMD method can decompose the non-stationary signal into a series of inherent modal functions, which can be used to reconstruct the non-stationary signal to form envelope correlation for envelope extraction [34]. However, it heavily depends on the extreme points of non-stationary signals. Therefore, as Table 3 shows, the performance of the EMD method becomes poor under low SNR conditions.

Conclusions
In this paper, an improved phase-derived range method based on high-order multiframe TBD is proposed to conquer the difficulties in distinguishing the warheads from decoys. Specifically, the micro-Doppler range sequence is regarded as a state, and the HO-MF-TBD is used to obtain the envelope. In addition, a new phase unwrapping method is proposed to improve the accuracy of micro-Doppler curves extraction. This method reduces the computational complexity by limiting the scope of dynamic programming. The proposed method can estimate the micro-Doppler curves with high accuracy at a sub-wavelength level even when the SNR is lower than −15 dB. Extensive simulation results validate the anti-noise performance of the proposed method. The directions of future research could be in distinguishing the electronic false targets.
Author Contributions: Conceptualization, N.Z.; methodology, N.Z. and S.X.; software, N.Z.; validation, N.Z. and X.F.; investigation, N.Z. and S.X.; data curation, N.Z. and S.X.; writing-original draft preparation, N.Z. and W.W.; writing-review and editing C.L., J.H., S.X. and Z.C.; supervision, S.X. and Z.C. All authors contributed to writing the final manuscript. All authors have read and agreed to the published version of the manuscript.