Clutter Suppression and Rotor Blade Feature Extraction of a Helicopter Based on Time–Frequency Flash Shifts in a Passive Bistatic Radar

: This paper presents a passive bistatic radar (PBR) conﬁguration using a global navigation satellite system as an illuminator of opportunity for the rotor blade feature extraction of a helicopter. Aiming at the strong ﬁxed clutter in the surveillance channel of the PBR, a novel iteration clutter elimination method-based singular-value decomposition approach is proposed. Instead of the range elimination method used in the classic extended cancellation algorithm, the proposed clutter elimination method distinguishes the clutter using the largest singular value and by remove this value. At the same time, the fuselage echo of the hovering helicopter can also be suppressed along with the ground clutter, then the rotor echo of this can be obtained. In the micro-motion feature extraction, the mathematic principle of the ﬂash generation process in the time–frequency distribution (TFD) is derived ﬁrst. Next, the phase compensation method is applied to achieve the time–frequency ﬂash shift in the TFD. After this, the center frequencies of the standard ﬂashes in the TFD are compared with the standard frequency dictionary. The mean l1 norm is utilized to estimate the feature parameters of the helicopter rotor. In the experiments, the scattering point model and the physical optics facet model demonstrate that the proposed method can obtain more accurate parameter estimation results than some classic algorithms.


Introduction
Based on their illuminators, radar systems can be categorized as active or passive. Active radar systems utilize transmitters with known locations, while passive systems exploit illuminators of opportunity (IO) with limited knowledge [1,2]. In comparison to the active radar systems, the passive systems are low in cost, size, weight, and power; they exhibit lower manufacturing and operational costs due to the lack of emitting elements [3,4]. Furthermore, depending on the relative location between the transmitter and receiver, the radar systems can also be categorized as monostatic or bistatic [5]. The employed location of the IO is always known, but the receiver is always employed on the basis of the functional requirements, meaning a passive radar is always a passive bistatic radar (PBR). In recent years, the interest in exploiting global navigation satellite systems (GNSSs) for remote sensing has been growing on account of the advantages of its high coverage of the entire Earth [6]. For hovering targets, such as helicopters and drones at low altitude, it is difficult for ground radar systems to detect these targets, and airborne-based radars may act as a supplement for the blind areas to ground-based radars [7][8][9]. The exploitation of the GNSS for target recognition and imaging has been reported in [10]; however, the main limitation is the weak target echo in the surveillance channel.
Coherent integration (CI) is always implemented to enhance the target echo energy to achieve target detection, which means that the PBR contains two receiving channels: the implemented to deal with the clutter-suppressed surveillance signal. Afterwards, the short-time Fourier transform (STFT) is utilized to obtain the TFD. Then, the minimum mean l1 norm between the center frequency of the standard flashes in the TFD and the standard frequency is the principle used for extracting the feature parameters of the helicopter. In addition, the simulation experiment results are shown.
The remainder of paper is organized as follows. Section 2 introduces the GNSS-based PBR geometry, the signal model, and the mathematic mechanism of the flash generation. The proposed clutter suppression method and feature extraction algorithm based on the TFFS are described in Section 3. The simulation experiments are shown in Section 4. Lastly, Section 5 concludes the paper.

The Geometry of GNSS-Based PBR
The observation geometry of the GNSS-based PBR system is shown in Figure 1, where the transmitter and receiver are separated. The transmitter is one of the GNSS satellites in the earth orbit, and the receiver is located in origin O of the Cartesian coordinate system. The receiver consists of two channels, namely the reference channel and surveillance channel. t R and r R are the transmitter-to-target range and the receiver-to-target range, respectively. 0 R is the direct range between the transmitter and receiver. The solid bistatic angle between the helicopter and radar system is defined as . To remove the effect of the satellite motion, the Doppler frequency from the satellite motion can be estimated and compensated during the signal capture procedure, and the concrete procedure can be seen in [1]. After this, the transmitter of the bistatic radar can be regarded as fixed. As shown in Figure 1, not only does the surveillance channel acquire the target echo, but also the ground clutter and noise. The principal transmitted signal from the GNSS satellite is a code division multiple access (CDMA) code consisting of pseudo-random noise (PRN).
Of course, any known GNSS signal can be the signal basis for this PBR detection model. For clarity, this paper mainly concentrates on the Global Positioning System (GPS) L1 signal at 1.5 GHz, and the period of the continuous wave is 1 ms. In this space, the Tp = 1ms echo signal is intercepted as a whole for processing. Therefore, the continuous wave signal is modelled as the approximate pulse signal with a pulse width of Tp and pulse repetition frequency of prf = 1/Tp = 1 kHz [18]. The transmitted signal is set as s(t) and the sampled surveillance signal can be written as: To remove the effect of the satellite motion, the Doppler frequency from the satellite motion can be estimated and compensated during the signal capture procedure, and the concrete procedure can be seen in [1]. After this, the transmitter of the bistatic radar can be regarded as fixed. As shown in Figure 1, not only does the surveillance channel acquire the target echo, but also the ground clutter and noise. The principal transmitted signal from the GNSS satellite is a code division multiple access (CDMA) code consisting of pseudo-random noise (PRN).
Of course, any known GNSS signal can be the signal basis for this PBR detection model. For clarity, this paper mainly concentrates on the Global Positioning System (GPS) L1 signal at 1.5 GHz, and the period of the continuous wave is 1 ms. In this space, the T p = 1ms echo signal is intercepted as a whole for processing. Therefore, the continuous wave signal is modelled as the approximate pulse signal with a pulse width of T p and pulse repetition frequency of prf = 1/T p = 1 kHz [18]. The transmitted signal is set as s(t) and the sampled surveillance signal can be written as: S sur (n, m) = S r (n, m) + S c (n, m) + N(n, m) Atmosphere 2022, 13, 1214 4 of 20 where n and m denote the intercepted fast-time and slow-time domains, respectively; S r (n, m) denotes the target echo, which is weaker than the white Gaussian noise N(n, m); S c (n, m) is the ground clutter: where ρ, c, and λ denote the reflectivity, light velocity, and wavelength, respectively. Here, R ci denotes the bistatic range of the ground clutter; p r (n, m) represents the envelope function after the match filter. Because there is no relative motion between the ground and radar, the R ci is not varied from the slow-time domain. However, the amplitude of the ground clutter is the strongest in these three parts. Therefore, before the mD feature extraction, the ground clutter needs to be removed from the surveillance channel. One type of mD analysis method adapted to low signal-to-noise ratios (SNRs) has great significance.
According to the target model, the target echo S r (n, m) can be modelled as: where σ i (σ i << ρ), R bi are the reflectivity and bistatic-range of the ith scatterer within the target. The direct signal from reference channel can be modelled as: This paper discusses a situation where the helicopter is hovering. Only the rotor is rotating, and the other parts are nearly stationary. At this time, the bistatic range of the scatterers in the fuselage does not vary over time. Therefore, the fuselage echo can be modelled as: From this, it can be noted that the fuselage echo model is pretty similar to the ground clutter. Especially in the IO of the GNSS, the resolution cell is larger than the fuselage size, which means the whole fuselage echo will appear in one range cell. The fuselage echo can be written as follows: where γ denotes the composited scattering coefficient of the fuselage echo. Because these echoes are all from fixed scatterers, the ground clutter and fuselage echo are both named the fixed clutter in this paper. The fuselage echo can be eliminated together with the ground clutter using the fixed clutter elimination procedure. Figure 2 presents a view of the geometry of the rotor blades relative to the bistatic radar. Here, we set the velocity of the blade tip to V tip , with ϕ to the angle between V tip and the bisector of the bistatic angle β; w = 2π f denotes the rotation angular frequency and f r denotes the rotation frequency. For a rigid-body target, the rotation of the rotor blades leads to a variety of bistatic ranges R bi from the slow-time domain m:

The Rotor Echo of The Helicopter
Atmosphere 2022, 13, 1214 5 of 20 where R bi0 is the bistatic range of the rotation center, l i is the distance from the ith scatterer to the rotation center, and θ 0 is the initial rotation angle. and r f denotes the rotation frequency. For a rigid-body target, the rotation of the rotor blades leads to a variety of bistatic ranges bi R from the slow-time domain m: 0 0 cos cos cos 2 bi bi i R m R l wm (7) where 0 bi R is the bistatic range of the rotation center, i l is the distance from the ith scatterer to the rotation center, and 0 is the initial rotation angle. After integrating the echoes from all scattering points, the rotor echo can be expressed as [19]: cos cos cos 2 k k L m wm (9) where L is the length of the blade, 0 A denotes the amplitude determined by the radar cross-section of the kth blade, and sinc sin / x x x denotes the sinc function. We assume that a rotor consists of K blades that have K different initial angles The mD frequency shift can be obtained for each blade from Equation (9) (10) Owing to the rotation of the rotor, it is obvious that the Doppler frequency is modulated by the sinusoidal function. The maximum mD frequency is: max cos cos 22 After integrating the echoes from all scattering points, the rotor echo can be expressed as [19]:

Analysis of Flash Generation
where L is the length of the blade, A 0 denotes the amplitude determined by the radar cross-section of the kth blade, and sin cx = sin πx/πx denotes the sinc function. We assume that a rotor consists of K blades that have K different initial angles θ k = θ 0 + 2πk/K with k = 0, 1, 2, · · · , K − 1.
The mD frequency shift can be obtained for each blade from Equation (9) as: Owing to the rotation of the rotor, it is obvious that the Doppler frequency is modulated by the sinusoidal function. The maximum mD frequency is:

Analysis of Flash Generation
Without a loss of generality, from Equation (8), the normalized signal of the first blade is: It can be seen from Equation (12) that the signal profile includes several sinc peaks; namely, there are several flashes in the slow-time domain. To further analyze the variety of Doppler frequencies, the short-time Fourier transform (STFT) is applied to process the rotor echo. According to the characteristics of the Fourier transform, we know that: where FT{ } denotes the Fourier transform and ⊗ denotes the convolution operation.
We set the parameters as in Section 4.1, while the TFDs of the profile, the phase, and the entire signal of Equation (13) are shown in Figure 3. Due to the fact that sinc(t) and rect(f ) are a Fourier transform pair, Figure 3a shows that there are several flash stripes located at the zero Doppler frequency, which occupies a certain bandwidth. The time width of the stripe is caused by the low time resolution of the STFT. On the other hand, it can be seen from the electromagnetic scattering theorem that the specular reflectivity generates these  Figure 3b, it is clear that the Doppler frequency is modulated by the sinusoid function, as mentioned above. As described in Equation (13), Figure 3c is the convolution result of Figure 3a,b in the Doppler frequency domain. Instead of being located in the zero-Doppler frequency, the positive frequency flashes centered in the positive frequency and the negative frequency flashes centered in the negative frequency are alternating here.
where FT{ } denotes the Fourier transform and denotes the convolution operation. We set the parameters as in Section 4.1, while the TFDs of the profile, the phase, and the entire signal of Equation (13) are shown in Figure 3. Due to the fact that sinc(t) and rect(f) are a Fourier transform pair, Figure 3a shows that there are several flash stripes located at the zero Doppler frequency, which occupies a certain bandwidth. The time width of the stripe is caused by the low time resolution of the STFT. On the other hand, it can be seen from the electromagnetic scattering theorem that the specular reflectivity generates these flashes. From Figure 3b, it is clear that the Doppler frequency is modulated by the sinusoid function, as mentioned above. As described in Equation (13), Figure 3c is the convolution result of Figure 3a, b in the Doppler frequency domain. Instead of being located in the zero-Doppler frequency, the positive frequency flashes centered in the positive frequency and the negative frequency flashes centered in the negative frequency are alternating here. Furthermore, the mathematical principle for a pair of symmetric blades is similar to that for one blade. The Fourier transform is: Furthermore, the mathematical principle for a pair of symmetric blades is similar to that for one blade. The Fourier transform is: Differing from Equation (13), the phase item of Equation (14) includes a positive phase and a negative phase, as shown in Figure 3d. Hence, the positive frequency flashes and negative frequency flashes appear at the same time in the TFD, as shown in Figure 3e. In short, these flash stripes are produced using the sinc function, but the locations are determined by the echo phase.

The Iteration Elimination
In view of the fact that the clutter (including the direct signal, multipath clutter, and fuselage echo) is much stronger than the rotor echo and noise, the target echo is attributed to noise and the composite approximate noise is constructed. Then, the surveillance signal can be summarized as follows: where the composite approximate noiseN(n, m) is: Then, the constructed surveillance matrix S sur can be written as: where → s i represents the vector of the surveillance echo. The singular values of the surveillance matrix can respond to the weights of each component. It can be deduced that the larger singular values represent the principal component, while the remaining smaller singular values denote a weak rotor echo and noise. As shown in Equation (15), the surveillance signal is composed of clutter and the composite approximate noise, and the clutter power is greater than the composite noise, meaning the clutter can be regarded as the principal component of the surveillance matrix. The surveillance matrix S sur can be decomposed to obtain the singular values and related eigenvector matrixes, as shown in Equation (18): where [·] H represents the conjugate transpose. Assuming that the rank (S sur ) = P, then Σ = diag(η 1 , η 2 , η 3 , · · · , η P ), where η l denotes all non-zero singular values of the surveillance matrix S sur and satisfies η 1 > η 2 > η 3 > · · · > η P . According to the singularvalue decomposition theorem, the unitary matrix U is composed of all eigenvectors (u l , l = 1, 2, · · · , M) of S sur S H sur , and all eigenvectors (d l , l = 1, 2, · · · , N) of S H sur S sur form the unitary matrix D. Therefore, the surveillance matrix can be decomposed by: The first Q larger singular values η 1 , η 2 , η 3 , · · · , η Q represent the clutter component, and the clutter suppression can be achieved by eliminating the component relative to the larger singular value from the surveillance matrix S sur . It is known that the rotor echo power is far lower than the clutter, but the amount of multipath clutter and the number target characteristics are unknown. As a consequence, it is difficult to acquire the number of larger singular values in the initial stage. During concrete processing, the iteration elimination is implemented to remove the clutter, as shown in Equation (20): where [] (l) represents the iterations and the initial iteration matrix is set as S (0) sur = S sur . After the matrix decomposition, all singular values are sorted from largest to smallest. The search for the maximum singular value ratio is applied to determine the initial number of iterationsQ, as shown in Equation (21). The proportional coefficient coe is used to give a sufficient iteration allowance, which should be greater than 1.

The Optimal Clutter Elimination Method
Considering that the rotor echo is much lower than the noise, CI is usually implemented to improve the rotor echo power after the clutter suppression. In the processing, two-dimensional pulse compression is implemented to achieve CI. To achieve the optimal clutter elimination effectiveness, the information entropy process is applied to characterize the orderly state or chaos of the CI result, and the mathematical expression is written as in Equation (22). The smaller entropy corresponds to the less and more distinct peaks in the CI results. Firstly, during the ground clutter elimination procedure, the entropy of the CI result is decreased gradually. However, for helicopter target, when the ground clutter is removed completely, the entropy of the CI result begins increasing, since the rotor rotation leads to the sinusoid modulated Doppler frequency, namely the wideband Doppler frequency. There are several peaks in the Doppler frequency domain; hence, the entropy of the CI result will become larger-even larger than the origin. It should be emphasized that this method is applied to suppress the clutter and detect targets. Therefore, the optimal number of iterations is located on the turning point where the entropy of the CI result moves from down to up. The coefficient coe in Equation (21) causes the whole variation process of the entropy.
where En(0) represents the entropy of the CI result of the initial surveillance matrix. Additionally, the clutter suppression procedure is given in Figure 4.

Flash Shift in the Time-Frequency Domain
As mentioned in Section 2.3, the variety of flash locations is determined by the sinefunction-modulated phase of the rotor echo. In contrast, if the phase of the specific blade echo were compensated completely, those flash stripes of the blade would be centered on a certain Doppler frequency (zero-frequency for an asymmetrical rotor). Next, the analysis of the Doppler center frequency of the flashes is presented.

The Asymmetrical Rotor
The compensation phase is constructed as:

Flash Shift in the Time-Frequency Domain
As mentioned in Section 2.3, the variety of flash locations is determined by the sinefunction-modulated phase of the rotor echo. In contrast, if the phase of the specific blade echo were compensated completely, those flash stripes of the blade would be centered on a certain Doppler frequency (zero-frequency for an asymmetrical rotor). Next, the analysis of the Doppler center frequency of the flashes is presented.

The Asymmetrical Rotor
The compensation phase is constructed as: where (α, Ω, ξ) denote the amplitude, angular frequency, and initial angle, respectively. The problem of parameter estimation is transformed into an optimization problem. According to Equation (11), we set Θ = L/λ · cos(β/2) cos(ϕ), and Θ = 2 f max mD /w is obtained. Hence, it can be seen that the amplitude α is relative to the angular frequency Ω. B 1 = 2 f max mD denotes the bandwidth of the flash stripes, and it can be estimated in the TFD; thus, α = πΘ = πB 1 /w. Then, the three-dimensional parameter space (α, Ω, ξ) is decreases to (Ω, ξ), which can improve the parameter estimation efficiency. Furthermore, the phase-compensated equation (Equation (12)) is: It can be seen that when w = Ω and θ 1 = ξ, the phase of Equation (24) is compensated completely and the parameters ŵ,θ 1 of the blade are obtained. At the same time, Equation (24) is changed to: As noted earlier, the TFD of Equation (25) presents a series of flash stripes located in the zero Doppler frequency. In this manner, the flashes of the specific blade are all centered on the zero Doppler frequency after the phase compensation. Then, the rotor blade number K can be obtained directly from the flash number between nearest two zero-frequency flashes in the TFD.

The Symmetrical Rotor
Because of the symmetry of the rotor, each flash of the specific symmetric blade pair, including the positive frequency flash and negative frequency flash, is called a "double flash". The compensation procedure is as follows: where B 2 = 4 f max mD denotes the bandwidth of the double flash. At this time, when w = Ω and θ 1 = ξ, Equation (26) is changed to: The phase of the first item is compensated completely in Equation (26); thus, the flashes are all centered on the zero-Doppler frequency. However, the phase of the second item in Equation (26) is changed to the double initial, ignoring the effect of the minus sign. In the second item of Equation (27), the mD frequency is twice the initial value; hence, the flash stripes are centered on 2 f max mD or −2 f max mD . It can be noted that the two flashes of these two items are connected together, without overlap along the Doppler frequency axis. Then, the composited double flash of the specific symmetrical blade pair can only be centered on f max mD or − f max mD , and these two center frequencies are alternating. As before, all double flashes of the specific symmetrical blade pair after the phase compensation can be sought out by the centered frequency. Then, the number of symmetrical blade pairs can be obtained, and the blade number of the rotor can be obtained.

Optimization Criterion
The standard frequency (SF) that is the center frequency of the completely compensated flash for the different rotor structure is analyzed in Section 3.2. We found that the blade number of a rotor usually does not exceed 8, meaning the SF dictionary can be constructed as in Table 1. Table 1. The standard frequency dictionary.

Standard Flashes
The horizontal axis represents the standard flashes and their SF, while the vertical axis represents the blade number. For simplicity, the SF of double flash is described as a function of B 2 . For one blade, every flash (D 1 , D 2 , D 3 , · · ·) is regarded as a standard flash, and the SFs after the TFFS all equal zero. For two blades, all double flashes (D 1 , D 2 , D 3 , · · ·) are identical standard flashes, but the SF alternates between −B 2 /4 and B 2 /4. For three blades, there are two non-zero-frequency flashes of the other blades between the nearest two standard flashes (D 1 , −, −, D 4 , −, −, D 7 , −, · · ·), and the SFs are all zero. As for other cases, the distributions of the standard flashes and the SFs are shown in Table 1.
To obtain the optimal compensation result, the center frequencies (C 1 , C 2 , C 3 , · · · C I ) of all standard flashes are compared with the SF dictionary (D 1 , D 2 , D 3 , · · · D I ). The mean l 1 norm for the parameter space (Ω, ξ) is: where k denotes the blades number in the dictionary, I denotes the number of standard flashes, and | | denotes the l 1 norm. Since M(k, Ω, ξ) is a matrix that has the same dimensions as the parameter space (Ω, ξ) for each k, it is a three-dimensional matrix. Then, the blade numberK and ŵ,θ can be obtained by searching for the minimum of M(k, Ω, ξ) from the parameter space (k, Ω, ξ). Finally, the blade length can be calculated as: Additionally, the flowchart of the proposed mD feature extraction method is given in Figure 5.

The Scattering Point Model Simulation Experiments
To demonstrate the effectiveness of the proposed method, a helicopter was modelled for the GPS-based bistatic radar. After the satellite motion compensation phase during signal capture, the GPS satellite is regarded as stationary. The surveillance channel consists of the target echo, direct signal, and five paths of ground clutter, and the ground clutter parameters are shown in Table 2. The helicopter is a model of a fuselage and rotor; the fuselage is modelled as uniform scatters and the rotor is modelled as three blades with a length L = 7.3 m, as well as a perfect electrical conductor. The blades are assumed to rotate at a rotation frequency fr = 4.8 r/s and to be illuminated by a simulated GPS L1 signal in the location where the bistatic factor cos 2 cos 0.2767 , as shown in Figure 1. The sampling rate of the employed receiver is set to 2.046MHz. Finally, the simulation experiments are implemented under SNR = −25 dB, namely the SNR of the raw helicopter target echo before the CI is −25 dB. After sampling signal, the raw surveillance signal is constructed as a surveillance matrix, and the CI result of the surveillance matrix is shown in Figure 6. It is obvious that there is a strong peak indicating the direct signal in (0.4399 Km, 0.6676 Hz). For the sake of observation, the direct signal is shifted in the bi-range of 0.4399 Km in this simulation

The Scattering Point Model Simulation Experiments
To demonstrate the effectiveness of the proposed method, a helicopter was modelled for the GPS-based bistatic radar. After the satellite motion compensation phase during signal capture, the GPS satellite is regarded as stationary. The surveillance channel consists of the target echo, direct signal, and five paths of ground clutter, and the ground clutter parameters are shown in Table 2. The helicopter is a model of a fuselage and rotor; the fuselage is modelled as uniform scatters and the rotor is modelled as three blades with a length L = 7.3 m, as well as a perfect electrical conductor. The blades are assumed to rotate at a rotation frequency f r = 4.8 r/s and to be illuminated by a simulated GPS L1 signal in the location where the bistatic factor cos(β/2) cos ϕ = 0.2767, as shown in Figure 1. The sampling rate of the employed receiver is set to 2.046 MHz. Finally, the simulation experiments are implemented under SNR = −25 dB, namely the SNR of the raw helicopter target echo before the CI is −25 dB. After sampling signal, the raw surveillance signal is constructed as a surveillance matrix, and the CI result of the surveillance matrix is shown in Figure 6. It is obvious that there is a strong peak indicating the direct signal in (0.4399 Km, 0.6676 Hz). For the sake of observation, the direct signal is shifted in the bi-range of 0.4399 Km in this simulation experiment. The Doppler frequency of the direct signal would be zero after the signal capture and satellite motion compensation phases, but the spectrum leakage brings the Doppler frequency to 0.6676 Hz. From this, it is known that only the direct signal can be detected from the CI result for the raw surveillance matrix, and the rotor echoes are submerged in strong cluttering and noise. Therefore, the clutter suppression is necessary before the micro-motion feature extraction.
Atmosphere 2022, 13, x FOR PEER REVIEW 13 of 21 experiment. The Doppler frequency of the direct signal would be zero after the signal capture and satellite motion compensation phases, but the spectrum leakage brings the Doppler frequency to 0.6676 Hz. From this, it is known that only the direct signal can be detected from the CI result for the raw surveillance matrix, and the rotor echoes are submerged in strong cluttering and noise. Therefore, the clutter suppression is necessary before the micro-motion feature extraction. Figure 6. The CI of raw echoes.
Next, the raw surveillance matrix is processed using the proposed iteration elimination method. Firstly, the singular decomposition phase is implemented, and the change in the singular value ratio is shown in Figure 7. For a clearer view, only the first 50 singular ratio values are shown here. From the image, it is clear that the singular value ratio presents a procedure of increasing and then decreasing, so the inflection point is the focus. The 7th singular ratio value is the maximum, which means that the variation between the 7th and the 8th singular values is the most dramatic. Namely, the first seven larger singular values denote the strong clutter component that needs to be eliminated, and the other smaller singular values indicate the rotor echo component and noise. For enough iteration margin, the factor coe is set to 2 in this experiment, and then the initial iterations Q can be set to 14. As mentioned above, the strong clutter will be eliminated by the iterations. The optimal number of iterations is determined by the entropy variation of the CI results, and the entropy variation during the iterations in the experiment is shown in Figure 8. It can be seen that entropy variation presents a decreasing and then increasing trend with the increase in iterations. From the 1st to 6th iterations, the entropy is decreasing gradually, then the entropy reaches the minimum at the 6th iteration. The reason is that the fixed clutter is eliminated gradually during this process, whereby the clutter peaks are cut down. After 7 iterations, the fixed clutter is removed completely, and the surveillance signal only includes the rotor echo. Differing from the normal motion target, the rotor includes several rotating blades. The Doppler frequency of the rotor echo is modulated by the sinusoid function, and there are several peaks in the Doppler frequency domain after the CI. Hence, the entropy of the CI results of the rotor echo is much larger. However, the Next, the raw surveillance matrix is processed using the proposed iteration elimination method. Firstly, the singular decomposition phase is implemented, and the change in the singular value ratio is shown in Figure 7. For a clearer view, only the first 50 singular ratio values are shown here. From the image, it is clear that the singular value ratio presents a procedure of increasing and then decreasing, so the inflection point is the focus. The 7th singular ratio value is the maximum, which means that the variation between the 7th and the 8th singular values is the most dramatic. Namely, the first seven larger singular values denote the strong clutter component that needs to be eliminated, and the other smaller singular values indicate the rotor echo component and noise. For enough iteration margin, the factor coe is set to 2 in this experiment, and then the initial iterationsQ can be set to 14.
Atmosphere 2022, 13, x FOR PEER REVIEW 13 of 21 experiment. The Doppler frequency of the direct signal would be zero after the signal capture and satellite motion compensation phases, but the spectrum leakage brings the Doppler frequency to 0.6676 Hz. From this, it is known that only the direct signal can be detected from the CI result for the raw surveillance matrix, and the rotor echoes are submerged in strong cluttering and noise. Therefore, the clutter suppression is necessary before the micro-motion feature extraction. Figure 6. The CI of raw echoes.
Next, the raw surveillance matrix is processed using the proposed iteration elimination method. Firstly, the singular decomposition phase is implemented, and the change in the singular value ratio is shown in Figure 7. For a clearer view, only the first 50 singular ratio values are shown here. From the image, it is clear that the singular value ratio presents a procedure of increasing and then decreasing, so the inflection point is the focus. The 7th singular ratio value is the maximum, which means that the variation between the 7th and the 8th singular values is the most dramatic. Namely, the first seven larger singular values denote the strong clutter component that needs to be eliminated, and the other smaller singular values indicate the rotor echo component and noise. For enough iteration margin, the factor coe is set to 2 in this experiment, and then the initial iterations Q can be set to 14. As mentioned above, the strong clutter will be eliminated by the iterations. The optimal number of iterations is determined by the entropy variation of the CI results, and the entropy variation during the iterations in the experiment is shown in Figure 8. It can be seen that entropy variation presents a decreasing and then increasing trend with the increase in iterations. From the 1st to 6th iterations, the entropy is decreasing gradually, then the entropy reaches the minimum at the 6th iteration. The reason is that the fixed clutter is eliminated gradually during this process, whereby the clutter peaks are cut down. After 7 iterations, the fixed clutter is removed completely, and the surveillance signal only includes the rotor echo. Differing from the normal motion target, the rotor includes several rotating blades. The Doppler frequency of the rotor echo is modulated by the sinusoid function, and there are several peaks in the Doppler frequency domain after the CI. Hence, the entropy of the CI results of the rotor echo is much larger. However, the As mentioned above, the strong clutter will be eliminated by the iterations. The optimal number of iterations is determined by the entropy variation of the CI results, and the entropy variation during the iterations in the experiment is shown in Figure 8. It can be seen that entropy variation presents a decreasing and then increasing trend with the increase in iterations. From the 1st to 6th iterations, the entropy is decreasing gradually, then the entropy reaches the minimum at the 6th iteration. The reason is that the fixed clutter is eliminated gradually during this process, whereby the clutter peaks are cut down. After 7 iterations, the fixed clutter is removed completely, and the surveillance signal only includes the rotor echo. Differing from the normal motion target, the rotor includes several rotating blades. The Doppler frequency of the rotor echo is modulated by the sinusoid function, and there are several peaks in the Doppler frequency domain after the CI. Hence, the entropy of the CI results of the rotor echo is much larger. However, the rotor echo also fades gradually after 8 iterations. Because the rotor consists of many scatterers, the first several iterations will not have an evident impact on the CI result, and the entropy remains nearly constant after the 7th iteration. To guarantee the integrity of the rotor echo, the optimal number of iterations is selected as 7. Then, the clutter-suppressed surveillance matrix is processed via CI, and the result is shown in Figure 9. There is a peak at about 7 Km that denotes the rotor echo in Figure 9a, and the strong clutter at 0.4399 Km is removed. Figure 9b indicates the Doppler frequency-amplitude plane of the CI result, where the rotor echo occupies a certain bandwidth. As discussed above, the wideband signal is produced by the rotational rotor. Therefore, after clutter suppression, the helicopter or the rotor can be detected and the rotor echo in the slow-time domain can be obtained.
rotor echo also fades gradually after 8 iterations. Because the rotor consists of many scatterers, the first several iterations will not have an evident impact on the CI result, and the entropy remains nearly constant after the 7th iteration. To guarantee the integrity of the rotor echo, the optimal number of iterations is selected as 7. Then, the cluttersuppressed surveillance matrix is processed via CI, and the result is shown in Figure 9. There is a peak at about 7 Km that denotes the rotor echo in Figure 9a, and the strong clutter at 0.4399 Km is removed. Figure 9b indicates the Doppler frequency-amplitude plane of the CI result, where the rotor echo occupies a certain bandwidth. As discussed above, the wideband signal is produced by the rotational rotor. Therefore, after clutter suppression, the helicopter or the rotor can be detected and the rotor echo in the slowtime domain can be obtained. To illustrate the effectiveness of the clutter elimination method proposed in this paper, the detection property is compared with the ECA. Firstly, each pulse is processed using the ECA with the direct signal. Then, the CI result is shown in Figure 10. The birange of the target is 7.24 Km, but the bi-range of multipath 2 exceeds 12 Km. Accordingly, the ECA is applied to suppress the fixed clutter within 15 Km in this experiment. From Figure 10a, it can be seen that the fixed clutter within 15 Km is suppressed entirely. Figure  10b shows the Doppler frequency-amplitude plane, although there is no wideband spectral illustration as in Figure 9b. It is obvious that the whole rotor echo is also removed. The ECA cannot be used for helicopter detection in this experiment. rotor echo also fades gradually after 8 iterations. Because the rotor consists of many scatterers, the first several iterations will not have an evident impact on the CI result, and the entropy remains nearly constant after the 7th iteration. To guarantee the integrity of the rotor echo, the optimal number of iterations is selected as 7. Then, the cluttersuppressed surveillance matrix is processed via CI, and the result is shown in Figure 9. There is a peak at about 7 Km that denotes the rotor echo in Figure 9a, and the strong clutter at 0.4399 Km is removed. Figure 9b indicates the Doppler frequency-amplitude plane of the CI result, where the rotor echo occupies a certain bandwidth. As discussed above, the wideband signal is produced by the rotational rotor. Therefore, after clutter suppression, the helicopter or the rotor can be detected and the rotor echo in the slowtime domain can be obtained. To illustrate the effectiveness of the clutter elimination method proposed in this paper, the detection property is compared with the ECA. Firstly, each pulse is processed using the ECA with the direct signal. Then, the CI result is shown in Figure 10. The birange of the target is 7.24 Km, but the bi-range of multipath 2 exceeds 12 Km. Accordingly, the ECA is applied to suppress the fixed clutter within 15 Km in this experiment. From Figure 10a, it can be seen that the fixed clutter within 15 Km is suppressed entirely. Figure  10b shows the Doppler frequency-amplitude plane, although there is no wideband spectral illustration as in Figure 9b. It is obvious that the whole rotor echo is also removed. The ECA cannot be used for helicopter detection in this experiment. To illustrate the effectiveness of the clutter elimination method proposed in this paper, the detection property is compared with the ECA. Firstly, each pulse is processed using the ECA with the direct signal. Then, the CI result is shown in Figure 10. The bi-range of the target is 7.24 Km, but the bi-range of multipath 2 exceeds 12 Km. Accordingly, the ECA is applied to suppress the fixed clutter within 15 Km in this experiment. From Figure 10a, it can be seen that the fixed clutter within 15 Km is suppressed entirely. Figure 10b shows the Doppler frequency-amplitude plane, although there is no wideband spectral illustration as in Figure 9b. It is obvious that the whole rotor echo is also removed. The ECA cannot be used for helicopter detection in this experiment. In brief, the ECA has an intrinsic limitation. Because the ECA achieves target detection by removing the short bi-range ground clutter, the ECA cannot achieve helicopter target detection when the bi-range of the target is smaller than that of the clutter. However, the proposed iteration elimination method is aimed at the echo energy, In brief, the ECA has an intrinsic limitation. Because the ECA achieves target detection by removing the short bi-range ground clutter, the ECA cannot achieve helicopter target detection when the bi-range of the target is smaller than that of the clutter. However, the proposed iteration elimination method is aimed at the echo energy, meaning it has a better effect.

The Micro-Motion Feature Extraction
After the fixed clutter suppression, the rotor echo in the slow-time domain can be extracted from the clutter-suppressed surveillance matrix. Figure 11 gives the range profile for this, which shows that the rotor echo mainly consists of several sinc peaks, as described in Section 2.2. In other words, there are several iterations of specular reflection between the target and bistatic radar. For a further analysis, Figure 12 presents the TFD of the extracted rotor echo by utilizing the STFT. It is different from Figure 11 in that there are alternating positive flashes and negative flashes in the TFD. The bandwidth B 1 of each flash is about 286.6 Hz. As discussed above, the asymmetry of the flashes in the TFD is caused by the asymmetry of the rotor. Therefore, it can be seen that the rotor is comprised of odd blades. The dimension of the SF dictionary is decreased to four (1, 3, 5, 7). The concrete dictionary for this experiment is shown in Table 3. From this, the different standard flashes are utilized for the different blade numbers.
(a) (b) Figure 10. The CI after the ECA: (a) the range-amplitude plane; (b) the Doppler frequencyamplitude plane.
In brief, the ECA has an intrinsic limitation. Because the ECA achieves target detection by removing the short bi-range ground clutter, the ECA cannot achieve helicopter target detection when the bi-range of the target is smaller than that of the clutter. However, the proposed iteration elimination method is aimed at the echo energy, meaning it has a better effect.

The Micro-Motion Feature Extraction
After the fixed clutter suppression, the rotor echo in the slow-time domain can be extracted from the clutter-suppressed surveillance matrix. Figure 11 gives the range profile for this, which shows that the rotor echo mainly consists of several sinc peaks, as described in Section 2.2. In other words, there are several iterations of specular reflection between the target and bistatic radar. For a further analysis, Figure 12 presents the TFD of the extracted rotor echo by utilizing the STFT. It is different from Figure 11 in that there are alternating positive flashes and negative flashes in the TFD. The bandwidth B1 of each flash is about 286.6 Hz. As discussed above, the asymmetry of the flashes in the TFD is caused by the asymmetry of the rotor. Therefore, it can be seen that the rotor is comprised of odd blades. The dimension of the SF dictionary is decreased to four (1, 3, 5, 7). The concrete dictionary for this experiment is shown in Table 3. From this, the different standard flashes are utilized for the different blade numbers.    Figure 13 presents the minimum ( , , ) Mk values for different blade numbers k (1, 3, 5, 7). In this image, a minimum can be found when K = 3, which means the rotor consists of 3 blades. The variety of (3, , ) M values is given in Figure 14. As presented in Section 3.3, the (3, , ) M is the same dimensional matrix as the parameter space. From the Figure 12. The TFD of the raw echo. Table 3. The standard frequency dictionary in this experiment.  3, 5, 7). In this image, a minimum can be found when K = 3, which means the rotor consists of 3 blades. The variety of M(3, Ω, ξ) values is given in Figure 14. As presented in Section 3.3, the M(3, Ω, ξ) is the same dimensional matrix as the parameter space. From the image, it is obvious that M(3, Ω, ξ) reaches the minimum value 2.405 whenf r = 4.8 r/s andθ = 5.76 rad. Based on the estimated parameters, the length of the blades can be calculated as in Equation (29), andL = 7.24 m. Otherwise, the TFD after the complete phase compensation procedure is as shown in Figure 15. There are a series of flashes centering on the zero frequency, as illustrated by the white arrows, which belong to the same blade. It can be seen that there are two non-zero-frequency flashes from other blades between the nearest two zero-frequency flashes, and consequently the rotor structure includes 3 blades. The result from Figure 15 is coincident with that of M(k, Ω, ξ). The above coincident estimation results demonstrate the effectiveness of the proposed method. Figure 13 presents the minimum ( , , ) Mk values for different blade numbers k (1,3,5,7). In this image, a minimum can be found when K = 3, which means the rotor consists of 3 blades. The variety of (3, , ) M values is given in Figure 14. As presented in Section 3. 3 . Based on the estimated parameters, the length of the blades can be calculated as in Equation (29), and m 7.24 L . Otherwise, the TFD after the complete phase compensation procedure is as shown in Figure 15. There are a series of flashes centering on the zero frequency, as illustrated by the white arrows, which belong to the same blade. It can be seen that there are two non-zero-frequency flashes from other blades between the nearest two zero-frequency flashes, and consequently the rotor structure includes 3 blades. The result from Figure 15 is coincident with that of ( , , ) Mk

Blade Number
. The above coincident estimation results demonstrate the effectiveness of the proposed method.   In order to validate the anti-noise ability of the proposed method, 100-iteration Monte Carlo simulation experiments are conducted, and Figure 16 gives the root mean square errors (RMSEs) of the estimated parameters with the increasing SNRs. In Figure  16, the simulation parameters are the same as above, but the SNR increases from −50 dB to 10 dB with the step size of 3 dB. For comparison, the processing results for the classic  In order to validate the anti-noise ability of the proposed method, 100-iteration Monte Carlo simulation experiments are conducted, and Figure 16 gives the root mean square errors (RMSEs) of the estimated parameters with the increasing SNRs. In Figure  16, the simulation parameters are the same as above, but the SNR increases from −50 dB In order to validate the anti-noise ability of the proposed method, 100-iteration Monte Carlo simulation experiments are conducted, and Figure 16 gives the root mean square errors (RMSEs) of the estimated parameters with the increasing SNRs. In Figure 16, the simulation parameters are the same as above, but the SNR increases from −50 dB to 10 dB with the step size of 3 dB. For comparison, the processing results for the classic OMP, Hough transform, and GLRT algorithm are also presented. Figure 15. The optimal shifted TFD.
In order to validate the anti-noise ability of the proposed method, 100-iteration Monte Carlo simulation experiments are conducted, and Figure 16 gives the root mean square errors (RMSEs) of the estimated parameters with the increasing SNRs. In Figure  16, the simulation parameters are the same as above, but the SNR increases from −50 dB to 10 dB with the step size of 3 dB. For comparison, the processing results for the classic OMP, Hough transform, and GLRT algorithm are also presented. Firstly, it can be seen from Figure 16a that these four methods can all obtain accurate the blade number for the rotors when the SNR is greater than −40 dB. Even the TFFS, HT, and GLRT methods can obtain the correct blade number when the SNR is as low as −47 dB. However, Figure 16a illustrates that the OMP algorithm cannot obtain the true blade number when the SNR is smaller than −40 dB. Accordingly, the proposed TFFS algorithm Firstly, it can be seen from Figure 16a that these four methods can all obtain accurate the blade number for the rotors when the SNR is greater than −40 dB. Even the TFFS, HT, and GLRT methods can obtain the correct blade number when the SNR is as low as −47 dB. However, Figure 16a illustrates that the OMP algorithm cannot obtain the true blade number when the SNR is smaller than −40 dB. Accordingly, the proposed TFFS algorithm is more robust than the OMP algorithm at low SNRs. Figure 16b,c show the RMSEs of the rotation frequency and blade length, respectively, with increasing SNRs. In Figure 16b, it is obvious that the purple diamond broken curve is lower than other three curves for all SNRs. The rotation frequency RMSEs of the TFFS and HT methods are pretty close. These three estimation results are nearly identical when the SNR is greater than −20 dB. It is believed that the estimation results for the GLRT are slightly better than the HT and the proposed TFFS method for the rotation frequency is better at low SNRs. Figure 16c shows that the red square line is lower than other three curves, and the purple diamond broken line is close to 1 m. Therefore, it is agreed that the proposed TFFS method can obtain accurate blade length values, but the blade length RMSE of the GLRT method is greater. The main reason is the effect of the flashes in the slow-time domain, whereby the GRLT cannot estimate the true blade length. Overall, the proposed TFFS algorithm can accurately extract the mD feature, and it is more robust than the classic OMP, HT, and GLRT algorithms.

The PO Facet Model Experiments
To validate the effectiveness of the proposed mD feature extraction method, the TFFS method is tested using the electromagnetic calculation data. Hence, this experiment assumes that the fixed clutter has been eliminated. The physical optics (PO) method is applied to generate the rotor echoes in the bistatic radar, and this experiment will consider the shield and rebound of the electromagnetic wave. Therefore, the PO facet model is designed to simulate the work situation of the real passive bistatic radar.
Based on the PO facet radar cross-section model, a rotor model can be represented by the arrays of triangular facets. The scatterer center of each triangle is assumed to be the geometric centroid of its triangle vertices, and the simulated three-blade rotor model is shown in Figure 17. The bistatic angle of the simulated model is 150 • , which is the sum of 70 • (incident elevation angle) and 80 • (reception elevation angle). The larger size of the rotor will lead to greater calculation complexity. To reduce the calculation complexity, the length and width of the target are changed to 3.9 m and 0.2 m, respectively, and the carrier frequency is also 1.5 GHz.
Based on the PO facet radar cross-section model, a rotor model can be rep by the arrays of triangular facets. The scatterer center of each triangle is assumed geometric centroid of its triangle vertices, and the simulated three-blade rotor m shown in Figure 17. The bistatic angle of the simulated model is 150°, which is th 70° (incident elevation angle) and 80° (reception elevation angle). The larger siz rotor will lead to greater calculation complexity. To reduce the calculation comple length and width of the target are changed to 3.9 m and 0.2 m, respectively, and th frequency is also 1.5 GHz.   Figure 18a. The reason is that the scatter i varies in different parts in the rotor. Otherwise, as a result of the incline of the obs model, the peak intervals are also not homogeneous. In the TFD, it can be seen flash intensity and interval are also varied, as presented in Figure 18b. A observation can reveal that the six flashes can be regarded as a rotation cycle that two strong flashes and four weak flashes, as indicated in the white circle in Fig  The processing results of the proposed TFFS method are given in Figure 19. Fi demonstrates that the estimated rotation frequency is 4.8 r/s. Furthermore, th length can be obtained as written in Table 4. The optimal TFD after the TFFS pr is shown in Figure 19b. As indicated in Figure 19b, the white arrow marks t  Figure 18 presents the simulated echo from the PO electromagnetic calculation method. Differing from the scattering point model, the scatter intensity of the peak varies severely from the slow-time domain in Figure 18a. The reason is that the scatter intensity varies in different parts in the rotor. Otherwise, as a result of the incline of the observation model, the peak intervals are also not homogeneous. In the TFD, it can be seen that the flash intensity and interval are also varied, as presented in Figure 18b. A closer observation can reveal that the six flashes can be regarded as a rotation cycle that includes two strong flashes and four weak flashes, as indicated in the white circle in Figure 18b. The processing results of the proposed TFFS method are given in Figure 19. Figure 19a demonstrates that the estimated rotation frequency is 4.8 r/s. Furthermore, the blade length can be obtained as written in Table 4. The optimal TFD after the TFFS processing is shown in Figure 19b. As indicated in Figure 19b, the white arrow marks the zero-frequency flashes belonging to the same blade. It is obvious that there are two non-zero-frequency flashes from other blades between the nearest two zero-frequency flashes. Therefore, it can be concluded that the rotor consists of three blades. Furthermore, the strength of the nearest two zero-frequency flashes shows a great change. The reason is that the reflectivity values of the blade moving toward and that of the same blade going away are not identical in the bistatic radar.
Atmosphere 2022, 13, x FOR PEER REVIEW 19 of 21 frequency flashes belonging to the same blade. It is obvious that there are two non-zerofrequency flashes from other blades between the nearest two zero-frequency flashes. Therefore, it can be concluded that the rotor consists of three blades. Furthermore, the strength of the nearest two zero-frequency flashes shows a great change. The reason is that the reflectivity values of the blade moving toward and that of the same blade going away are not identical in the bistatic radar.

2
For comparison, the processing results for the OMP, HT, and GLRT methods are also given in Figure 20, and the parameter estimation results are all written in Table 4. It can be seen that there are several sidelobes around the three main lobes in Figure 20a. In addition, the sidelobes and main lobe are almost located at an identical blade length or initial angle. Thus, we can ignore these sidelobes and only consider the three main lobes in the OMP parameter estimation results. The rotor consists of three blades. Only two obvious peaks have stronger reflectivity values that can be observed in Figure 20b; namely, the HT method processing result only detects two blades, and the rotation frequency is about 4.8 r/s. Additionally, the processing result for the GLRT method also includes two strong peaks, as shown in Figure 20c. In other words, the GLRT method only strong flashes weak flashes The same blade Figure 19. The processing result of the TFFS: (a) parameter estimation results of the TFFS; (b) the optimal shifted TFD. For comparison, the processing results for the OMP, HT, and GLRT methods are also given in Figure 20, and the parameter estimation results are all written in Table 4. It can be seen that there are several sidelobes around the three main lobes in Figure 20a. In addition, the sidelobes and main lobe are almost located at an identical blade length or initial angle. Thus, we can ignore these sidelobes and only consider the three main lobes in the OMP parameter estimation results. The rotor consists of three blades. Only two obvious peaks have stronger reflectivity values that can be observed in Figure 20b; namely, the HT method processing result only detects two blades, and the rotation frequency is about 4.8 r/s. Additionally, the processing result for the GLRT method also includes two strong peaks, as shown in Figure 20c. In other words, the GLRT method only detects two blades of the rotor. The principal reason for the worse processing results is the inhomogeneous nature of the scatter intensity. According to Table 4, only the proposed TFFS and OMP methods can detect the three blades of the rotor, and the estimation result of the TFFS method is more accurate than that of the OMP method.
Atmosphere 2022, 13, x FOR PEER REVIEW 20 of 21 detects two blades of the rotor. The principal reason for the worse processing results is the inhomogeneous nature of the scatter intensity. According to Table 4, only the proposed TFFS and OMP methods can detect the three blades of the rotor, and the estimation result of the TFFS method is more accurate than that of the OMP method. With these classic methods, the number of peaks in the parameter search result determines the blade number. The inhomogeneous reflectivity of the rotor can generate the different acuity peaks; hence, the acuity of the peaks can also affect the parameter estimation accuracy. However, the parameter search result for the proposed method only includes one valley, and the blade number can be obtained from the optimal TFD directly. Consequently, in the PO facet model simulation experiments, the proposed TFFS can also obtain more accurate parameter estimation results. With these classic methods, the number of peaks in the parameter search result determines the blade number. The inhomogeneous reflectivity of the rotor can generate the different acuity peaks; hence, the acuity of the peaks can also affect the parameter estimation accuracy. However, the parameter search result for the proposed method only includes one valley, and the blade number can be obtained from the optimal TFD directly. Consequently, in the PO facet model simulation experiments, the proposed TFFS can also obtain more accurate parameter estimation results.

Conclusions
Focusing on the micro-motion feature extraction process for the GNSS-based PBR, this paper proposed a fixed clutter suppression method based on iteration elimination and a feature extraction method based on the TFFS method for a helicopter rotor. In the proposed method, the SVD is utilized to disassemble the clutter component and the rotor echo component in the surveillance signal, and the iterations of the clutter component corresponding to the larger singular values are eliminated gradually. In comparison to the ECA, the iteration elimination method can be sued for helicopter detection and rotor echo extraction when the bi-range of the target is shorter than that of the clutter. Next, the TFFS algorithm was applied to shift the standard flashes to a certain frequency. Then, the feature parameters of the rotor were extracted by minimizing the mean l 1 norm M(k, Ω, ξ). Finally, the scattering point model and PO facet model simulation experiments were performed to validate the effectiveness and robustness of the new method. Compared with the classic OMP, HT, and GLRT methods, this method can achieve more accurate and robust rotor feature extraction results for a helicopter.