Instantaneous Frequency Estimation Based on Modified Kalman Filter for Cone-Shaped Target

Abstract

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 different IF components, and the existence of noise. Compared with the general parameterized time-frequency (GPTF), the traditional Kalman filter can perform better when the energy of different signal components is close. Nevertheless, the traditional Kalman filter 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 modified Kalman filter (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 defined 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 effectively 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 filter. 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).

1. 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, 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.

2. Procession Model of the Cone-Shaped Target and SET

2.1. Procession Model

As shown in Figure 1, the reference coordinate system $\left(O-x^{\prime }{y}^{\prime }{z}^{\prime }\right)$ and target coordinate system $\left(O-xyz\right)$ are established and the precession center of target is located at the origin $O$. When the target does precession, the precession angle is $\theta$, ${\omega }_s$ and ${\omega }_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. $\gamma$ is the sight angle which is the angle between the $RLOS$ and $OC$. $\beta$ 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

$$\{\begin{array}{l}{r}_1(t)=R_0-(H-h)\cos \beta (t)\\ r_2(t)=R_0+hcos\beta (t)-r\sin \beta (t)\end{array}$$

(1)

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 $\beta$ can be expressed as

$$\cos \beta (t)=\cos \gamma \cos \theta -\sin \gamma \sin \theta \cos (w_{c}t)$$

(2)

The micro-Doppler formula of scatters can be described as

$$f_1\left(t\right)=2\omega b(h-H)\sin (w_{c}t)/\lambda$$

(3)

$$f_2\left(t\right)=-\frac{2}{\lambda }r\left(a+b\cos (w_{c}t)\right)\times \frac{w_{c}b\sin (w_{c}t)}{\sqrt{1-{\left(a+b\cos (w_{c}t)\right)}^2}}-\frac{2}{\lambda }{w}_{c}hb\sin (w_{c}t)$$

(4)

where $a=\cos \theta \cos \gamma$, $b=\sin \theta \sin \gamma$. 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.

2.2. Synchroextracting Transform

The STFT of a given signal $x(t)=A(t)e^{i\varphi (t)}$ is defined by

$$S_f^g(t,\omega )={\displaystyle {\int }_{-\infty }^{+\infty }x(\tau )g^{*}}(\tau -t)e^{-j\omega (\tau -t)}d\tau$$

(5)

where $g(t)$ is a real-valued window function with $g(t)=g^{*}(t)$.

The expression (5) can be rewritten as [15]

$$S_f^g(t,\omega )=A(t)e^{j\varphi (t)}\widehat{g}(\omega -{\varphi }^{\prime }(t))$$

(6)

where $\widehat{g}(\omega )$ is the Fourier transform of $g(t)$.

The two-dimensional IF estimation of the STFT result (6) can be obtained by

$$\widehat{\omega }(t,\omega )=\mathrm{Re}\left(-j\cdot \frac{{\partial }_t{S}_f^g(t,\omega )}{2\pi S_f^g(t,\omega )}\right)$$

(7)

In the framework of SET, the TF coefficient in the IF trajectory $\omega =\widehat{\omega }(t,\omega )$ is extracted to obtain a more concentrated TF representation:

$$Te(t,\omega )=S_f^g(t,\omega )\delta (\omega -\widehat{\omega }(t,\omega ))$$

(8)

where the synchroextracting operator (SEO) $\delta (\omega -\widehat{\omega }(t,\omega ))$ can be interpreted as

$$\delta (\omega -\widehat{\omega }(t,\omega ))=\{\begin{array}{cc}1& \omega =\widehat{\omega }(t,\omega )\\ 0& \omega \ne \widehat{\omega }(t,\omega )\end{array}$$

(9)

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.

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

$$\{\begin{array}{l}{X}_k^i=A{X}_{k-1}^i+{\omega }_{k-1}^i,{\omega }_{k-1}^i\sim N(0,R_{\omega })\\ Y_k^i=H{X}_k^i+e_k^i,e_k^i\sim N(0,R_e)\end{array}$$

(10)

where $X_k^i$ denotes state vector of $i\mathrm{th}$ trajectory at time step $k$, $Y_k^i$ denotes the measurement of $i\mathrm{th}$ trajectory at time step $k$. $A$ represents the state transition and $H$ represents the measurement matrix. ${\omega }_{k-1}^i$ and $e_k^i$ are mutually uncorrelated process noise with covariance $R_{\omega }$ 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 state vector is $X_k^i={[x_k^i,{\dot{x}}_k^i]}^T$, where $x_k^i$ denotes the position, ${\dot{x}}_k^i=\frac{x_k^i-x_{k-1}^i}{\mathsf{\Delta }t}$ denotes the velocity, $\mathsf{\Delta }t$ is the time lag. The velocity is approximately constant within a short time. The state transition matrix $A=\left[\begin{array}{cc}1& \mathsf{\Delta }t\\ 0& 1\end{array}\right]$ and the measurement matrix $H=\left[\begin{array}{cc}1& 0\\ 0& 0\end{array}\right]$ are assumed to be invariant with time. The covariance of measurement noise and the covariance of process noise are assumed to be $R_e=\mathrm{diag}\left(\right[r_1{r}_2\left]\right)=\left[\begin{array}{cc}{r}_1& 0\\ 0& r_2\end{array}\right]$ and $R_{\omega }=q\left[\begin{array}{cc}\frac{1}{3}△t^3& \frac{1}{2}△t^2\\ \frac{1}{2}△t^2& △t\end{array}\right]$.

The flowchart of MKF is shown in Figure 3. The detailed MKF algorithm can be found in Algorithm 1.

Algorithm 1. The MKF Algorithm.

Input: measurement $Y_{1:T}=S_f^g(t,\omega )$, window length $l$, threshold of trajectory survival rate ${\epsilon }_0$, number of signal component $I$.

for each time step $k=1:T$ do

update the trajectory survival rate ${\epsilon }_k^i=m_{k:k+l}^i/l$ where $m_{k:k+l}^i$ denotes the inlier

number of $i\mathrm{th}$ trajectories in a short time window.

if ${\epsilon }_k^i<{\epsilon }_0$

generate new trajectories.

end if

for each trajectory do

${\widehat{X}}_k^i=A{X}_{k-1}^i$.

compute inlier set $S$.

if $S=\varnothing$

$y_k^i=H{\widehat{X}}_k^i$.

else

$y_k^i=mean(S)$.

end if

update $X_k^i$ using $y_k^i$ by Kalman filter.

update the trajectories $T_i=[T_i,y_k^i]$.

end for

end for

Output: trajectories set $\left\{T_1,\dots ,T_I\right\}$

The key procedures of algorithm are explained in detail as follows:

1.

One-step Kalman predictor is used to obtain the prediction ${\widehat{X}}_k^i={[{\widehat{x}}_k^i,{\widehat{\dot{x}}}_k^i]}^T$ of the state vector $X_k^i$ at $k$ time step for each trajectory. Then update the state vector using the Kalman filter. The updated state vector $X_k^i$ is estimated as follows:

$$X_k^i={\widehat{X}}_k^i+G_k^i(y_k^i-H{\widehat{X}}_k^i)$$

(11)

The Kalman prediction gain $G_k^i$ is obtained by:

$$G_k^i={\widehat{P}}_k^i{H}^T{(H{\widehat{P}}_k^i{H}^T+R_e)}^{-1}$$

(12)

The estimation error covariance matrix $P_k^i$ is:

$$P_k^i={\widehat{P}}_k^i-G_k^i(H{\widehat{P}}_k^i{H}^T+R_e)G_k^i{}^T$$

(13)

and the prediction ${\widehat{P}}_k^i$ is

$${\widehat{P}}_k^i=A{P}_{k-1}^i{A}^T+R_{\omega }$$

(14)

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

$$S=\left\{y_{kj}|\left|y_{kj}-{\widehat{x}}_k^i\right|<\frac{{\tau }_f}{2},\left|\frac{y_{kj}-{\widehat{x}}_k^i}{△t}-{\widehat{\dot{x}}}_k^i\right|<{\tau }_v\right\}$$

(15)

where $y_{kj}$ is the location of $j\mathrm{th}$ trajectory at the time step $k$. ${\tau }_f$ and ${\tau }_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, ${\tau }_f$ is equal to the frequency resolution and ${\tau }_v={\tau }_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, $H{\widehat{x}}_k^i$ 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 ${\epsilon }_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 vector ${\widehat{X}}_k^i$ is obtained. In this case, the inlier point set obtained by the false state vector ${\widehat{X}}_k^i$ does not contain the points correspond to the true IF, and the trajectory survival rate will drop rapidly.

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

$$RRMSE=sqrt\left({\displaystyle \sum _{n=0}^N{\left(f_n-{\widehat{f}}_n\right)}^2}/\left(N{\displaystyle \sum _{n=0}^N{\left(f_n\right)}^2}\right)\right)$$

(16)

$$MAE={\displaystyle \sum _{n=0}^N\left|f_n-{\widehat{f}}_n\right|}/N$$

(17)

$$R^2=1-{\displaystyle \sum _{n=0}^N{\left(f_n-{\widehat{f}}_n\right)}^2}/{\displaystyle \sum _{n=0}^N{\left(f_n-{\displaystyle \sum _{n=0}^N{f}_n/N}\right)}^2}$$

(18)

where $f_n$ denotes the true IF and ${\widehat{f}}_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. Some Parameters.

Parameter	Value	Parameter	Value
Radar carrier frequency $f_c$	$10\mathrm{GHz}$	center height $h$	$1\mathrm{m}$
Pulse repletion frequency $prf$	$1000\mathrm{H}z$	Precession angle $\theta$	$5^{\circ }$
Target height $H$	$3\mathrm{m}$	Target bottom radius $r$	$0.5\mathrm{m}$

. The measurement window length is set to be $l=10$ and the threshold of trajectory survival rate is set to be ${\epsilon }_0=0.6$. The parameters of measurement noise and process noise are set as $r_1=r_2=0.2$ and $q=0.1$.

Example 1. 

In this experiment, sight angle, precession rate and dwell time are set as$\gamma ={45}^{\circ }$,${\omega }_c=\pi \mathrm{rad}/\mathrm{s}$,$\mathrm{T}=0.7\mathrm{s}$respectively. At the moment, the scatter$P_3$is invisible. The signal contains two well-separated and disjointed IF components.

Figure 4 and Figure 5 show the TF spectrums and IF estimation results for different methods with $\mathrm{SNR}=20\mathrm{dB}$ and $\mathrm{SNR}=5\mathrm{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.

In order to make quantitative analysis, the RRMSEs, MAEs, and $R^2$ with different SNRs are given in

Table 2. Performance evaluation of IF Estimation for example 1 based on the relative root mean squared error (RRMSE), the mean absolute error (MAE) and the determination coefficient ($R^2$ ).

Method			MKF	MTT	CSRDI- MGPTF	ACO	MD-STFT	MD-SET
RRMSE (×10−4)	5 dB	1st	8.7561	39.450	8.6603	44.102	9.2366	8.9419
		2nd	14.922	77.059	32.826	99.001	18.858	17.218
	10 dB	1st	8.4385	29.061	8.4151	41.743	8.8461	8.5507
		2nd	13.486	45.515	16.540	86.649	14.760	14.139
	15 dB	1st	8.3951	28.818	8.3115	41.652	8.7069	8.4868
		2nd	12.890	42.530	13.093	76.140	13.619	12.901
	20 dB	1st	8.3288	24.105	8.2771	39.214	8.4426	8.4245
		2nd	12.860	34.257	12.878	74.629	12.884	12.870
Mean RRMSE (×10−4)	5 dB		11.839	58.254	20.743	71.551	14.047	13.080
	10 dB		10.962	37.288	12.477	64.196	11.803	11.345
	15 dB		10.642	35.674	10.702	58.896	11.162	10.693
	20 dB		10.594	29.181	10.578	56.921	10.663	10.647
MAE	5 dB	1st	1.0866	4.8195	1.0275	7.3464	1.1284	1.0924
		2nd	1.1220	6.6894	2.3035	8.5916	1.6371	1.4947
	10 dB	1st	1.0498	3.1736	1.0157	6.1121	1.0830	1.0505
		2nd	1.1196	3.9473	1.5978	7.2762	1.6044	1.3369
	15 dB	1st	1.0025	3.1702	0.9869	6.0521	1.0613	1.0501
		2nd	1.1193	2.6299	1.1435	5.5773	1.4803	1.3269
	20 dB	1st	0.9938	2.7358	0.9702	5.9137	1.0242	0.9964
		2nd	1.1134	2.4103	1.0968	5.5770	1.1697	1.1347
Mean MAE	5 dB		1.1043	5.7544	1.6635	7.9690	1.3828	1.2936
	10 dB		1.0847	3.5605	1.3068	6.6942	1.3435	1.1937
	15 dB		1.0609	2.9001	1.0652	5.8147	1.2708	1.1885
	20 dB		1.0536	2.5731	1.0335	5.7454	1.0970	1.0641
$R^2$	5 dB	1st	0.9871	0.7377	0.9874	0.6722	0.9856	0.9865
		2nd	0.9651	0.0697	0.8312	−0.5355	0.9443	0.9536
	10 dB	1st	0.9880	0.8577	0.9881	0.7063	0.9868	0.9877
		2nd	0.9715	0.6755	0.9571	−0.1762	0.9659	0.9687
	15 dB	1st	0.9881	0.8600	0.9884	0.7076	0.9872	0.9879
		2nd	0.9740	0.7166	0.9731	0.0918	0.9709	0.9739
	20 dB	1st	0.9883	0.9021	0.9885	0.7408	0.9880	0.9880
		2nd	0.9741	0.8162	0.9740	0.1275	0.9740	0.9741
Mean $R^2$	5 dB		0.9761	0.4037	0.9093	0.0684	0.9650	0.9701
	10 dB		0.9798	0.7666	0.9726	0.2651	0.9764	0.9782
	15 dB		0.9811	0.7883	0.9808	0.3997	0.9791	0.9809
	20 dB		0.9812	0.8592	0.9813	0.4342	0.9810	0.9811

. 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 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$\gamma ={60}^{\circ }$,${\omega }_c=1.8\pi \mathrm{rad}/\mathrm{s}$,$\mathrm{T}=0.7\mathrm{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.

The TF representations and IF estimation results of different algorithms are displayed in Figure 6 and Figure 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 Figure 6d and Figure 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 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. Performance evaluation of IF Estimation for example 2.

Method			MKF	MTT	CSRDI- MGPTF	MD-STFT	MD-SET
RRMSE (×10−4)	5 dB	1st	10.870	259.63	280.72	262.86	24.192
		2nd	9.0626	433.09	517.35	9.4629	9.0738
	10 dB	1st	8.9880	255.43	321.23	262.01	22.260
		2nd	9.0490	425.01	448.46	9.0990	9.0599
	15 dB	1st	8.6752	252.02	313.33	236.00	17.929
		2nd	9.0410	420.05	468.54	9.0485	9.0507
	20 dB	1st	8.2833	251.09	333.50	233.47	12.818
		2nd	9.0338	418.61	418.23	9.0366	9.0366
Mean RRMSE (×10−4)	5 dB		9.9663	346.36	399.04	136.16	16.633
	10 dB		9.0185	340.22	384.85	135.55	15.660
	15 dB		8.8581	336.04	390.94	122.52	13.490
	20 dB		8.6585	334.85	375.86	121.25	10.927
MAE	5 dB	1st	1.8319	60.544	64.331	60.160	4.4248
		2nd	1.1534	69.018	51.118	1.1722	1.1559
	10 dB	1st	1.4906	60.351	75.898	60.904	4.2954
		2nd	1.1487	60.615	43.959	1.1687	1.1521
	15 dB	1st	1.4774	59.545	78.797	55.762	3.2361
		2nd	1.1236	59.909	49.650	1.1605	1.1508
	20 dB	1st	1.4630	59.326	72.361	55.549	3.0424
		2nd	1.1130	59.907	38.553	1.1585	1.1134
Mean MAE	5 dB		1.4927	64.781	57.725	30.666	2.7904
	10 dB		1.3197	60.483	59.929	30.631	2.7238
	15 dB		1.3005	59.727	64.224	28.461	2.1935
	20 dB		1.2880	59.617	55.457	28.354	2.0779
$R^2$	5 dB	1st	0.9992	0.5197	0.4385	0.5076	0.9958
		2nd	0.9994	−0.3386	−0.9162	0.9994	0.9994
	10 dB	1st	0.9994	0.5351	0.2647	0.5108	0.9965
		2nd	0.9994	−0.2891	−0.4353	0.9994	0.9994
	15 dB	1st	0.9995	0.5474	0.3005	0.6031	0.9977
		2nd	0.9994	−0.2592	−0.5667	0.9994	0.9994
	20 dB	1st	0.9995	0.5507	0.2074	0.6116	0.9988
		2nd	0.9994	−0.2506	−0.2483	0.9994	0.9994
Mean $R^2$	5 dB		0.9993	0.0906	−0.2389	0.7535	0.9976
	10 dB		0.9994	0.1230	−0.0853	0.7551	0.9980
	15 dB		0.9995	0.1441	−0.1331	0.8013	0.9986
	20 dB		0.9995	0.1501	−0.0205	0.8055	0.9991

. Note that desirable performance of the MKF algorithm is achieved even for $\mathrm{SNR}=5\mathrm{dB}$. 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.

Example 3. 

In this experiment, sight angle, precession rate and dwell time are set as$\gamma ={40}^{\circ }$,${\omega }_c=1.8\pi \mathrm{rad}/\mathrm{s}$,$\mathrm{T}=0.7\mathrm{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 Figure 8 and Figure 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 $\mathrm{SNR}=20\mathrm{dB}$. Nevertheless, in the case of $\mathrm{SNR}=5\mathrm{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 $\mathrm{SNR}=5\mathrm{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 $\mathrm{SNR}=20\mathrm{dB}$, shown in Figure 8b,d,f, have an obvious estimation error near the intersection. Moreover, when $\mathrm{SNR}=5\mathrm{dB}$, the MTT and MD-STFT algorithms have worse estimation results for the second IF component.

Subsequently,

Table 4. Performance evaluation of IF Estimation for example 3.

Method			MKF	MTT	CSRDI- MGPTF	MD-STFT	MD-SET
RRMSE (×10−4)	5 dB	1st	7.6616	51.084	4.2201	10.028	7.6658
		2nd	9.0178	699.85	217.78	359.20	13.202
	10 dB	1st	6.2131	10.467	3.9868	9.9708	6.2733
		2nd	8.1456	189.40	179.95	18.671	11.420
	15 dB	1st	6.0647	10.460	3.8160	8.9094	6.2253
		2nd	7.9485	47.192	7.1505	18.091	10.927
	20 dB	1st	6.0621	8.8946	3.8157	8.8951	6.0794
		2nd	7.9413	45.011	7.1457	12.068	10.912
Mean RRMSE (×10−4)	5 dB		8.3397	375.47	111.00	184.61	10.434
	10 dB		7.1793	99.934	89.968	14.366	8.8470
	15 dB		7.0066	28.826	5.4832	13.500	8.5762
	20 dB		7.0017	26.953	5.4807	10.482	8.4957
MAE	5 dB	1st	1.0548	3.7475	0.6992	1.8990	1.1437
		2nd	1.0475	89.413	33.090	50.117	1.5365
	10 dB	1st	0.9444	1.8816	0.6550	1.8882	1.0880
		2nd	0.8429	26.426	25.107	2.6050	1.3602
	15 dB	1st	0.7779	1.8808	0.5776	1.1179	0.9576
		2nd	0.6790	2.5844	0.4577	2.5241	1.0789
	20 dB	1st	0.7731	1.1195	0.5772	1.1179	0.7814
		2nd	0.6710	2.2578	0.4524	1.1437	1.0617
Mean MAE	5 dB		1.0512	46.580	16.896	26.008	1.3401
	10 dB		0.8937	14.154	12.881	2.2466	1.2241
	15 dB		0.7285	2.2326	0.5177	1.8210	1.0183
	20 dB		0.7221	1.6887	0.5148	1.1308	0.9216
$R^2$	5 dB	1st	0.9996	0.9814	0.9998	0.9993	0.9996
		2nd	0.9994	−2.4969	0.6614	0.0788	0.9688
	10 dB	1st	0.9997	0.9992	0.9998	0.9993	0.9997
		2nd	0.9995	0.7439	0.7688	0.9975	0.9991
	15 dB	1st	0.9997	0.9992	0.9998	0.9994	0.9997
		2nd	0.9995	0.9841	0.9995	0.9977	0.9991
	20 dB	1st	0.9997	0.9994	0.9998	0.9994	0.9997
		2nd	0.9995	0.9885	0.9996	0.9990	0.9992
Mean $R^2$	5 dB		0.9995	−0.7578	0.8306	0.5391	0.9842
	10 dB		0.9996	0.8716	0.8843	0.9984	0.9994
	15 dB		0.9996	0.9917	0.9997	0.9986	0.9994
	20 dB		0.9996	0.9940	0.9997	0.9992	0.9995

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\mathrm{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\mathrm{dB}$. However, in low SNR, i.e., $\mathrm{SNR}\le 10\mathrm{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.

Example 4. 

In this experiment, sight angle, precession rate and dwell time are set as$\gamma ={124}^{\circ }$,${\omega }_c=1.8\pi \mathrm{rad}/\mathrm{s}$,$\mathrm{T}=0.7\mathrm{s}$respectively. At the moment, the scatter$P_3$can be observed. The signal contains three intersecting IF components and the energy of the second IF component and third IF component is close.

A more complex condition is considered in the experiment. The intuitional results are shown in Figure 10 and Figure 11. Although the CSRDI-MGPTF has the best estimation performance of the first IF component, it is not suitable for the second IF component and the third IF component with close energy. However, in the meanwhile, the MKF algorithm has a better estimation result for the second IF component compared with the MTT, CSRDI-MGPTF, MD-STFT, and MD-SET algorithms. Besides, at $\mathrm{SNR}=20\mathrm{dB}$, the MKF still has a better IF estimation, whereas all of the IF estimation methods are unavailable at $\mathrm{SNR}=5\mathrm{dB}$ due to the third IF component is almost flooded by noise.

The performance comparisons with the RRMSEs, MAEs, and $R^2$ are given in

Table 5. Performance evaluation of IF Estimation for example 4.

Method			MKF	MTT	CSRDI- MGPTF	MD-STFT	MD-SET
RRMSE (×10−4)	5 dB	1st	10.683	13.817	3.8870	11.409	10.982
		2nd	32.993	653.12	296.50	378.49	288.59
		3rd	102.68	505.15	574.11	515.47	504.04
	10 dB	1st	10.677	13.568	3.7387	10.904	10.887
		2nd	19.985	644.98	299.67	371.65	25.870
		3rd	50.283	482.72	563.52	525.41	76.212
	15 dB	1st	10.423	13.436	3.6201	10.832	10.663
		2nd	10.861	26.514	367.45	367.59	20.102
		3rd	8.1118	16.169	585.27	426.64	16.100
	20 dB	1st	10.394	13.069	3.4951	10.611	10.656
		2nd	9.8399	21.054	129.24	654.14	18.491
		3rd	5.3313	8.9482	614.39	406.13	14.085
Mean RRMSE (×10−4)	5 dB		48.785	390.70	291.50	301.79	267.87
	10 dB		26.982	380.42	288.98	302.65	37.656
	15 dB		9.7986	18.706	318.78	268.35	15.622
	20 dB		8.5217	14.357	249.04	356.96	14.651
MAE	5 dB	1st	0.7446	1.0510	0.5179	0.8501	0.7541
		2nd	3.7771	79.599	36.814	35.131	37.033
		3rd	19.018	115.84	134.84	122.47	101.12
	10 dB	1st	0.7440	1.0267	0.4956	0.7662	0.7522
		2nd	2.3255	72.228	35.820	31.787	2.5958
		3rd	13.264	97.334	140.65	128.60	20.104
	15 dB	1st	0.7192	0.9776	0.4940	0.7592	0.7506
		2nd	1.5056	2.9666	51.089	31.089	2.5420
		3rd	2.1398	3.1612	144.39	94.386	2.7469
	20 dB	1st	0.7181	0.9497	0.4845	0.7505	0.7502
		2nd	1.0299	2.4010	15.143	86.170	2.5009
		3rd	1.0425	1.5168	137.27	86.608	2.4104
Mean MAE	5 dB		7.8466	65.497	57.391	52.817	46.302
	10 dB		5.4445	56.863	58.989	53.718	7.8173
	15 dB		1.4549	2.3685	65.324	42.078	2.0132
	20 dB		0.9302	1.6225	50.966	57.843	1.8872
$R^2$	5 dB	1st	0.9992	0.9986	0.9999	0.9991	0.9991
		2nd	0.9923	−2.0197	0.3777	−0.0141	0.4104
		3rd	0.9253	−0.8086	−1.3361	−0.8833	−0.8007
	10 dB	1st	0.9992	0.9987	0.9999	0.9992	0.9992
		2nd	0.9972	−1.9449	0.3643	0.0222	0.9953
		3rd	0.9821	−0.6516	−1.2508	−0.9566	0.9588
	15 dB	1st	0.9992	0.9987	0.9999	0.9992	0.9992
		2nd	0.9992	0.9950	0.0442	0.0435	0.9971
		3rd	0.9995	0.9981	−1.4278	−0.2901	0.9982
	20 dB	1st	0.9992	0.9988	0.9999	0.9992	0.9992
		2nd	0.9993	0.9969	0.8818	−2.0291	0.9976
		3rd	0.9998	0.9994	−1.6755	−0.1691	0.9986
Mean $R^2$	5 dB		0.9723	−0.6173	0.0138	0.0339	0.2029
	10 dB		0.9928	−0.5326	0.0378	0.0216	0.9844
	15 dB		0.9993	0.9973	−0.1279	0.2509	0.9958
	20 dB		0.9994	0.9984	0.0687	−0.3997	0.9985

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_k^i$ over time in example 2 are given in Figure 12. It can be seen that the proposed MKF algorithm has a good convergence.

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