Instantaneous Frequency Estimation Based on Modiﬁed Kalman Filter for Cone-Shaped Target

: The instantaneous frequency (IF) is a vital parameter for the analysis of non-stationary multicomponent signals, and plays an important role in space cone-shaped target recognition. For a cone-shaped target, IF estimation is not a trivial issue due to the proximity of the energy of the IF components, the intersections among di ﬀ erent IF components, and the existence of noise. Compared with the general parameterized time-frequency (GPTF), the traditional Kalman ﬁlter can perform better when the energy of di ﬀ erent signal components is close. Nevertheless, the traditional Kalman ﬁlter usually makes association mistakes at the intersections of IF components and is sensitive to the noise. In this paper, a novel IF estimation method based on modiﬁed Kalman ﬁlter (MKF) is proposed, in which the MKF is used to associate the intersecting IF trajectories obtained by the synchroextracting transform (SET). The core of MKF is the introduction of trajectory correction strategy in which a trajectory survival rate is deﬁned to judge the occurrence of association mistakes. When the trajectory survival rate is below the predetermined threshold, it means that an association mistakes occurs, and then the new trajectories generated by the random sample consensus algorithm are used to correct the wrong associations timely. The trajectory correction strategy can e ﬀ ectively obviate the association mistakes caused by the intersections of IF components and the noise. The windowing technique is also used in the trajectory correction strategy to improve computational speed. The experimental results based on the electromagnetic computation data show that the proposed method is more robust and precise than the traditional Kalman ﬁlter. Moreover, the proposed method has great performance advantages compared with other methods (i.e., the multiridge detection, the ant colony optimization, and the GPTF methods) especially in the case of low signal noise ratio (SNR).


Introduction
Relative motion of nonrigid parts of complex targets can induce micro-Doppler effect [1]. Micro-Doppler effect widely exists in the radar targets, for example, the rotating helicopter blades, precession cone-shaped targets, and the walking people. Moreover, the micro-Doppler effect appears as sidebands around the central Doppler frequency. Since different kinds of targets usually induce different types of micro-Doppler modulations, which can be regarded as a unique feature, there has been a lot of research applying the micro-Doppler effect to the accurate target identification applications [2][3][4][5].
Because of the time-varying properties of micro-Doppler, time-frequency (TF) analysis of non-stationary radar signals can provide information about time-dependent Doppler frequency variations and allow to extract classification features. In the last decade, several radar automatic target recognition techniques and activity classification algorithms, which are based on TF analysis (TFA) have been presented. After analyzing the differences between the micro-Doppler signatures from the three kinds of ground moving targets, a three-dimensional TF feature vector is extracted in [6] for categorizing ground moving targets. The paper [7] proposes a novel algorithm for the estimation of blade length and rotations per second of drone using TFA. In [8], the estimation method of spin rate and nutation angle based on the TF spectrum for free rigid targets is proposed.
For the multicomponent signal, instantaneous frequency (IF) of each component will appear as an IF trajectory which is determined by the micro-dynamic characteristics in the TF spectrum. Estimation of their IF characteristics for multicomponent signals is a widely investigated problem. The most IF estimation methods contain two critical steps: the first step is to obtain the concentrated TF spectrum and the second step is about tracking and associating each IF trajectory, i.e., the peak points of individual signal component in the TF spectrum. Viterbi algorithm (VA), as a representative method introduced in [9,10], has been applied in the IF estimation by utilizing the idea of dynamic programming. Another IF estimation approach based on an adaptive short-time Fourier transform (STFT) has been proposed in [11]. However, restricted by the Heisenberg uncertainty principle or undesired cross-terms, the IF estimation methods based on the classical TF analysis algorithms, such as Wigner-Ville distribution and STFT have poor TF resolution and cannot estimate the IF precisely. In order to achieve the performance of ideal TFA (ITFA) [12], Daubechies and Maes propose synchrosqueezing transform (SST) [13], which sharpens the time-scale representation given by continuous wavelet transform to improve the TF resolution. Subsequently, an extension of SST to the short-time Fourier transform is proposed in [14]. Compared with SST, which squeezes all TF coefficients into the IF trajectory, the synchroextracting transform (SET) [15] only retains the most valued TF information of STFT and has better TF energy concentration.
For a cone-shaped target with micro-motion, there are multiple IF components. Two scattering points are usually visible and the other one only can be observed at particular sight angles due to the occlusion. Each IF component is a standard sine form or in the form of the multistage superimposed sine series, and these IF components are always intersecting. Since there are many peak points at each time in the TF spectrogram, it is not easy to identify which peak points correspond to which components, especially in the presence of the intersections and noise. Thus, how to track and associate the IF trajectory for a cone-shaped target is also a crucial issue. A popular multiridge detection (MD) algorithm is proposed in [16,17] by calculating the minimum value of the energy functional to obtain the estimation of the IF. It is worth noting that the MD algorithm only considers the magnitudes and the absolute frequency variations, and does not account for the variation directions of the IF trajectory. Thus, the MD algorithm is hard to deal with the intersecting IF component. An IF estimation approach based on the ant colony optimization (ACO) is proposed in [18] by taking into account that edges represent image segments with high energy to overcome the high noise impact. In practice, the local convergence problem reduces the stability of the algorithm and the algorithm will be invalid when IF of each component is intersecting. An IF estimation algorithm based on Kalman filter is proposed in [19]. The Kalman filter is a time-saving association algorithm, which can predict the locations of missing points and filter out the spurious points caused by noise. Nevertheless, the wrong association may occur because of the intersections of IF components. Afterwards, a method named the coherent single range Doppler interferometry-modified general parameterized TF (CSRDI-MGPTF) [20] solves the problem of intersection with a multicomponent signal-separating operator, which reduces the impact of each component. However, when the energy of different signal components is close, this method is inapplicable because the TF ridge, which is obtained by modified GPTF transform, is the aliasing of multiple signal components. Besides, the residual component energy by the multicomponent signal-separating operator can cause serious estimation mistakes involving high noise.
In this paper, a novel IF estimation method based on the modified Kalman filter (MKF) is proposed. The initial intersecting IF trajectories are obtained by the SET and the MKF is proposed to associate the intersecting IF trajectories. The traditional Kalman filter is an association algorithm, which can predict the locations of missing points and filter out the spurious points caused by noise. However, Remote Sens. 2020, 12, 2766 3 of 22 potential association mistakes usually cannot be avoided at the intersections of IF trajectories especially in the case of low signal to noise ratio (SNR). The core of the proposed MKF is the introduction of trajectory correction strategy, which can correct the incorrect associations. The proposed trajectory correction strategy defines a trajectory survival rate to judge the occurrence of incorrect associations drawing on the idea of random sample consensus (RANSAC) [21]. The predicting locations of points will gradually deviate from the real IF trajectories when an association error occurs and there is a decline in the trajectory survival rate. The RANSAC algorithm then is used to generate a new trajectory to correct the wrong association obtained by the Kalman filter when the trajectory survival rate is below the predetermined threshold. The proposed trajectory correction strategy can effectively avoid the association mistakes caused by the intersections of IF components and the noise. In addition, considering that the RANSAC algorithm has a heavy computational burden, the windowing technique is also used in the trajectory correction strategy to improve computational speed of the RANSAC algorithm. Besides the intuitive IF estimation results in the experiment, the quantitative evaluation results are also attained in terms of the relative root mean squared error (RRMSE), the mean absolute error (MAE), and the determination coefficient. Experimental results based on the electromagnetic computation data show that the proposed method has stronger robustness and higher accuracy than the traditional Kalman filter. Furthermore, the proposed method yields great performance advantages compared with other methods (i.e., the MD, the ACO, and the general parameterized TF (GPTF) methods), especially in the case of low SNR.
The plan of the work is as follows. After deriving the precession model of the cone-shaped target and SET in Section 2, this paper introduces MKF algorithm in Section 3. Numerical examples with the estimation results are presented to verify the performance of the proposed method in Section 4. Finally, conclusions are described in Section 5.

Procession Model
As shown in Figure 1, the reference coordinate system (O − x y z ) and target coordinate system (O − xyz) are established and the precession center of target is located at the origin O. When the target does precession, the precession angle is θ, ω s and ω c denote the angular velocity of the space target rotating around the axis Oz and OC respectively. Because of the symmetry of the target, it is generally assumed that the radar line-of-sight (RLOS) lies in the yOz plane. γ is the sight angle which is the angle between the RLOS and OC. β is the angle between the RLOS and axis Oz and this angle is time-varying. The smooth cone-shaped target contains three effective scattering points P 1 , P 2 , P 3 . Usually, the scatter P 3 is occlusive. According to [19], the distance between two scatters and radar can be obtained by where R 0 is the distance between the radar and the origin O, H, and r is the height and radius of the target respectively, h is the distance between the bottom surface and the origin O. The variation of β can be expressed as cos β(t) = cos γ cos θ − sin γ sin θ cos(w c t) The micro-Doppler formula of scatters can be described as where a = cos θ cos γ, b = sin θ sin γ. It is observed that the micro-Doppler formula of P 1 is a standard sine form and the micro-Doppler formula of P 2 is relatively complex which can be regarded as the superposition of the sine series of infinite order.
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 22 where cos cos a θ γ = , It is observed that the micro-Doppler formula of 1 P is a standard sine form and the micro-Doppler formula of 2 P is relatively complex which can be regarded as the superposition of the sine series of infinite order.

Synchroextracting Transform
The STFT of a given signal where ( ) g t is a real-valued window function with The expression (5) can be rewritten as [15] ( )( where ˆ( ) g ω is the Fourier transform of ( ) g t .
The two-dimensional IF estimation of the STFT result (6) can be obtained by In the framework of SET, the TF coefficient in the IF trajectory = ( , ) t ω ω ω is extracted to obtain a more concentrated TF representation: where the synchroextracting operator (SEO) Figure 2 shows the electromagnetic simulation results induced by the precession for the coneshaped target. The TF representations of the signal acquired by STFT and SET are shown in Figure  2a,b respectively. It can be seen that the result of SET has an obviously higher energy concentration than that of STFT. The more energy-concentrated TF result denotes the better ability of the TF location and is more favorable to IF estimation. Meanwhile, different signal components in the TF representation are intersecting, which brings a huge challenge for IF estimation. Thus, based on the TF representation obtained by SET, this paper proposes a MKF algorithm, which will be introduced in detail in the next section.

Synchroextracting Transform
The STFT of a given signal where g(t) is a real-valued window function with g(t) = g * (t).
The two-dimensional IF estimation of the STFT result (6) can be obtained bŷ In the framework of SET, the TF coefficient in the IF trajectory ω =ω(t, ω) is extracted to obtain a more concentrated TF representation: where the synchroextracting operator (SEO) δ(ω −ω(t, ω)) can be interpreted as Figure 2 shows the electromagnetic simulation results induced by the precession for the cone-shaped target. The TF representations of the signal acquired by STFT and SET are shown in Figure 2a,b respectively. It can be seen that the result of SET has an obviously higher energy concentration than that of STFT. The more energy-concentrated TF result denotes the better ability of the TF location and is more favorable to IF estimation. Meanwhile, different signal components in the TF representation are intersecting, which brings a huge challenge for IF estimation. Thus, based on the TF representation obtained by SET, this paper proposes a MKF algorithm, which will be introduced in detail in the next section.

Instantaneous Frequency Estimation
The most IF estimation methods contain two critical steps: the first step is to obtain the concentrated TF spectrum and the second step is about finding the peak point sequences that correspond to individual signal component in the TF spectrum. The second step can be seen as tracking and associating the IF trajectory. In a real scenario, there are often many peak points in the TF spectrum at each time, and their number often varies. In such circumstances it can be not easy to judge which peak points correspond to which components, and which are spurious points caused by noise. The traditional Kalman filter is an association algorithm, which can predict the locations of missing points and filter out the spurious points caused by noise. However, potential association mistakes usually cannot be avoided at the intersections of IF trajectories especially in the case of low SNR. Therefore, in this paper, the traditional IF estimation problem is transformed into the trajectory association problem and a novel IF estimation method is proposed based on the MKF.
To describe the trajectory association problem, the discrete system with the linear Gaussian dynamic and measurement model is given by The echo signal in a short-time window can be approximately a chirp signal and the frequency varies linearly with time at the moment [22]. In this paper, the constant velocity (CV) model is adopted. The state vector is The velocity is approximately constant within a short time. The state transition matrix and the measurement matrix are assumed to be invariant with time. The covariance of measurement noise and the covariance of process noise are assumed to be The flowchart of MKF is shown in Figure 3. The detailed MKF algorithm can be found in Algorithm 1.

Instantaneous Frequency Estimation
The most IF estimation methods contain two critical steps: the first step is to obtain the concentrated TF spectrum and the second step is about finding the peak point sequences that correspond to individual signal component in the TF spectrum. The second step can be seen as tracking and associating the IF trajectory. In a real scenario, there are often many peak points in the TF spectrum at each time, and their number often varies. In such circumstances it can be not easy to judge which peak points correspond to which components, and which are spurious points caused by noise. The traditional Kalman filter is an association algorithm, which can predict the locations of missing points and filter out the spurious points caused by noise. However, potential association mistakes usually cannot be avoided at the intersections of IF trajectories especially in the case of low SNR. Therefore, in this paper, the traditional IF estimation problem is transformed into the trajectory association problem and a novel IF estimation method is proposed based on the MKF.
To describe the trajectory association problem, the discrete system with the linear Gaussian dynamic and measurement model is given by where X i k denotes state vector of ith trajectory at time step k, Y i k denotes the measurement of ith trajectory at time step k. A represents the state transition and H represents the measurement matrix. ω i k−1 and e i k are mutually uncorrelated process noise with covariance R ω and measurement noise with covariance R e .
The echo signal in a short-time window can be approximately a chirp signal and the frequency varies linearly with time at the moment [22]. In this paper, the constant velocity (CV) model is adopted.
The flowchart of MKF is shown in Figure 3. The detailed MKF algorithm can be found in Algorithm 1.    The key procedures of algorithm are explained in detail as follows: 1.
One-step Kalman predictor is used to obtain the predictionX i T of the state vector X i k at k time step for each trajectory. Then update the state vector using the Kalman filter. The updated state vector X i k is estimated as follows: Remote Sens. 2020, 12, 2766 The Kalman prediction gain G i k is obtained by: The estimation error covariance matrix P i k is: and the predictionP i 2.
For filtering out the spurious points which are generated by noise and avoiding the association mistakes at the intersection, this paper introduces the inlier set S which is define as where y k j is the location of jth trajectory at the time step k. τ f and τ v are the threshold based on the actual situation. Because of the limitation of TF resolution, the points in the same frequency cell belong to one trajectory. Therefore, τ f is equal to the frequency resolution and τ v = τ f / t. The setting of velocity threshold is to avoid the association error at the intersection where the states of two trajectories are close but the velocities of the trajectories have an obvious difference. There may be more than one points in the inlier set due to existence of the spurious points caused by noise. If S is not empty, the mean value of the inlier set is regarded as the location of estimated trajectory. If S is empty, Hx i k is regarded as the location of estimated trajectory.

3.
Trajectory correction strategy is performed for each time step to correct the association mistakes caused by the intersections and noise. The RANSAC algorithm is utilized to generate trajectories when the trajectory survival rate is lower than the threshold ε 0 . The basic idea of RANSAC is to form simple data association hypotheses from a batch of data and verify it to the data.
As an iterative algorithm, RANSAC contains two sections: hypothesis generation and hypothesis evaluation [21]. In the hypothesis generation phase, RANSAC chooses a subset of data at random, and then estimates the parameters from the samples. Many assumptions are generated during the iteration. In the hypothesis evaluation phase, the most likely hypothesis is selected according to the maximum inlier candidates. A subset of data is considered as the inlier candidate, whose error is assumed to be within a predefined threshold. RANSAC is popular until now because it is easy to implement. The signal can be approximately regarded as the linear chirp signal in a short time, i.e., the trajectories in the short time window can be regarded as a simple straight line approximately. The windowing technique is used in the trajectory correction strategy to improve computational speed of the RANSAC algorithm. If the wrong association occurs because of the intersections and noise, subsequently, the false predicted state vectorX i k is obtained. In this case, the inlier point set obtained by the false state vectorX i k does not contain the points correspond to the true IF, and the trajectory survival rate will drop rapidly.

Performance Evaluation with Numerical Results
In this section, this paper focuses on the comparisons between the proposed method and other advanced IF estimation methods for the cone-shaped target. The experimental electromagnetic computation data is obtained by FEldberechnung bei Körpern mit beliebiger Oberfläche (FEKO) using the physical optics (PO) method. The performance of the proposed method with the following TF-based IF estimation methods is compared: • multiple target tracking (MTT) algorithm [19] which is based on the STFT; • CSRDI-MGPTF algorithm [20] which is based on the MGPTF; • ACO algorithm [18], which is based on the STFT; • multiridge detection (MD) algorithm [17], which is based on the STFT and SET.
For multicomponent signals, which are considered in our experiment, the cross-terms of Wigner distribution (WD) can confuse the auto-terms seriously and cause failures of ACO algorithm. Therefore, in this paper, the ACO algorithm is performed on the basis of STFT.
Quantitative evaluation criteria used in the experiments include the relative root mean squared error (RRMSE), the mean absolute error (MAE) and the determination coefficient (R 2 ) between the true IF and estimated IF, which are formulated as follows: where f n denotes the true IF andf n denotes the estimation IF, N is the length of signal. The determination coefficient R 2 is a statistical measure of how well the predictions approximate the real data. R 2 = 1 indicates that the predictions perfectly fit the real data. The negative value of R 2 indicates that the mean of the data provides a better prediction result than the original prediction. The Monte Carlo simulations are performed to calculate the RRMSE, MAE, and R 2 . Meanwhile, the performance is evaluated versus different SNR. Some simulation and target parameters are listed in Table 1. The measurement window length is set to be l = 10 and the threshold of trajectory survival rate is set to be ε 0 = 0.6. The parameters of measurement noise and process noise are set as r 1 = r 2 = 0.2 and q = 0.1. In this experiment, sight angle, precession rate and dwell time are set as γ = 45 • , ω c = π rad/s, T = 0.7 s respectively. At the moment, the scatter P 3 is invisible. The signal contains two well-separated and disjointed IF components. 5 show the TF spectrums and IF estimation results for different methods with SNR = 20 dB and SNR = 5 dB respectively. It can be seen that the two worst estimation results are exhibited by the ACO algorithm and MTT algorithm. For the MTT algorithm, the use of constant false alarm rate (CFAR) technology and plots centroid technology may cause the deviation of estimated IF especially in low SNR. Although the ACO algorithm is robust to high noise influence, the resolution of STFT is the performance-limiting factor. exhibited by the ACO algorithm and MTT algorithm. For the MTT algorithm, the use of constant false alarm rate (CFAR) technology and plots centroid technology may cause the deviation of estimated IF especially in low SNR. Although the ACO algorithm is robust to high noise influence, the resolution of STFT is the performance-limiting factor.  In order to make quantitative analysis, the RRMSEs, MAEs, and R 2 with different SNRs are given in Table 2. For the CSRDI-MGPTF algorithm, the first IF component can be modeled as sinusoid, which lead to the RRMSEs, MAEs, and R 2 of the first IF component are minimal for all SNRs. However, the second IF component is in the form of the complex multistage superimposed sine series. The short dwell time affects the estimation accuracy of model parameter especially in low SNR. On the contrary, in low SNR situation, the proposed MKF method has a better performance on the second IF component. Compared with the MD algorithm, the MKF algorithm can filter out the spurious points caused by noise and has more precise estimation results. From the Mean evaluation criteria, comprehensive performance analysis on total IF components, in Table 2, it can be obviously found that the proposed method outperforms all the other the IF estimation methods when the SNR is not higher than 15 dB and the performance advantages of the proposed algorithm are more obvious in the case of low SNR. Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 22 In order to make quantitative analysis, the RRMSEs, MAEs, and 2 R with different SNRs are given in Table 2 Table 2, it can be obviously found that the proposed method outperforms all the other the IF estimation methods when the SNR is not higher than 15dB and the performance advantages of the proposed algorithm are more obvious in the case of low SNR.  Example 2. In this experiment, sight angle, precession rate and dwell time are set asγ = 60 • , ω c = 1.8π rad/s, T = 0.7 s respectively. At the moment, the scatter P 3 is invisible. The signal contains two intersecting IF components and the energy of the two components is close.

Figures 4 and
The TF representations and IF estimation results of different algorithms are displayed in Figures 6 and 7. Obviously, the proposed MKF algorithm can effectively estimate two IF components, especially at the intersections of IF trajectories. The CSRDI-GPTF algorithm is invalid because the TF ridge obtained by MGPTF transform is the aliasing of two signal components at the moment. In Figures 6d and 7d, association errors occur at the intersections for the MTT algorithm because of the lack of the correction mechanism. The MD algorithms tries to look for an IF trajectory which goes through positions with large energy. However, MD algorithms do not consider the variation directions of IF trajectory, which may lead to the failure for estimating the IFs of intersecting components. The MD-STFT algorithm is easy to converge to an incorrect local extremum point and leads to wrong estimation of the first IF component when the energy of overlapped two components is close, as shown in Figures 6f and 7f. From Figures 6g and 7g, it can be seen that the MD-SET algorithm has significant estimation errors around the intersections. For ACO algorithm, the estimation cannot be performed because it is impossible to identify the regions of components and estimate IF of each region separately for the multicomponent signal containing overlapping components. Thus, in the experiment of example 2, example 3, and example 4, ACO algorithm is not adopted for comparison. of the lack of the correction mechanism. The MD algorithms tries to look for an IF trajectory which goes through positions with large energy. However, MD algorithms do not consider the variation directions of IF trajectory, which may lead to the failure for estimating the IFs of intersecting components. The MD-STFT algorithm is easy to converge to an incorrect local extremum point and leads to wrong estimation of the first IF component when the energy of overlapped two components is close, as shown in Figure 6f and Figure 7f. From Figure 6g and Figure 7g, it can be seen that the MD-SET algorithm has significant estimation errors around the intersections. For ACO algorithm, the estimation cannot be performed because it is impossible to identify the regions of components and estimate IF of each region separately for the multicomponent signal containing overlapping components. Thus, in the experiment of example 2, example 3, and example 4, ACO algorithm is not adopted for comparison.  The quantitative evaluations, RRMSEs, MAEs and R 2 , are provided in Table 3. Note that desirable performance of the MKF algorithm is achieved even for SNR = 5dB. From the column 5 and column 6 it can be observed the MTT algorithm and CSRDI-GPTF algorithm cannot provide accurate estimation results for both the first IF component and second IF component. Even more, the negative values of R 2 on the second IF component indicate that the estimation results are worse than the average of the original data. It can be concluded that the proposed IF estimation method has a superior mean estimation performance compared with other methods due to the trajectory correction strategy. Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 22     Example 3. In this experiment, sight angle, precession rate and dwell time are set as γ = 40 • , ω c = 1.8π rad/s, T = 0.7 s respectively. At the moment, the scatter P 3 is invisible. The signal contains two intersecting IF components and the energy of one of them is much bigger than that of the other one.
In Figures 8 and 9, (a) and (b) give the TF representations obtained by the STFT and SET, (c)~(g) show the IF estimation results of the proposed method and other compared methods. It is noted that the CSRDI-MGPTF algorithm has the best estimation result at SNR = 20 dB. Nevertheless, in the case of SNR = 5 dB, the residual component energy by the multicomponent signal separating method in the CSRDI-MGPTF algorithm may be larger than the other component energy, and this can cause serious estimation mistake for the other component. As shown in Figure 9e, the second IF is not accurately estimated. By contrast, the MKF algorithm still works well in the case of SNR = 5 dB, which denotes that the proposed algorithm has a good noise robustness. For the MTT, MD-STFT, and MD-SET algorithms, the IF estimation results at SNR = 20 dB, shown in Figure 8b,d,f, have an obvious estimation error near the intersection. Moreover, when SNR = 5 dB, the MTT and MD-STFT algorithms have worse estimation results for the second IF component. estimated. By contrast, the MKF algorithm still works well in the case of SNR=5dB, which denotes that the proposed algorithm has a good noise robustness. For the MTT, MD-STFT, and MD-SET algorithms, the IF estimation results at SNR=20dB , shown in Figure 8b,d,f, have an obvious estimation error near the intersection. Moreover, when SNR=5dB, the MTT and MD-STFT algorithms have worse estimation results for the second IF component.  Subsequently, Table 4 lists performance evaluation results. It is observed that the CSRDI-MGPTF algorithm has more precise estimation results than other algorithms when SNR is higher than 10 dB. It is worth noting that the MKF algorithm and the CSRDI-MGPTF algorithm give quite close results about R 2 when SNR is higher than 10 dB. However, in low SNR, i.e., SNR ≤ 10 dB, the proposed MKF has a better performance in terms of mean evaluation criteria. Moreover, in high SNR, the proposed MKF has an approximate performance with the CSRDI-MGPTF algorithm. Remote Sens. 2020, 12, x FOR PEER REVIEW 16 of 22 Subsequently, Table 4 lists performance evaluation results. It is observed that the CSRDI-MGPTF algorithm has more precise estimation results than other algorithms when SNR is higher than 10dB.
It is worth noting that the MKF algorithm and the CSRDI-MGPTF algorithm give quite close results about 2 R when SNR is higher than 10dB. However, in low SNR, i.e., SNR 10dB ≤ , the proposed MKF has a better performance in terms of mean evaluation criteria. Moreover, in high SNR, the proposed MKF has an approximate performance with the CSRDI-MGPTF algorithm.     The performance comparisons with the RRMSEs, MAEs, and 2 R are given in Table 5 for different values of SNR. It can be concluded from Figure 10 that the proposed IF estimation method shows a better accuracy compared with other methods in the more complex scenario. The performance comparisons with the RRMSEs, MAEs, and R 2 are given in Table 5 for different values of SNR. It can be concluded from Figure 10 that the proposed IF estimation method shows a better accuracy compared with other methods in the more complex scenario.  In the end, a brief analysis for the convergence of the estimation process is performed. The values of diagonal elements of covariance matrix P i k over time in example 2 are given in Figure 12. It can be seen that the proposed MKF algorithm has a good convergence.

Conclusions
In this paper, a novel method is proposed which combines the Kalman filter algorithm and the RANSAC algorithm to gain a more precise IF estimation with an acceptable computation cost. For the traditional Kalman filter, the advantage is its speed. However, potential association mistakes usually cannot be avoided at the intersections of IF trajectories especially in the case of low SNR. Moreover, for the traditional RANSAC algorithm, the batch trajectory processing is used to increase the accuracy and the noise robustness but brings the huge computation burden in the meanwhile.

Conclusions
In this paper, a novel method is proposed which combines the Kalman filter algorithm and the RANSAC algorithm to gain a more precise IF estimation with an acceptable computation cost.
For the traditional Kalman filter, the advantage is its speed. However, potential association mistakes usually cannot be avoided at the intersections of IF trajectories especially in the case of low SNR. Moreover, for the traditional RANSAC algorithm, the batch trajectory processing is used to increase the accuracy and the noise robustness but brings the huge computation burden in the meanwhile. How to guarantee precision with less computational cost in the proposed method is a difficulty we facing with. Therefore, we introduce the trajectory correction strategy to achieve the combination between the Kalman filter algorithm and the RANSAC algorithm. Experimental results based on electromagnetic computation data indicate that the proposed method can obtain good estimation performance and show the superiority in processing signal with high noise. It is worth mentioning that the MKF algorithm is a post-processing method. Therefore, the proposed method is applicable to other multicomponent signals and TF representations. IF of space cone-shaped target is determined by the micro-dynamic characteristics. Therefore, how to retrieve the physical structure and motion information of the target from estimated IF, and subsequent space cone-shaped target recognition need to be deeply investigated in future works. Besides, the IF estimation problem for the space targets of other shapes will also be investigated.

Conflicts of Interest:
The authors declare no conflict of interest.