Radar Maneuvering Target Motion Parameter Estimation Based on Hough Transform and Polynomial Chirplet Transform

This paper considers the problem of estimating the motion parameters of a low-observable high-speed maneuvering target with high-order motions, where the target’s complex motions will result in the range migration (RM) and Doppler frequency migration (DFM) effects within the observation time interval, which will significantly affect the radar performance. Based on the Hough transform (HT) and the polynomial Chirplet transform (PCT), a novel method for the motion parameter estimation, which we coin as HPCT, is proposed in this paper. It first employs the trajectory feature extracted from the image of range profiles by the HT to correct the first-order RM and resolve the Doppler ambiguity. Then, it compensates the residual range curvature and the DFM with the motion parameters estimated from the time-frequency ridge, which is extracted from the highly concentrated time-frequency representation (TFR) generated by the PCT. After that, the energy of radar returns can be coherently integrated along the target’s moving trajectory. In addition, we evaluate the parameter estimation performance, integration performance and computational cost of the proposed method via theoretical analysis or numerical experiments. The results show that in terms of the anti-noise ability, cross-term interference avoidance, measurement accuracy and computational complexity, the HPCT is superior to some common methods for the motion parameter estimation. Finally, we extend it to the scenario of multiple targets with different motion orders.


I. INTRODUCTION A. MOTIVATION AND RELATED WORKS
Under low signal-to-noise ratio (SNR) environment, robust and effective motion parameter estimation of moving target is one of the fundamental and important subjects in the field of radar signal processing, since the estimated motion information is very meaningful for the target detection, tracking, imaging and identification [1]- [11]. It is well known that the SNR of radar returns can be improved significantly by means of the coherent integration among multiple pulses [1], which is extremely helpful for the motion parameter estimation. However, for a maneuvering target, the target's complex motion usually leads to the range migration (RM) The associate editor coordinating the review of this manuscript and approving it for publication was Wei Liu . and Doppler frequency migration (DFM) [8]- [13]. Especially for the high-speed maneuvering target with complex motions, these two effects may inevitably cause the energy of radar returns spread in both the time and Doppler frequency domains. If these two effects cannot be compensated precisely and simultaneously, then the radar performance will be seriously affected because the echo energy could not be efficiently exploited to improve the SNR of radar returns.
In the past decade, the issue of eliminating the effects of RM and DFM simultaneously has attracted considerable attention [8]- [34], where the target's movement within the observation time interval is generally modeled as a polynomial function. In order to get better fidelity, the order of the motion model discussed in different literatures is getting higher in recent years. It improves from the uniform velocity model [14]- [18] to the uniform acceleration model [19]- [23], then to the uniform jerk model, and even to the arbitrary highorder motion model. Focus on the high-speed maneuvering target whose motion order is no less than jerk, the relevant studies could be roughly categorized into the following four categories.
The first category utilizes the amplitude information or the similarity/correlation between the envelope range profiles to measure the target's motion parameters, e.g., the range centroid method [24], envelope cross-correlation method [3]. This kind of method is straightforward and computationally efficient. However, in low SNR situation it could not accurately estimate the target's motion parameters because the abovementioned similarity/correlation would be corrupted by the huge noise and interference.
The second category utilizes the phase characteristic of moving target's echoes and the coupling relationship among target's motion parameters, instantaneous range and echoes' phase to reduce (or iteratively reduce) the order of RM and DFM and then to estimate the target's motion parameters. The existing method includes adjacent across correlation function (ACCF) method [25], [26], symmetric autocorrelation function method [23], [27] sequence/time-reversing transform (SRT/TRT) [28], [29], frequency axis reversal transform (FRT) [30], etc. In theory, this kind of method does not need a searching procedure for the motion parameters. Hence, its computational burden is comparatively low. However, it requires that the signal is strict symmetry either in the time domain or in the frequency domain [23] so that it would not work well under low SNR environment. Furthermore, it would inevitably introduce plenty of cross-term interferences that are troublesome for the maneuvering target detection and motion parameter estimation in realistic and practical application.
The third category utilizes the image feature extracted from the range profile image and/or the time-frequency image via the image processing and feature extraction techniques to estimate the target's motion parameters. For example, B. D. Carlson et al. [31] realized long-time non-coherent integration via Hough transform (HT) in a high-dimensional data space. However, they did not consider the DFM effect of the maneuvering target. In order to solve this problem and realize the long-time coherent integration, X. Li et al. [32] combined the Radon transform (RT) with the generalized de-chirp process (GDP) to estimate the target's motion parameters. However, its computational cost is huge because an exhaustive searching procedure is needed in GDP. P. HUANG et al [33] combined the generalized high-order ambiguity function (GHAF) with the HT to detect both the range walk trajectory and the time-frequency trajectory. However, the cross-term interferences will inevitably appear due to the nonlinear transform property [33] of GHAF.
The fourth category achieves the long-time coherent integration along the target's moving trajectory via jointly searching in the target's motion parameter space, e.g., the generalized Radon Fourier transform (GRFT) [13], generalized Keystone transform-generalized de-chirp process (GKT-GDP) [12] and Radon-S transform [34]. This kind of method can obtain desirable estimation and integration performances. However, for a maneuvering target with complex motions, its main drawback is that the exhaustive searching in high-dimensional parameter space is often computationally very expensive, which may limit its practical applications.

B. SUMMARY OF THE CONTRIBUTIONS
This study considers the problem of motion parameter estimation for a low-observable high-speed maneuvering target with high-order motions. The contributions of this work can be summarized as follows: • Based on the Hough transform (HT) [31] and the polynomial Chirplet transform (PCT) [35], a novel motion parameters estimation method, which we coin as HPCT, is proposed in this paper. It first adopts the HT to detect the target's moving trajectory from the image of range profiles and measure the target's average velocity within the observation time duration, which is used for the first order range migration (a.k.a. range walk) compensation (FRMC) and the Doppler ambiguity resolution (DAR). Then, it employs the PCT to obtain the time-frequency representation (TFR) of the signal that is extracted from the most representative range sampling cell, and estimates the remaining motion parameters, including the residual velocity, acceleration, jerk, and so on, with the time-frequency ridge feature extracted from the TFR. After that, it re-corrects the RM and re-estimates the initial slant range from the radar platform to the target. Finally, it coherently integrates the energy of radar returns along the target's moving trajectory that is rebuilt by the estimated motion parameters.
• Then, the parameter estimation performance, integration performance and computational complexity of the proposed method are investigated via theoretical analysis or numerical experiments. The results demonstrate that the proposed method has many advantages when compared with some common methods for motion parameter estimation, such as excellent anti-noise performance, free from cross-term interference, high measurement accuracy and low computational complexity. Therefore, it seems to be promising and attractive for the motion parameter estimation of the high-speed maneuvering target under low SNR environment.
• After that, the proposed method is extended to the scenario of multiple targets, and its performance is evaluated with numerical experiment. The results demonstrate that the proposed method can estimate the motion parameters of the maneuvering targets with different motion order accurately and integrate the signal energy of each target along its moving trajectory, which is helpful to the target detection.

C. PAPER ORGANIZATION
The remainder of this paper is organized as follows. The transmitted and received signal models are briefly introduced in Section II. In Section III, the motion parameter estimation method via HPCT is proposed. Section IV evaluates the performance of the proposed method via several numerical experiments. Finally, we draw the conclusions and present its possible future research direction in Section V. Beside, although every mathematical notation is defined in the paper prior to its use, for the sake of convenience, the major symbols used in this paper are summarized in TABLE 1.

II. SIGNAL MODEL
Suppose that the radar transmits a train of identical linear frequency modulated pulses [1] for simplicity, i.e., where rect (x) = 1 for |x| ≤ 1/2 0 otherwise, m = 0, 1, . . . , M − 1 and M is the number of transmitted pulses during the observation time interval T obs , T p is the pulse width, T r denotes the pulse repetition interval, f c is the carrier frequency, K r = B T p is the frequency modulated rate and B represents the signal bandwidth. Suppose there is a moving target in the radar field of view, and the geometry of the radar and a moving target can be depicted by FIGURE 1. The target is moving from its initial location at (a 0 , 0) with a heading angle [3] of θ, where the X coordinate is along with the radar-target line of sight. At time instant t, the instantaneous slant range from the radar platform to the target is R (t). According to the Weierstrass approximation principle [36] and the ''stop-andhop'' assumption [1], R (t) can be approximated with a Q order discrete-time polynomial function R (t m ) within the finite observation time interval, i.e., where t m = mT r represents the slow-time, a 0 is the initial slant range between the radar and the target, a 1 , a 2 . . . , a Q are the projected motion parameters of target along the 35180 VOLUME 9, 2021 radar-target line of sight (LOS), such as the radial velocity, acceleration, jerk, etc.
In accordance with the pulse repetition interval T r and the range sampling interval t s (range sampling rate f s = 1 t s ), the radar echo signals after quadrature demodulation within the observation time interval T obs and the range gate r g can be discretized and then be rearranged as a two-dimensional M × N matrix (ignoring the noise components for convenience): where m = 0, 1, . . . , M − 1, n = −N 2, −N 2 + 1, . . . , N 2 − 1, M is the number of received pulses, N is the number of range sampling bins within the range gate r g , A r is the amplitude of received baseband echo signals, t n is the discrete version of the fast-time t = t −mT r −2R c c with the range sampling interval t s , λ = c f c is the carrier wavelength, R (t m ) = R (t m ) − R c denotes the relative range from the scene range center R c , which corresponds to the center of range gate, to the target in the direction of the LOS.
After pulse compression (a.k.a. range compression), the compressed signal in the slow time-range frequency (t m − f r ) domain can be expressed as where f r is the range frequency corresponding to the fast time t n and A 1 is the amplitude of S c (t m , f r ). By simply taking the inverse Fourier transform (IFT) of (4) with respect to the range frequency f r , the so-called complex range profiles in the t m − t n domain can be obtained: where sin c (x) = sin (πx) (πx), A 1 is the amplitude of the complex range profiles. According to (5), the envelope peak of s c (t m , t n ) locates at t n = 2 R (t m ) c and the instantaneous Doppler frequency at time t m is f d (t m ) = −2 λ · d R (t m ) dt m . Obviously, both the location of envelope peak and the instantaneous Doppler frequency vary with the slow-time t m , which are caused by the target's movement. If the envelope peak of range profile moves beyond a range resolution cell, the effect of range migration occurs. If the Doppler frequency varies beyond a Doppler frequency resolution cell, the effect of Doppler frequency migration occurs. As a result, the energy of radar returns may spread across both the range coordinate and the Doppler coordinate, which would deteriorate the radar performance seriously.

III. TARGET MOTION PARAMETERS ESTIMATION VIA HPCT
In this section, a method based on the Hough transform (HT) and the polynomial Chirplet transform (PCT) is proposed to estimate the motion parameters of a low-observable highspeed maneuvering target with high-order motions. Its block diagram is given in FIGURE 2.

A. TRAJECTORY DETECTION VIA HT
The Hough transform (HT) is an efficient feature detector for detecting lines, curves, ellipses, etc., in noisy images [37]. In principle, it is designed to detect the parameterized curves that are not required to be continuous or derivable [31]. These properties are very useful in detecting the target's moving trajectory under low SNR environment. For example, the prior knowledge of the target's motion can be easily introduced and the gaps in the moving trajectory caused by the outlier data or noise can be tolerated to some extent.
According to (2), the target moving trajectory can be approximately represented by a Q order polynomial. Although the generalized HT can detect this type of polynomial curves, its computational load will go up significantly as the value of Q is increased [31]. Fortunately, in the socalled coarse range migration compensation procedure [3], we only need to align the prominent scatterer into the same range resolution cell across the range profiles. Therefore, we can approximate the target's moving trajectory with a piecewise linear function, and adopt the standard HT to detect the trajectory to reduce the computational burden.
Before applying the HT to detect the target's moving trajectory, it is convenient to eliminate part of noise via a primary thresholding transform as in the following: with where max n (|s c (t m , t n )|) represents the maximum amplitude of the mth envelope range profile along the t n coordinate and VOLUME 9, 2021 A max denotes the sample mean of max n (|s c (t m , t n )|), λ A is a scaling factor used for adjusting the threshold. As we have discussed, the target's moving trajectory can be approximated by a piecewise linear function for the coarse range alignment. Each straight line in t m − t n space can be defined by the angle θ of its normal and the distance ρ from the origin to this line. The HT maps all the threshold-crossing points in the t m − t n space into a set of sinusoids in the ρ − θ space by [31] If a straight line does exist in t m − t n space, it will be represented in ρ − θ space as an intersection point of all of the corresponding mapped sinusoids. In contrast, the measurements attributed to the noise are randomly distributed over the entire ρ − θ space. Therefore, a two-dimensional accumulator array can be used to detect the lines with the secondary threshold [31]. In addition, a masking operation based on the prior knowledge about the target's motion can also be applied to further reduce the clutter components. With all this in place, a point denoted by (ρ c , θ c ) in ρ − θ space with the strongest straight line return can be extracted as the medoid (the most representative point) [38], [39] to estimate the target's average velocity in T obs , which can be expressed in terms of θ c :â

B. FRMC AND DAR
In order to correct the first-order range migration (a.k.a. range walk) caused by the target's velocity, with the estimated average velocityâ 1 we can construct a compensation function in the slow time-range frequency (t m , f r ) domain as follows: After multiplying (4) with (9) and performing the IFT along the range frequency f r domain, we have obviously, loc (t m ) is the location of the envelope peak of the mth range profile after FRMC. It can be seen from (10) that the first order range migration has been eliminated. Nevertheless, the high-order effects of the range and Doppler migrations are still present. If the variation of loc (t m ) within the observation time interval is less than the range resolution, then we can say that all the range profiles are almost aligned, otherwise the value of T obs may need to be adjusted to reduce the effect of the high-order motions. The index of range sampling cell where the target locates at can be estimated bŷ with where p r (n) can be understood as the probability of the target locates at the index n. It is worth pointing out that if the target is moving with a constant velocity within the observation time interval, we can directly compute the indexn with (ρ c , θ c ).
The signal at thenth range sampling cell along the slowtime t m coordinate can be approximated as Because the Doppler frequency of the high-speed target is usually higher than the pulse repetition frequency (PRF) when the low-PRF scheme is employed to increase the range unambiguous observation scope, the radar will be ambiguous (or extremely ambiguous) in Doppler. With the estimated average velocityâ 1 , a compensation function used for resolving the Doppler ambiguity can be constructed as follows: Multiplying (12) with (13), we have with where ε i a , i = 0, 1, . . . , Q, is the residual errors. Obviously, s d (t m ) is a polynomial frequency modulated signal. Its instantaneous frequency (IF) can be expressed as If the maximal value of |IF (t m )| is less than half of PRF, we can say that the Doppler ambiguity has been resolved.

C. TIME-FREQUENCY REPRESENTATION VIA PCT
The Polynomial Chirplet transform (PCT) [35], [40] is a parameterized time-frequency analysis method. Essentially, it approximates the time-frequency trajectory of a given signal with a polynomial kernel function. By selecting a set of appropriate characteristic parameters for the polynomial kernel, the PCT is able to produce a time-frequency representation (TFR) with excellent energy concentration. In addition, the PCT is free from the problem of cross-term interferences [35], [40]. In order to get at the time-frequency signatures of s d (t m ), motivated by the characteristics of PCT, we employ the PCT to obtain the TFR of s d (t m ), which can be stated as where R c (τ ) and S c t m , τ are the frequency rotating operator and the frequency shifting operator, respectively, κ c (t m ) represents the polynomial kernel of PCT with parameters c = c 2 , c 3 , . . . , c Q , and g σ τ − t m denotes the Gaussian window function centered at t m with time bandwidth σ .
According to (17) and [35], [40], if the kernel characteristic parameters c = c 2 , c 3 , . . . , c Q accurately match the residual errors ε a = ε 2 a , ε 3 a , · · · , ε Q a at all specific time instant, then the DFM effect caused by the target's complex motions (i.e. acceleration, jerk, and so on) can be eliminated completely. With this in hand, the PCT can generate a TFR with a superior energy concentration whose frequency resolution is only determined by the bandwidth of Gaussian window 1 σ . This is particularly helpful for extracting the time-frequency ridge feature from the TFR to estimate the residual errors ε i a , i = 1, . . . , Q. Therefore, the determination of proper kernel characteristic parameters is critical for the motion parameter estimation.

D. TIME-FREQUENCY RIDGE FEATURE EXTRACTION
As we have discussed, in order to generate the highly concentrated TFR, we need to determine a set of appropriate characteristic parameters for the polynomial kernel. From the generalized matched filtering point of view, if the number of the kernel characteristic parameters is comparatively small (such as 1 or 2) and the appropriate prior information about the target's motion can be used to limit the searching range, then the algorithms based on the brute-force exhaustive search paradigm [39] may be a good choice for finding the optimal set of kernel parameters. Nevertheless, for the maneuvering target with complex motions, these two requirements are usually not satisfied. As a consequence, this kind of algorithm is computationally very expensive, which may not be the best choice in a radar system for the requirement of real-time and online processing.
As we can see from (16), the instantaneous frequency (IF) of s d (t m ) is a polynomial function with unknown coefficients ε i a , i = 1, 2, . . . , Q. Because this function can be approximated by the time-frequency ridge extracted from the TFR, obviously, we can use it to estimate the residual errors ε i a via the polynomial regression [38] and then set the kernel characteristic parameters c. Based on the principle of time-frequency feature approximate [41], the procedure for extracting the time-frequency ridge feature and estimating the residual errors without searching in high-dimensional parameter space can be stated as follows: 1) Initialize the order of the approximating polynomial with a list Q ploy , the Gaussian window size with N win , the initial kernel characteristic parameters with c = (0, 0, . . . , 0), and set the maximum number of passes over the dataset (epochs) to N iter ; 2) Utilize (17) to generate the TFR of s d (t m ) with c and then extract the time-frequency ridge by locating the peaks of TFR at each time instant; 3) Fit the time-frequency ridge with a set of polynomial functions of order G i poly (G i poly ∈ Q ploy ) to determine the value of Q by comparing the corresponding mean squared errors (MSEs), and estimate the residual errors ε i a , i = 1, 2, · · · , Q; 4) If N iter is not reached, update c with the estimated residual errorsε a = ε 2 a ,ε 3 a , · · · ,ε Q a from step 3 and go to step 2, otherwise outputε i a , i = 1, 2, · · · , Q. The reason we initialize c = (0, 0, . . . , 0) is that there is no prior knowledge about ε a can be used to set the kernel characteristic parameters c at the beginning. It is also important to mention that the right choice of the hyper-parameter Q in step 3 is crucial to find a good balance between over-fitting and under-fitting [39]. Since the value of Q should be a relatively small positive integer for the high-speed maneuvering target in a finite observation time interval, the additional computational cost is relatively cheap.

E. REMAINING PARAMETERS ESTIMATION AND RMC
According to (15), the estimated motion parameters can be updated by The range migration compensation function can then be updated as With (19), after re-correcting the range migration, the range profiles can be rewritten as with (10), it can be seen that both the range walk and the range curvatures caused by the target's high-order motions have been corrected with the estimated motion parametersâ i , i = 1, 2, . . . , Q. Substitute s ra (t m , t n ) with s ra (t m , t n ) in (11), the index of range sampling cell where the target locates at can be re-estimated aŝ n , and the initial slant range between the radar platform and the target can be calculated bŷ

F. DFMC AND ENERGY INTEGRATION
Using the estimated motion parametersâ i , i = 0, 1, . . . , Q, the Doppler frequency migration compensation function can also be re-constructed as With the re-estimatedn , we extract the signal at then th range sampling cell along the slow-time domain and multiply it with H d (t m ), the Doppler frequency migration compensated signal can be written as It can be seen from (23) that the high-order DFMs have been eliminated with the help ofâ i . After that, the energy of radar returns can be coherently integrated along the target's moving trajectory by performing the Fourier transform (FT) on s d (t m ): In addition, if we need the two-dimension (2-D) image of target, it is also convenient to generate it with the estimated motion parametersâ i , as in the following: where S c (t m , f r ) is the range-compressed signal in t m − f r domain, FT t m (·) stands for the Fourier transform operation with respect to t m and IFT f r (·) denotes the inverse Fourier transform operation with respect to f r . It is worth pointing out that in the case of extended target theâ i , i = 0, 1, . . . , Q, used in (25) should be the motion parameters of the most representative point in the target, such as the centroid point [3].

G. DETAILED PROCEDURE
Based on the above derivation, the detailed flowchart of the proposed method is summarized in FIGURE 3. The specific procedures can be stated as follows: 1) Take the range compression to obtain the complex range profiles s c (t m , t n ); 2) Eliminate part of noise via the primary thresholding preprocessing; 3) Perform Hough transform on the threshold-crossing envelope range profiles s c (t m , t n );

4)
Apply the secondary thresholding processing and the masking operation to further suppress the clutter, and then estimate the target's average radial velocity based on the extracted medoid point (ρ c , θ c ); 5) Correct the first-order range migration with the estimated average radial velocityâ 1 ; 6) Extract the signals at the medoid range sampling cell along the pulses (slow time) domain based on the probability of the location of the target; 7) Resolve the Doppler ambiguity withâ 1 ; 8) Perform PCT on s d (t m ) to generate the TFR; 9) Extract the time-frequency ridge feature from the TFR; 10) Estimate the residual errors ε i a , i = 1, 2, · · · , Q, via the polynomial regression; 11) Estimate and refine the remaining motion parameters; 12) Re-correct the range migration and estimate the initial radial range, and output the results,â i , i = 0, 1, . . . , Q; 13) DFMC and energy Integration, for motion parameter estimation this step is optional, which is not included in FIGURE 3.

H. COMPUTATIONAL COST
In what follows, we will analyze the computational cost (complexity) of the major steps of HPCT. As mentioned in the previous section, the numbers of range sampling cells and received pulses are N and M , respectively. In the Hough transform, the ρ − θ space is quantized in the ρ and θ dimensions with each cell of size ρ × θ, the corresponding numbers of grids are N ρ and N θ respectively. In the PCT, the Gaussian window length is N win , the number of epochs is N iter and the polynomial order is Q.
First, since the computational cost of a 1-D FFT of N points is O N log 2 (N ) [37], the cost (by separability) for the range compression is O 2MN log 2 (N ) + MN . Second, as to detecting the target moving trajectory and estimating the target's average velocity via HT, the computational cost in this step mainly depends on the number of the primary threshold crossing data points (N e ) and the length N l of the lines formed in the t n − t m space, the computational cost is in the order of O (N e N l ) [37]. Third, for the first-order RM compensation and the Doppler ambiguity resolution, the cost is O MN + MN log 2 N + M . Fourth, for the PCT operation, the computational cost is O N iter · NN win log 2 N win + NN win . Fifth, for the time-frequency ridge feature extraction and the polynomial regression, the cost is O N iter Q 3 + MQ 2 2 4 [38]. Sixth, as to re-correcting the range migration and estimating the remaining motion parameters, the computational cost is  On the other hand, the computational cost for GRFT is O MNN v N a N j [12], [13], [32] where M , N , N v , N a , N j denote the number of echo pulses, range sampling cells, searching velocity, searching acceleration and searching jerk, respectively. Suppose that N = N v = N a = N j = M for the sake of simplicity, then the computational burden of GRFT is O M 5 . We compare the computational costs of HPCT with that of GRFT by FIGURE 4. It is evident that the computational cost of GRFT is significantly larger than that of HPCT (the proposed method). For example, when M = 256, the computational cost of HPCT is about 3.9 × 10 −5 of the GRFT; and when M = 1024, it is about 6.5 × 10 −7 of the GRFT.

IV. NUMERICAL EXPERIMENTS AND PERFORMANCE ANALYSIS
In this section, several numerical experiments are presented to demonstrate the effectiveness of the proposed method for the motion parameters estimation of a low-observable, highspeed maneuvering target with jerk motion in the Gaussian white noise background, where the key parameters of radar and the target are listed in TABLE 3.
In order to reduce the side-lobe of the compressed signal, the hamming weighting is adopted, the corresponding signalto-noise ratio loss and main-lobe broadening factor are about 1.36 dB and 1.47 [42], respectively. It should also be noted that, in this paper, we do not adopt the Rayleigh range resolution of r = c 2B (i.e., the main-lobe −4 dB width), but employ the −10 dB width metric [42] to characterize the range resolution ( r ≈ 1.5 · c 2B). After windowing operation, the r will be further broadened to 1.5 · k w · c 2B with the broadening factor k w . Consequently, the hamming weighted Rayleigh range resolution is about 7.4 m, and the corresponding range resolution defined by the main-lobe −10 dB width is about 11.0 m; the gain of range compression is about 26 dB. In addition, the unambiguous range, unambiguous velocity and observation time are equal to 75 km, 100 m/s, and 0.512 s, respectively. It must be pointed out that the dynamic ranges of the following two-dimensional (2-D) images in this paper are all set to 30 dB for showing the image contents clearly.

A. ENERGY DISTRIBUTION OF RADAR RETURNS
In order to get a better feeling of the effects of the range and Doppler frequency migrations, we first analyze the energy distribution of the range-compressed radar returns in the range-pulses number domain with FIGURE 5, where the SNR of the received radar echo signals before pulse compression is set to 0 dB for showing the energy distribution clearly. Because of the target's relative motion to the radar, it can be observed that the slant range between radar and the target has changed from 30.00 km to 29.47 km within the observation time interval. In other words, the signal energy spreads across 48 range resolution cells, which is defined by the main-lobe −10 dB width as mentioned at the beginning of this section.
By taking the Fourier transform of the range-compressed signal s c (t m , t n ) along the slow-time t m domain, we obtain the energy distribution of the radar returns in the range-Doppler frequency domain, as shown in FIGURE 6. It can be clearly observed again that the energy of radar returns spreads across both the range and Doppler coordinates severely. Therefore,  particular attention must be paid to this issue. Otherwise, it is difficult to improve the radar performance even the dwell time of antenna beam can be prolonged, because only part of signal energy can be exploited.

B. MOVING TRAJECTORY DETECTION VIA HT
As mentioned in Section III.A, the Hough transform has some valuable advantages for detecting the target's trajectory. In this subsection, we qualitatively evaluate its ability to detect the moving trajectory of a high-speed maneuvering target in complex Gaussian white noise background by FIGURE 7, where the SNR before pulse compression is set to −23 dB. The result after range compression is shown in FIGURE 7 (a). For simplicity, the scaling factor λ A in (6) is set to 0.5 and the result of the primary threshold is shown in FIGURE 7 (c). As a point of reference, the instantaneous slant ranges measured from FIGURE 7 (a) and (c) via the envelope cross-correlation function (CCF) algorithm, which is a simple yet popular algorithm that is frequently used for range alignment [3], are illustrated in FIGURE 7 (b) and (d), respectively. In the figure, the ''GT'' label denotes the ground truth-value of instantaneous slant range. It can be observed from FIGURE 7 (b) and (d) that there are a large number of outliers in the measured data, even though a large part of noise has been omitted via the threshold transform. Intuitively, we can think of the reason as the correlation between range-profiles would be corrupted when the SNR is low to 35186 VOLUME 9, 2021 Furthermore, after compensating the first-order range migration (a.k.a., range walk) with the estimated average velocityâ 1 , we obtain FIGURE 7 (g) and its zoomed-in version FIGURE 7 (h). The estimated probability of the location of the target is illustrated in FIGURE 7 (i). It can be observed that the range walk has been removed and the target returns are concentrated at 30.003 km. However, the nonlinear range migration (range curvature) induced by the target's higher order motions still exists, which will be processed in the following subsection.

C. TFR AND TF-RIDGE FEATURE EXTRACTION
As described in Section III, We adopt the PCT to produce the time-frequency representation (TFR) of s d (t m ), and estimate the remaining motion parameters, including the residual velocity, acceleration, jerk, and so on, with the time-frequency ridge feature extracted from the TFR. The related results are presented in FIGURE 8. For comparison, the TFRs generated by the short-time Fourier transform (STFT), Wigner-Ville distribution (WVD) and PCT are shown in FIGURE 8 (a), (b) and (c), respectively. It can be seen that the STFT cannot achieve highly concentrated TFR, and the WVD suffers from the undesired cross-term interference. In contrast, the PCT can produce a cross-termfree TFR with an excellent concentration, which is favorable for extracting the time-frequency ridge feature. The estimated instantaneous frequency (IF) from the ridge of TFR is given in FIGURE 8 (d). The change of IF versus the number of epochs is given in FIGURE 8 (e), which is evaluated by As we can see in FIGURE 8 (e), the IF difference between the fourth epoch and the third epoch is almost equivalent to zero. In other words, in this case the estimated IF has reached convergence after the third epoch. Furthermore, with the refined motion parameters, the RM re-corrected range profiles and its zoomed-in version, and the re-estimated probability of the location of the target are shown in FIGURE 8 (f), (g) and (h), respectively. Compare FIGURE 8 (g) with FIGURE 7 (h), it can be seen that the residual nonlinear range migration (range curvature) induced by the target's acceleration and jerk has been removed. After that, the 1-D focused result S focused,1D (f d ) is presented in FIGURE 8 (i), from which we can see that the noise has been effectively suppressed and the energy of radar returns can be efficiently integrated by HPCT.

D. PARAMETERS ESTIMATION PERFORMANCE
The estimation errors of target's motion parameters will affect the performance of compensating the effects of range and Doppler frequency migrations, and impact the coherent integration ability to accumulate the signal energy along the target's moving trajectory. In this subsection, the motion parameters estimation performance of the proposed method is analyzed via the Monte Carlo experiment under different SNR levels. For simplicity, the variance of the additive complex Gaussian white noise is fixed to 1, and the input SNR is changed from −40 dB to +20 dB by adjusting the power of the radar echo signals. If the SNR is in the interval [-36 dB, −34 dB], then the corresponding step size is set to 0.2 dB and is set to 1 dB otherwise. For each SNR value, 200 times Monte Carlo trials are performed to calculate the root mean squared errors (RMSEs), which is defined as where N MC stands for the number of Monte Carlo trials,ŷ i and y i respectively denote the estimated value and the ground truth value of the motion parameter in the ith Monte Carlo trial. The measured RMSEs in decibels for the initial slant range, velocity, acceleration and jerk via the CCF and the HPCT are shown in FIGURE 9 and FIGURE 10, respectively. It can be  observed that the input SNR threshold region [2] of CCF for parameter estimation is nearly range from −23 dB to −18 dB. In contrast, the corresponding region of HPCT is nearly from −35.2 dB to −35 dB, which is much lower and narrower than that of CCF. In other words, the parameter estimation performance of HPCT is significantly better than that of CCF. The RMSEs at the right boundary points (RBPs) of the above two SNR threshold regions are listed in TABLE 4. It should be pointed out that we do not compare the RMSEs of HPCT with those of GRFT via Monte Carlo experiment. The main reason is that the computational cost of GRFT for the high-speed maneuvering target with complex motions is so heavy that evaluating the RMSEs of GRFT is a tedious job. In addition, we think it is meaningless to evaluate the RMSEs using the prior information heavily to reduce the time cost. As an alternative, for the sake of simplicity, we compare the RMSEs of HPCT with the searching intervals of GRFT, which may be regarded as the RMSE bounds of GRFT in low-SNR case. With the −10 dB main-lobe width metric mentioned in the previous section, the range resolution is about r ≈ k w · 1.5 · c 2B and the velocity resolutions is about v ≈ k w · 1.5 · λ (2T obs ), where k w ≈ 1.47 is the corresponding main-lobe broadening factor. According to [13], the searching intervals of initial slant range and velocity are set to δr = r (about 11.03 m) and δv = v (about 0.22 m/s) respectively, and the steps of acceleration and jerk are set to δa = v T obs ≈ 0.42 m/s 2 and δj = v T 2 obs ≈ 0.82 m/s 3 . Compare these results with the RMSEs of HPCT in TABLE 4 at the right boundary point of the input SNR threshold region, it can be seen that the HPCT can achieve close parameter estimation performance as the GRFT algorithm when the input SNR is no less than −35 dB.

E. INTEGRATION FOR A MANEUVERING TARGET
In this subsection, we evaluate the coherent integration performance of HPCT for a high-speed maneuvering target with the motion parameters described in TABLE 3, where the additive complex Gaussian white noise is added to the radar returns and the SNR before pulse compression is manually set to −35 dB. The reason why we set the input SNR to −35 dB is that the right boundary point of the SNR threshold region of the HPCT is at −35 dB. The 1-D and 2-D coherent integration results can be obtained via (24) and (25) with the estimated motion parameters, as shown in FIGURE 11 and FIGURE 12, respectively. For comparison, the corresponding coherent integration results of the moving target detection (MTD), envelope across correlation function (CCF), Radon Fourier transform (RFT), Hough transform (HT) and Hough-modified discrete Chirp Fourier transform (HMD-CFT) are also provided in FIGURE 11 and FIGURE 12. In the figures, the HMDCFT denotes the method we proposed to estimate the target's velocity and acceleration by combining the Hough transform with the MDCFT [43] algorithm.
In the case of MTD, the motion parameters used in (24) and (25) are all set to zeroes because it does not need them. Its 1-D and 2-D integration results are shown in FIGURE 11 (a) and FIGURE 12 (a), respectively. It can be observed that the MTD outputs are submerged in the noise. The main reason for this is that the MTD does not consider the RM and DFM effects when using the Doppler filter bank to achieve the coherent integration. The envelope CCF method can compensate the RM and DFM effects simultaneously in high SNR situation. However, as mentioned in TABLE 4, because the SNR threshold of CCF for motion parameter estimation is about −18 dB, the estimated motion parameter is worthless when the input SNR is equivalent to −35 dB. Hence, the results of envelope CCF are submerged in the noise too, as shown in FIGURE 11 (b) and FIGURE 12 (b). The RFT and HT methods can estimate the target's velocity. Therefore, the first-order RM and Doppler ambiguity can be corrected and the peaks of target are higher than the noise level. However, the range curvature and highorder DFM caused by the target's acceleration and jerk  cannot be compensated so that the output SNRs of the RFT and HT are unacceptable, as shown in FIGURE 11 (c), (d) and FIGURE 12 (c), (d). Because the HMDCFT method can compensate the RM and DFM effects caused by the target's velocity and acceleration, the peaks of FIGURE 11 (e) and FIGURE 12 (e) are sharp, However, the effects caused by the target's jerk do not be corrected, which will cause some performance loss. As mentioned in Section IV.D, the HPCT can estimate the target's motion parameters accurately when the input SNR is no less than −35 dB. Therefore, the effects of RM and DFM are accurately compensated and the energy of radar returns is well focused as a sharp peak in the outputs. Compare FIGURE 11 (f)/FIGURE 12 (f) with other subfigures, it can be observed that the HCPT can achieve the sharpest main-lobe and the lowest background of noise.

F. COHERENT INTEGRATION PERFORMANCE
In this subsection, the integration performances of MTD, CCF, RFT, HT, HMDCFT and the proposed method (HPCT) are further investigated by Monte Carlo trials under different SNR levels (in the range [−50, 0] dB with a step of 1 dB) with the parameters described in TABLE 3. For each SNR value, the mean peak amplitude (MPA) of S focused,1D (f d ) is first calculated by where MD denotes the methods that are used to estimate the motion parameters and integrate the energy of radar returns, including MTD, CCF, RFT, HT, HMDCFT, HPCT and GT. The ''GT'' denotes that the ground truth-values of the target's motion parameters are directly used in (24) to integrate the energy of radar returns, which is used as a reference for the integration performance evaluation. The number of Monte Carlo trials N MC in (29) is set to 500. After that, we normalize the MPAs by max [MPA GT (snr)] and convert them to decibels, i.e.
The corresponding curves generated by MTD, CCF, RFT, HT, HMDCFT, HPCT and GT respectively are shown in FIGURE 13. It can be seen that the proposed method (HPCT) can achieve close integration performance as the ground truth-values of the target's motion parameters are directly used when the input SNR is no less than −35 dB, which may indicate that the estimated motion parameters is sufficiently accurate to integrate the energy of radar returns along the target's moving trajectory. The curves also illustrate that the integration performance of the proposed method is superior to that of HMDCFT method. This is because the radial jerk of target is not taken into consideration in the HMDCFT method and its peak amplitude would be decreased by about 2.5 dB  with respect to that of HPCT. It is also worth pointing out that the CCF-based method can achieve an acceptable integration performance when the input SNR is greater than −18 dB. It may also hint that the accuracy of the estimated motion parameters via CCF can reach an acceptable level in this case, which is consistent with our analysis about parameters estimation performance in Section IV.D. In contrast, because the RFT and HT can only compensate the first-order effects of RM and DFM and the MTD does not consider these two effects, accordingly, these three methods suffer from serious integration performance loss as shown in FIGURE 13.
From a quantitative analysis point of view, when the input SNR is equal to −35 dB, the normalized mean peak amplitude (NMPA) of HPCT is about −17.5 dBV that is denoted by a gray horizontal line in  After estimating the target's motion parameters and integrating the energy of radar returns along the target's moving trajectory, the radar target detection is pretty straightforward by combining the Neyman-Pearson detector or the constant false alarm ratio (CFAR) detector [2] with the proposed method. Because the target detection problem is very similar to the above discussion except the detection threshold is generally determined by a predetermined false alarm rate, we will not discuss it here.

G. EXTENSION TO MULTIPLE TARGETS SCENE
In the previous sections, the proposed method has been analyzed in the single target scene. However, it can also be extended to the scenario with multiple targets. In this subsection, we will analyze the parameters estimation performance and the coherent integration performance for multiple targets of the proposed method by an example, as shown in FIGURE 14. The key parameters of radar are the same as those in TABLE 3, the motion parameters of three maneuvering targets are listed in TABLE 6 and the SNRs of these three targets before pulse compression are all set to −23 dB for simplicity and for showing the results clearly. The targets in the scene have different motion orders. Moreover, the radial velocity of target A is set equal to that of target B, and the initial slant range of target C is the same as that of target A. According to TABLE 3, because we employ  the low-PRF scheme to increase the range unambiguous observation scope, the radar will be extremely ambiguous in Doppler for these three targets.
The result after range compression is shown in FIGURE 14 (a), and FIGURE 14 (b) shows the energy distribution of the radar returns in the range-Doppler frequency domain. The targets' movements not only make the moving trajectories of these three targets intersect in the range-pulses number domain but also cause the energy of radar returns spread in both the time and Doppler frequency domains severely. After taking the Hough transform, the moving trajectories of these three targets can be accumulated as three peaks in ρ − θ space as shown in FIGURE 14 (c). With the estimated average velocities of target A, B and C, the results after the first order range migration compensation (FRMC) are shown in FIGURE 14 (d), (e) and (f), respectively. The estimated probabilities of the locations of the targets are illustrated in FIGURE 14 (g), (h) and (i), respectively. It can be seen that there are two sharp peaks in FIGURE 14 (g) and (h) because target A and B have the same radial velocity. With the help of the estimated initial slant range by HT, a cyan area is set up to resolve this type of ambiguity. In addition, the gray horizontal line is a threshold used for controlling the false alarm further, which is set equal to √ 2 times the mean value of p r (n) for simplicity. After extracting the signals from the most representative range sampling cells of target A, B and C and resolving the Doppler ambiguity, the TFRs are shown in FIGURE 14 (j), (k) and (l), respectively. The time-frequency trajectories can clearly reflect the kinematic features of targets. Finally, the 2-D coherent integration results along with the estimated motion parameters of target A, B and C are shown in FIGURE 14 (m), (n) and (o), respectively. It can be observed that the estimated motion parameters agree well with the ground truth-values listed in TABLE 6, and the signal energy of each target is coherently accumulated as a sharp peak in the corresponding output, which is helpful for the target detection.

V. CONCLUSION
For the low-observable high-speed maneuvering target with high-order motions, based on the Hough transform (HT) and the polynomial Chirplet transform (PCT) we proposed a new method, i.e., HPCT, to estimate the target's motion parameters in this paper. This method can correct the range walk, range curvature and Doppler frequency migration simultaneously, and estimate the target's motion parameters accurately, by using the trajectory feature extracted from the image of range profiles and the time-frequency ridge feature extracted from the highly concentrated TFR. After that, the energy of radar returns can be coherently integrated along the target's moving trajectory. Moreover, we have evaluated the performance of the proposed method via several numerical experiments. The results showed that the HPCT can achieve a much lower and narrower input SNR threshold region than the widely used envelope CCF algorithm, that is, it has a much better noise resistance capability. The results also highlighted that the HPCT can avoid the cross-term interferences, and estimate the target's motion parameters accurately with a much lower computational cost compared with the GRFT. In addition, from the coherent integration perspective, under the experimental condition described in TABLE 3, the proposed method is superior to MTD, CCF, RFT, HT and HMDCFT by about 30 dB, 15 dB, 21 dB, 21 dB, 5 dB, respectively, thanks to its ability to deal with the high-order RM and DFM simultaneously. Finally, we extended the proposed method to the scenario of multiple targets that have different motion orders, and demonstrated its effectiveness with an example.
The proposed method also has its apparent limitation: although it was extended to the multiple targets scene in Section IV. G, we only considered the situation in which the targets are separable in the standard Hough parameter space. For the case when the targets' motion parameters are very similar (with the same initial slant range, same velocity, but different on acceleration and jerk), it may be computationally very expensive to separate the targets in higher dimension using the generalized Hough transform. Fortunately, the proposed method can avoid the cross-term interferences. It may be possible to separate the very similar targets in both the standard Hough parameter space and the time-frequency domain to reduce the computational cost and improve the parameter estimation accuracy. A possible future work might concern the extension of the proposed method to a more complicated situation with multiple targets.
CHAO ZENG received the B.S. degree from Yunnan University, China, in 1990, the M.S. degree from the Graduate School, China Academy of Engineering Physics, Mianyang, China, in 1993, and the Ph.D. degree from the Beijing Institute of Technology, in 2004. He is currently a Professor of Electronic Engineering with the Institute of Electronic Engineering, China Academy of Engineering Physics. His research interests include radar signal processing and system reliability.
HAI ZHANG received the B.S. degree from the Hefei University of Technology, Hefei, China, in 1990, and the M.S. and Ph.D. degrees from the University of Electronic Science and Technology of China, Chengdu, China, in 1994 and 2001, respectively. He is currently a Professor with the Institute of Electronic Engineering, China Academy of Engineering Physics. His research interests include radar signal processing, radar image processing, and the design of radar altimeter.
GE JIANG received the B.S. and M.S. degrees from the University of Electronic Science and Technology of China, Chengdu, China, in 2004 and 2007, respectively, and the Ph.D. degree in radio physics from the Graduate School, China Academy of Engineering Physics, Beijing, China, in 2017. He is currently an Associated Professor with the Institute of Electronic Engineering, China Academy of Engineering Physics. His research interests include signal processing and radar system design. VOLUME 9, 2021