Imaging and Motion Parameter Estimation of Flying Helicopter With Duo Airborne SARs in $X$-Band

Duo airborne synthetic aperture radars in the $X$-band are proposed to acquire images of a flying helicopter. The backscattered signals manifest entangled features attributed to the compound linear motion of fuselage and spin motion of rotors. A two-stage template matching method is proposed to meticulously untangle signals from moving fuselage and spinning main rotors and tail rotor. The velocity vector of fuselage and the spinning rate of main rotors can be accurately estimated, and the spinning rate of the tail rotor can be estimated after removing the strong signals from fuselage and main rotors. The acquired images in different scenarios are evaluated and verified with three state-of-the-art algorithms.


I. INTRODUCTION
H ELICOPTERS can sneak by different terrains via agile maneuvering. The detection and identification of flying helicopters have been studied for decades [1], [2], [3], by using passive radars [4], [5], [6], signals of opportunity [7], [8], infrared, synthetic aperture radar (SAR), and inverse SAR (ISAR) [3]. Conventional detection methods on moving targets mainly focused on their translational motion; it takes more technical challenges to retrieve the parameters of spinning rotors for the identification of a helicopter model. SAR has been widely used to acquire the image of a moving target and estimate its motion parameters [9], [10]. However, SAR imaging on flying helicopters was rarely discussed in the literature. In this article, duo airborne SARs in the X-band are proposed to acquire images of a flying helicopter and to estimate its velocity and spinning parameters.
The rotation or vibration of objects induces micro-Doppler signatures, widening the spectrum centered at the Doppler frequency [11], [12]. The micro-Doppler signatures can be utilized to identify objects with spinning parts, such as the fast spinning blades of drones [12]. The micro-Doppler features of spinning rotors make helicopters distinct from fix-winged aircraft [4], [5]. Radar signals backscattered from a flying helicopter's fuselage, rotor blades, rotor mast, and tail rotor [1] are typically dominated by that from the fuselage [3]. The fuselage is easier to conceal with stealth technology, but the micro-Doppler features are not easy to cover [13]. Images on different parts of a flying helicopter can be acquired from their corresponding backscattered signals, which can be meticulously extracted by exploiting their specific motion features.
The micro-Doppler features of rotors depend on the number of blades, spin rate, and the blade length [4]. Moving target indication radar used to detect aircraft via Doppler signatures can easily miss the relatively feeble micro-Doppler features of a flying helicopter [14]. To manifest the weak micro-Doppler features, stronger echoes from fuselage must be removed first, which can be done by using short-time Fourier transform [2], Wigner-Ville transform [11], Choi-Williams transform [11], wavelet transform [15], or CLEAN method [3]. In this article, the echoes from fuselage are removed by applying a series of inverse processes on the coarse SAR image after compensating the linear motion of the helicopter. Then, an angular spectrum is proposed for the first time to estimate the spin rates of the main rotors.
ISAR methods have been used to acquire the image of rotating objects [16], [17], [18], [19]. In [20], a modified cubic phase function was proposed for ISAR imaging with quadratic frequency-modulated signals, and the angular velocity of a slowly rotating target could be estimated. The range model of a spinning rotor blade and its micro-Doppler signature were analyzed in [21]. In [16], the effect of rotation motion on imaging was analyzed and compensated by using an optimization method. Methods of subaperture, subpatch, polar-formatting, and backprojection were applied to deal with a larger rotation angle [17]. In [18], a segmental pseudo Keystone transform was proposed to compensate the motion of fast rotating targets, where a target was modeled as a bunch of scatterers, and a spinning scatterer was focused by exploiting its spin rate and initial phase.
Various methods have been proposed to estimate the rotational motion parameters. In [22], an extended Hough transform was applied to detect the underlying sinusoid signature attributed to spinning rotors. In [19], a 3-D ISAR method on rapidly spinning targets was proposed, by developing a matched filter bank from the signal features to generate a series of 2-D image slices. In [23], Wigner-Ville distribution was used to estimate the spin rate of a target from its millimeter-wave returned signals.
Inspired by these works, we propose a configuration of duo airborne SAR to acquire the image of a flying helicopter and estimate its motion parameters. An effective procedure is developed and verified by simulations. is integrated with a range-frequency reversal transform (RFRT) to estimate two Doppler centroid frequencies from the duo SARs for estimating the target velocity. Strong signals backscattered from the fuselage need to be removed to reveal weak signals from rotors. A two-stage template matching method is proposed to meticulously separate signals backscattered from fuselage, main rotors, and tail rotor, respectively. The spin rate and the initial phase of the main rotors are estimated by using an angular spectrum, which is proposed for the first time. The difficulties and solutions to estimate the spin rate of the tail rotor for acquiring its image are elaborated and verified by simulations for the first time, which were neglected in the literature. With velocity vector and spinning parameters accurately estimated, the image can be well focused. The acquired SAR images in different scenarios are evaluated and verified with three state-of-the-art (SOTA) algorithms.
The rest of this article is organized as follows. The simulation scenarios, parameters, and backscattered signals from fuselage, main rotors, and tail rotor are elaborated in Section II. The estimation of linear motion parameters and velocity vector of the flying helicopter, as well as a coarse image of the fuselage, is presented in Section III. The estimation on spin parameters of main rotors and tail rotor is elaborated in Section IV. In Section V, the images of helicopter are acquired in scenarios with various combinations of the velocity vector and the spin rate. The images are evaluated and verified with three SOTA algorithms. Analysis on computational load and highlight of contributions are also presented. Finally, Section VI concludes this article. Fig. 1 shows the configuration of duo SAR for imaging and estimating the velocity vector of a point target. The three black arrows with starting points A 1 , A 2 , and A 3 mark the flight paths of the first airborne platform P 1 , the second airborne platform P 2 , and helicopter Q, respectively. The three gray arrows with starting points A 1 , A 2 , and A 3 are the projections of the three black arrows, respectively, on the ground. Both platforms P 1 and P 2 fly at designate altitudes, with their positions and velocities

A. Default Parameters
Both airborne SAR platforms P 1 and P 2 can reach an altitude of 12 km and the maximum speed of 0.82 mach (279 m/s) [24]. Table I lists the default parameters used in the simulations to acquire the SAR image of the flying helicopter and estimate its motion parameters. The X-band has been widely used for the SAR imaging of static or moving targets. Wideband and ultrawideband SAR systems were used to achieve fine spatial resolution, with a bandwidth of 400 MHz in TELAER [31], 800 MHz in F-SAR [33], and 3.6 GHz in an upgraded PAMIR (in 2012) [25], [30]. In this article, the center frequency is set to f 0 = 10 GHz and the bandwidth is set to B r = 1.25 GHz to achieve a range resolution of 12 cm. Referring to the PAMIR SAR, with the maximum range of 100 km [25], the slant range is set to R 0 = 30 km in this article. The pulsewidth is set to T r = 2.5 μs, compatible to 15 μs in PAMIR SAR [30], leading to a range chirp rate of K r = 500 THz/s. The pulse repetition frequency is set to F a = 2000 Hz, achieving a cross-range resolution of 10 cm, which is comparable to 200-16 000 Hz [31], [32] and 3000 Hz [30], respectively, used for high-resolution X-band SAR. Fig. 2 shows a scale model of Boeing-Sikorsky SB > 1 Defiant attack helicopter [34], [35] used in the simulations. It is propelled with two coaxial main rotors and one tail rotor. Each main rotor is composed of four blades, covering a spin diameter of 15 m [36]. The helicopter has a maximum speed of 463 km/h and can hover at altitude of 1830 m [26]. The main rotor spins at 120-400 r/min (2-6.7 Hz) [37], and the tail rotor spins about three to six times faster than the main rotors [38]. Table II lists eight cases of velocity vectors and spin rates of rotors to be verified and compared in the simulations.  To demonstrate the intermediate results in various stages of the proposed procedure, the helicopter is assumed to move with the default motion parameters specified in case 1 of Table II, and the signals received at platform P 1 , specified with the parameters listed in Table I, are processed with the proposed procedure. In this section, the backscattered signals from fuselage, main rotors, and tail rotor are simulated separately, and their rangecompressed counterparts are computed for further process in the next section.

C. Backscattered Signals From Fuselage
Each radar periodically emits linear frequency modulation (FM) pulse waveform given by where τ is the fast time, T r is the pulsewidth, K r is the range chirp rate, f 0 is the carrier frequency, and is a window function. The helicopter model is decomposed into a bunch of scatterers. The range between P 1 and the nth scatterer, located at (x n , y n , z n ) when η = 0 and moves at velocityv d = (v x , v y , 0), is given by where h 0 = z a , c is the speed of light, is the range at η = 0, is the Doppler centroid frequency, and is the azimuth FM rate. The backscattered signal from the fuselage after demodulation is given by where the fuselage is composed of N 1 scatterers. The signal s rb1 (τ, η) spreads over a wide belt in the τ -η plane, and is first compressed in τ as follows. Take the range Fourier transform of s rb1 (τ, η) to have where c 1n is a constant of integration when the principle of stationary phase technique is applied [39]. Then, multiply a range compression (RC) filter [40] H rc (f τ ) = e jπf 2 with S 11 (f τ , η) to obtain where T a is the aperture time or coherent processing interval. Next, take the inverse Fourier transform of S 12 (f τ , η) in range to have which is the range-compressed signal associated with the fuselage.  Fig. 3(a) shows the schematic of SAR imaging on a spinning point target S(η). During the aperture time, platform P 1 (η) moves in the y-direction at a constant speed V p , at altitude h 0 above the y-axis. The radar look angle is θ and the squint angle is fixed. The rotor mast Q(η) moves at velocity vectorv d , with Q(0) = (x r , y r , z r ). The point scatterer S(η) spins about Q(η), with spin radius a n , spin rate ω r , and initial phase of ψ n at η = 0. The range between P 1 (η) and S(η) is given by

D. Backscattered Signals From Main Rotors
Thus, the backscattered signal from a main rotor (modeled as N 2 scatterers) after demodulation is given by where A 2n is the backscattering coefficient from the nth scatterer. Next, apply RC by first taking the range Fourier transform of s rb2 (τ, η) to have where c 2n is a constant of integration. Then, multiply S 21 (f τ , η) with the RC filter in (4) to obtain followed by an inverse Fourier transformed in range to have which is the range-compressed signal associated with the upper main rotor. The backscattered signal s rb3 (τ, η) from the lower main rotor is range compressed in the same way to have s rc3 (τ, η).

E. Backscattered Signals From the Tail Rotor
shows the schematic of SAR imaging on a point scatterer W (η) spinning about the tail-rotor axle aligned with v d . The tail rotor is centered at rotor mast U (η), with U (0) = (x t , y t , z t ). The point W (η) spins about U (η) with radius b n , spin rate ω t , and initial phase φ n . The range between P 1 (η) and W (η) is given by The backscattered signal from the tail rotor (modeled as N 4 scatterers) after demodulation is given by where A 4n is the backscattering coefficient from the nth scatterer. Next, apply RC by first taking the range Fourier transform of s rb4 (τ, η) to have where c 4n is a constant of integration. Then, multiply the RC filter in (4) with S 41 (f τ , η) to obtain which is the range-compressed signal associated with the tail rotor.

III. VELOCITY VECTOR AND COARSE IMAGE
Different parts of a helicopter operate in different motion patterns, making the backscattered signals overly intricate. Fig. 4 shows an effective procedure to acquire the images and to estimate various motion parameters of a flying helicopter, where the breakdown flowcharts of modified range-Doppler algorithm (RDA) and inverse RDA are shown in Fig. 5. It takes two iterations to estimate all the motion parameters in order to acquire a focused image of the flying helicopter. In the first iteration, the spin rate and initial phase of the main rotors are accurately estimated, but their counterparts of the tail rotor are ambiguous. Thus, a second iteration is taken to prune off residual signals from fuselage and main rotors, leading to more accurate estimation on the spin rate and initial phase of the tail rotor.

A. Estimation of Linear Motion Parameters
Since the backscattered signal from the fuselage is much stronger than those from the main rotors and the tail rotor, we will first estimate the linear motion parameters of the helicopter from the total backscattered signal. The signal attributed to the fuselage is then pruned off, facilitating the estimation on spin parameters of the rotors. The received signal after demodulation is decomposed as which is dominated by the contribution from the fuselage. Fig. 6(a) shows the magnitude of s rb (τ, η) in case 1 of Table II, with v x = v y = 50 m/s. The signal traces are slanted due to nonzero v x . The horizontal streaks near the bottom are attributed to the sinusoidal terms in the range model of rotors.
Next, apply RC by first taking the range Fourier transform of   Table II, v x = v y = 50 m/s. as in (5). The spin motion parameters will be estimated from S 2 (f τ , η). An RFRT [41] can be used to simultaneously compensate range cell migration and range walk effect in which is then inverse Fourier transformed in range to have where cross-coupled terms of τ η that account for range cell migration and range walk in s rc1 (τ, η) have been compensated. Fig. 7 shows the magnitude of s rfrt (τ, η), where most signals are concentrated around τ = 0. Next, apply a PMT to s rfrt (τ, η) as with kernel which is modified from a generalized Radon-Fourier transform (GRFT) [42]. The parameters (α, β) are candidates of (f dc1 , K a1 ), with f dc1 min ≤ α ≤ f dc1 max and K a1 min ≤ β ≤ K a1 max . By substituting (11) into (12), we have Then, the motion parameters are estimated by solving the following optimization problem Fig . 8 shows the magnitude of X pm (α, β) at τ p = 0 s and τ p = 0.8 ns, respectively. The peaks in these two figures are located at (f dc1 ,K a1 ) = (3078.7, 30.65) and (3078.7, 30.37), respectively. The quasi-periodical wavy pattern in Fig. 8(a) is caused by the interference with rotor signals. To avoid such interference, the PMT is applied at multiple τ p 's near τ = 0, to acquire multiple peaks at (α p , β p ). The motion parameters are then estimated from their average as (f dc1 ,K a1 ) = (3, 178.277, 49.767), compared with the true values of (f dc1 , K a1 ) = (3, 177.556, 55.556).

B. Estimation of the Velocity Vector
The velocity vector of the helicopter will be estimated with duo platforms P 1 and P 2 , as specified in Table I. The principle of trigonometry suggests that two observing directions are needed to determine the 2-D velocity vector of a moving target, and a wider angle span between these two observing directions is expected to get more accurate estimation.
The Doppler centroid embedded in the signals received at P 2 is where R 2n0 is the range between P 2 and the nth scatterer at η = 0.
Equations (14) and (15) represent two curves in the v x v y plane, and their intersection is the velocity vectorv d = (v x , v y ). By applying the PMT to s rfrt (τ, η) derived from the backscattered signals at platform P 1 , we obtain an estimatedf dc1 . Similarly, an estimatedf dc2 is obtained by applying the PMT to s rfrt (τ, η) derived from the backscattered signals at platform P 2 . To solve (14) and (15) numerically, a 2-D grid in the v x v y plane is defined, at a uniform spacing of Δv. Then, the difference between the estimatedf dcα and the true f dcα (v d ) at each grid point is computed as with α = 1, 2. A grid search method is then applied to findv d that satisfies leading to an estimated velocity vector ofṽ d = (ṽ x ,ṽ y ). Table III lists the estimated velocity vector in cases 1-4 listed in Table II. The percentage error, 100 × |Δv d |/|v d | (%), with Δv d =ṽ d −v d , is less than 0.03% in all the four cases.

C. Coarse Image of Fuselage
After estimating the velocity vector, we continue the imaging process. First, design a motion compensation (MOCO) filter, using the estimated values off dc1 andK a1 , as which is multiplied with S 2 (f τ , η) in (10) to have where the cross-coupled term of f τ η 2 and the phase term e −j2π(f 0 +f τ )f dc1 η/f 0 accounting for range walk in S 2 (f τ , η) are compensated. Then, take the inverse Fourier transform of S 3 (f τ , η) in range to have Fig. 6(b) shows the magnitude of s moco (τ, η), featuring stripes parallel to the azimuth axis and sinusoid-like wavy pattern caused by spinning rotors. Next, apply azimuth compression (AC) to S 3 (f τ , η) by first taking its azimuth Fourier transform as which is inverse Fourier transformed in range to have where B a1 = K a1 T a is the azimuth bandwidth, and c 2n is a constant of integration. Then, multiply an AC filter [40] H ac (f η ) = e −jπf 2 (17) to have where the second-order phase term in S 5 (τ, f η ) is compensated. Finally, take the inverse Fourier transform of S 6 (τ, f η ) in azimuth to have where the argument of sinc function in η indicates focus in the azimuth direction. Fig. 6(c) shows the magnitude of s 7 (τ, η) after RC, MOCO, and AC, where the fuselage is discernable, while the rotors remain fuzzy. Fig. 9(a) shows the SAR image mapped from s 7 (τ, η) in Fig. 6(c) to an r-y (range-track) plane, with where τ c and η c are the fast time and slow time, respectively, of the image center, which are mapped to r = r 0 = cτ c /2 and y = y 0 = −(V p −ṽ y )η c . Fig. 9(b) shows a scan image obtained by projecting the 3-D scale model of fuselage in Fig. 2 to the r-y plane, viewed from the position of P 1 at η = 0. The scan image looks similar to the SAR image in Fig. 9(a).
Next, apply a threshold mask H th (τ, η) to prune off the fuselage signal from s 7 (τ, η), leaving s 8 (τ, η). An inverse RDA is applied on s 8 (τ, η) to restore an equivalent received signal s rbr (τ, η), which is approximately the backscattering signals from the two main rotors and the tailor rotor, namely where A 2n = c 4n c 3n c 2n c 1n . Fig. 10(a) shows the magnitude of s rbr (τ, η). The horizontal streaks in s rbr (τ, η) become brighter than their counterparts in s rb (τ, η) shown in Fig. 6(a). Then, apply the RC filter in (4) to s rbr (τ, η) to have Fig. 10(b) shows the magnitude of s rcr (τ, η), manifesting an obvious sinusoid-like wavy pattern caused by spinning rotors.

IV. ESTIMATION ON SPIN PARAMETERS OF ROTORS
In the last section, we pruned off most fuselage contribution from the original received signal s rb (τ, η) to obtain s rbr (τ, η), followed by RC to obtain s rcr (τ, η). In this section, a template matching method is applied on s rcr (τ, η) to estimate the spin parameters of the main rotors. Then, the contribution from the main rotors can be restored with these spin parameters and pruned off from the original received signal s rb (τ, η). This marks the completion of the first iteration.
In the second iteration, without the perturbation from the main rotors, the motion parameters of fuselage can be estimated more accurately; hence, the contribution from the fuselage can be pruned off more completely. The residual signal s rcr (τ, η) in the second iteration will be used to estimate the spin parameters of the main rotors and the tail rotor, respectively.

A. Spin Rate of Main Rotors
The two coaxial main rotors spin at the same angular frequency ω r , or spin rate f r = ω r /(2π), but in opposite senses to balance the torque exerted on the fuselage. To detect the spin rate, an angular spectrum is defined for the first time as which is shown in Fig. 11. The spacing between the main peak and the first side peak indicates |f r | = 5.859 Hz.

B. Initial Phase of the Upper Main Rotor
A template matching method is proposed to estimate the initial phase of the main rotors. Each main rotor is composed of four blades, and the tail rotor is composed of eight blades. The blade length L b from its tip to the rotor mast is related to the peak-topeak fast-time difference (Δτ s ) of the sinusoid-like pattern in |s rcr (τ, η)| shown in Fig. 10(b), which is Δτ s 0.192μs. The range R 1n0 from platform P 1 to the blade tip implies that whereτ c 200 μs by inspecting s 7 (τ, η) in Fig. 6(c). Thus, the blade length is estimated asL b 15.106 m, compared to the true value of L b = 15.089 m. The initial phase (at η = 0) of a main rotor will be searched from a grid of {ψ p }, with ψ p = (p − 1)Δψ, 1 ≤ p ≤ N p , and Δψ = (π/2)/(N p − 1). Since four identical blades are symmetrically hinged to the rotor mast, the search interval is reduced to [0, π/2].
The signal backscattered from the upper main rotor, assuming an initial phase ψ p , is labeled as s rbp and is computed in the same way as s rb2 (τ, η) in (7), namely By deducting s rbp (τ, η) from s rbr (τ, η), we obtain s rbp (τ, η). If ψ p is the true initial phase, then the contribution from the upper main rotor will be removed, leading to s rbp (τ, η) s rb3 (τ, η) + s rb4 (τ, η).
Next, apply RC and MOCO to s rbp (τ, η) by first taking the range Fourier transform of s rbp (τ, η) to have Then, multiply it with the RC filter in (4) and the MOCO filter in (16) which is inverse Fourier transformed in range to have If ψ p is the true initial phase, s mocop (τ, η) will contain no contribution from the upper main rotor, which implies the wavy pattern of s moco (τ, η) in Fig. 6(b) becomes weaker, only contributed by the lower main rotor and the tail rotor. In other words, the best matched ψ p is the one that yields the weakest wavy pattern or the most salient features of straight lines in s mocop (τ, η).
The straight-line feature can also be evaluated in terms of a variance defined as is the mean value of θ. Smaller variance implies stronger straight-line feature in |s mocop (τ, η)|. Table IV(a) lists the entropy and variance in matching the initial phase of the upper main rotor. Both the indices suggest that the initial phase is 45 • , with marginal difference from the other candidates. We will check how a second iteration will improve the match.

C. Initial Phase of the Lower Main Rotor
The same template matching method, with the same grid of search, is applied on the lower main rotor. The backscattered signal from the lower main rotor after demodulation, with initial phase ψ p , is labeled as s rbp and is computed as By deducting s rbp (τ, η) from s rb (τ, η), we obtain s rbp (τ, η), where s rb (τ, η) is the residual signal after removing the contribution of the best-matched upper main rotor from s rbr (τ, η). If Next, apply RC and MOCO to s rbp (τ, η) by first taking the range Fourier transform of s rbp (τ, η) to have S 31p (f τ , η) = F τ {s rbp (τ, η)} followed by applying the RC filter in (4) and the MOCO filter in (16) to obtain which is inverse Fourier transformed in range to have Next, apply the Radon transform to s mocop (τ, η) to obtain the entropy and variance, as listed in Table IV(b). Similar to  Table IV(a), both entropy and variance are insensitive to the initial phase, and the minimum values are barely definite. We conjecture the residual signal from fuselage is the cause.
A second iteration following the flowchart in Fig. 4(b) is then conducted to prune off fuselage signals more completely. The best-matched initial phase of each main rotor in the first iteration is used to restore the received signalss rb2 (τ, η) and s rb3 (τ, η) from the main rotors, which are deducted from the original received signal to yield s (2) which is closer to the contribution from fuselage than in the first iteration. A more focused fuselage signal s (2) 7 (τ, η) is acquired by applying the flowchart in Fig. 4(b) to process s (2) rb (τ, η). Next, apply a threshold mask H (2) th (τ, η) to retrieve the fuselage signals s (2) 8 (τ, η) from s (2) 7 (τ, η). The inverse RDA is then applied to s (2) 8 (τ, η) to restore the received signals from the fuselage,s (2) rb1 (τ, η). By deductings (2) rb1 (τ, η) from the original received signal, we have signals dominated by the rotors from which more accurate spin parameters of rotors can be estimated.

D. Enhancement of Tail Rotor Signals
Next, apply the previous template matching method on s (2) rbr (τ, η) to obtain the entropy and variance in matching the upper main rotor. The results are listed in Table IV(c), which become much more sensitive to the initial phase than their counterparts in Table IV(a). After matching the upper main rotor, its backscattered signals are pruned off from s (2) rbr (τ, η), to match the lower main rotor. Table IV(d) lists the entropy and variance in matching the lower main rotor, which are much more sensitive to the variation of initial phase compared to their counterparts listed in Table IV(b). By removing the backscattered signals from both main rotors, s (2) rb2 (τ, η) ands (2) rb3 (τ, η), we have s (2) rb4 (τ, η) = s (2) rbr (τ, η) −s (2) rb2 (τ, η) −s (2) rb3 (τ, η) s rb4 (τ, η) which is mainly contributed by the tail rotor.

E. Spin Rate of the Tail Rotor
Next, s (2) rb4 (τ, η) is range compressed to have The spin rate f t of the tail rotor can be estimated by defining an angular spectrum on s (2) rc4 (τ, η), similar to (20), as rc4 (τ, η)|} . Fig. 12(a) and (b) show the angular spectrum of the tail rotor in the first and second iterations, respectively. The angular spectrum in the first iteration is too fuzzy to estimate f t . After pruning  TABLE V  INDICES FOR MATCHING THE INITIAL PHASE OF THE TAIL ROTOR  IN ITERATION  off signals from fuselage and main rotors more completely in the second iteration, the magnitude of the spin rate is estimated as f t 32.223 Hz by observing the difference between the central spike and a side spike on the second sidelobe in Fig. 12(b).

F. Initial Phase of the Tail Rotor
Next, template matching is applied on the tail rotor signal to estimate its initial phase φ p with a grid of search as {φ p }. The backscattered signal from the tail rotor after demodulation, assuming initial phase of φ p at η = 0, is labeled as s rbp (τ, η) and is computed in the same way as s rb4 (τ, η) in (8), namely By deducting s rbp (τ, η) from s (2) rb4 (τ, η), we obtain s rbp (τ, η). If φ p is the true initial phase, then The condition of (22) means no residual signal at all. Table V shows that the entropy and variance of s mocop (τ, η), by processing s rbp (τ, η), have opposite trend compared to those for main rotors listed in Table IV. It is conjectured that in s (2) rb4 (τ, η), the residual signals accrued to the fuselage, and main rotors still outweigh the signals from the tail rotor in computing entropy or variance. Fig. 13 shows that the magnitude of s rbp (τ, η) with correct initial phase of φ p = 39 • is significantly smaller than that with an incorrect guess of φ p = 54 • , which confirms the null condition in (22). Thus, we define an alternative index of total energy as where E 0 = 10 −8 is a normalization constant for this case, and s rbp [m, n] are samples of s rbp (τ, η), with 1 ≤ m ≤ N r , 1 ≤ n ≤ N a . Table V shows that the total energy is an effective index to determine the initial phase of the tail rotor as φ p = 39 • . The corresponding received signals attributed to the tail rotor are restored ass (2) rb4 (τ, η). By implementing a third-order and higher order iteration in the same way as the second iteration, the sensitivity to the initial phase is improved, and the estimated initial phase remains the same.

V. IMAGING RESULTS AND DISCUSSIONS
After two iterations described in the last section, the velocity vector of the helicopter and the spin parameters of the rotors are accurately estimated. The received signals from fuselage and rotors, moving at the estimated velocity without spinning, can be restored and processed with the modified RDA depicted in the flowchart of Fig. 5(a) to acquire an SAR image at each platform. Fig. 14(a) and (b) shows the acquired SAR images |s 7 (τ, η)| and the corresponding scan images, respectively, viewed from P 1 . The shape of helicopter is recognizable in Fig. 14(a). The structural parts around the tail can be clearly identified, including tail axle, active rudder, and horizontal stabilizer. The backscattered signal from the two main rotors is weaker than that from the fuselage, and that from the tail rotor is weaker than that from the main rotors. The images of fuselage and rotors can be separately acquired, as shown in Fig. 14(c) and (e), respectively, which are discernible as compared with their counterpart scan images in Fig. 14(d) and (f), respectively. Fig. 15(a)-(d) shows the angular spectrum of s rcr (τ, η) in (19) for cases 5-8, respectively. The spacing between the main spike and the first side spike indicates the spin rate |f r | of the main rotors, which is |f r | = 4.883 Hz for cases 5 and 6, and |f r | = 6.836 Hz for cases 7 and 8, matching well with the default parameters listed in Table II. Fig. 16(a) shows the angular spectrum of s (2) rc4 (τ, η) in (21) for the tail rotor, in cases 5 and 7, after the second iteration. Similarly, Fig. 16(b) shows the angular spectrum of the tail rotor in cases 6 and 8, after the second iteration. The highest side spike looks ambiguous in all the four cases, possibly smeared by the residual fuselage signals even after two iterations.

A. Variation of Spin Rates
To confirm this speculation, Fig. 17(a) and (b) shows the counterpart angular spectrum of the tail rotor, in the absence of fuselage signals. The highest side spikes become obvious, and the spin rate of tail rotor can be accurately estimated. Although the motion parameters of tail rotor cannot always be estimated, the helicopter can still be identified from its fuselage and main rotors.

B. Effects of Noise
Up to this point, the simulations are used to verify the efficacy of the proposed procedure in estimating the motion parameters and acquiring the target image, in the absence of noise. Next,  the proposed procedure will be tested by inserting additive white Gaussian noise at a specific signal-to-noise ratio (SNR) level to the received signal s rb (τ, η) in (9).  Table VI summarizes the estimated parametersf dc1 andK a1 in the presence of noise. The peak position can still be determined in Fig. 18(d) at SNR = −15 dB, but becomes difficult in Fig. 18(c) at SNR = −20 dB. In the latter, the motion parametersf dc1 andK a1 are not accurately    is observed; thus, the spin rate of the main rotors cannot be estimated at all. In summary, the translational motion parameters can be accurately estimated, and high-quality SAR images of the fuselage can be acquired at SNR ≥ −15 dB. The estimation of the spin rate is more sensitive to noise, because the backscattered signals from the rotors are much weaker than those from the fuselage.
Next, the SNR is gradually increased from −15 dB to observe the effects of noise on the angular spectrum. Fig. 19(a) and (b) shows the magnitude of the received signal s rb (τ, η) and s 7 (τ, η) of fuselage, respectively, at SNR = 0 dB. Fig. 19(c) shows the  angular spectrum of s rcr (τ, η), with a ditch marked by red arrow. The spacing between the ditch and the main spike is used to estimate the spin rate |f r | of the main rotors. Table VII lists the indices for matching the initial phase of the upper main rotor in the first iteration, at SNR = 0 dB and 5 dB, respectively. It is observed that the entropy and the variance drop to minima at the correct initial phase ψ p = 45 • , indicating that the template matching is effective at SNR ≥ 0 dB. The final SAR image of the whole flying helicopter is shown in Fig. 19(d).

C. Comparisons With SOTA Algorithms
To further verify the efficacy of the proposed procedure, three SOTA algorithms are implemented for comparison. Fig. 20(a)-(d) shows the SAR images of the flying helicopter in case 1, acquired with the proposed procedure, conventional RDA, chirp scaling algorithm (CSA), and omega-K algorithm (omega-KA), respectively [40]. With the proposed procedure, the velocity vector of the flying helicopter is accurately estimated, the shape of helicopter is well focused, and the (τ , η) position of the target center is very close to the true value of (200 μs, 0 s).
In contrast, the fuselage shapes in Fig. 20(b)-(d) are discernible but too fuzzy for identifying the helicopter model. The (τ , η) position of the target is more or less deviated from the true one because the motion parameters are not accurately estimated.
To quantify the similarity between the acquired image with our procedure and those with the SOTA algorithms, a structural similarity (SSIM) index, a mean squared error (MSE), and a peak signal-to-noise ratio (PSNR) are adopted. The SSIM index between imagesx andȳ is defined as [44] SSIM are specified for comparing luminance, contrast, and structure, respectively; μ η and σ η are the mean and standard deviation, respectively, of imageη, with η = x, y, and σ xy is the covariance betweenx andȳ. The SSIM index lies between [0, 1], and is equal to 1 if two images are identical.
The MSE between imagesx andȳ, of size N r × N a , is defined as [45] MSE (x,ȳ) = 1 Lower MSE implies that both images are more similar.  The PSNR is the ratio between the maximum signal power and the noise power as [45] PSNR (x,ȳ) = 10 log 10 where the maximum signal power refers to the largest pixel intensity of an image and the noise power refers to the MSE between the two images. Higher PSNR indicates imagex has higher pixel intensity and is more similar to the reference imageȳ. Table VIII lists the SSIM indices, MSE, and PSNR between the image in Fig. 20(a) acquired with the proposed procedure and those in Fig. 20(b)-(d) acquired with the SOTA algorithms. The image with conventional RDA is most similar to that of the proposed procedure, judging from its highest SSIM index, lowest MSE, and highest PSNR. Among the four images in Fig. 20, the one acquired with the proposed procedure manifests the sharpest contour of the helicopter. The SSIM indices between Fig. 20(a) and Fig. 20(b)-(d) are lower than 0.9, suggesting that the images acquired by the SOTA algorithms may have lower quality than ours. Table IX lists the computational load of the constituent algorithms in the flowchart of Fig. 5, where N r and N a are the numbers of range bins and azimuth bins, respectively, N α and N β are the numbers of candidates for α and β, respectively, of the kernel K pm (α, β, η) in (13), and FT, IFT, and FFT are the abbreviations of Fourier transform, inverse Fourier transform, and fast Fourier transform, respectively. Table X lists the computational load of the constituent algorithms used for estimating the spin rate and template matching in the flowchart of Fig. 4. The overall computational load is proportional to O (N r N a log 2 (N r N a )).

D. Computational Load Analysis
The proposed procedure has the advantages of acquiring SAR images without knowing the motion parameters of the target. Instead, the target velocity is estimated by using the two Doppler centroid frequencies derived from the received signals at the The extra computational load required to estimate the target velocity is relatively small. As listed in Tables IX and X, the computational loads of the proposed procedure and the SOTA methods are compatible. In terms of CPU (elapse) time, the simulation results presented in this article are computed with MATLAB R2021b in Windows 10, on a PC desktop with Intel Core i7-11700 eight-core 2.5-GHz processor and 64-GB RAM. To acquire one SAR image of the fuselage, the proposed procedure takes CPU time of 29.74 s, while the conventional RDA takes 31.23 s, the CSA takes 38.98 s, and the omega-KA takes 452.31 s, in which the Stolt mapping takes 420.86 s. The proposed procedure takes 9.80 s to estimate the motion parameters, Doppler centroid frequency, and azimuth FM rate, which is about one-third of the total CPU time.

E. Highlight of Contributions
The novelties and contributions in this article are summarized as follows.
1) This is the first article to propose a complete and rigorous procedure to acquire the SAR images of a swiftly moving helicopter, which bears both the translational and spinning motions.
2) The signals backscattered from rotating parts are more difficult to untangle than those from linear moving parts. The images of rotating target, such as airplane can be focused with ISAR if its angular velocity is low (≤ 1 rad/s) [20], but rotors do not belong in this category. The micro-Doppler signatures and ISAR imaging of rotating helicopter blades were discussed [5], [6], without considering the blockage of fuselage and its translational motion. This is the first paper to acquire the SAR image of a whole moving helicopter, including fuselage with translational motion, and spinning rotors with both the translational and rotational motions.
3) The two-stage template matching method is proposed for the first time to meticulously separate backscattered signals from moving fuselage and spinning main rotors and tail rotor, respectively.
4) The PMT modified from the GRFT is integrated with the RFRT for the first time so that the phase terms related to Doppler centroid frequency and azimuth FM rate can be compensated in the range frequency-azimuth domain. 5) This is the first article to propose a duo airborne SAR for estimating the target velocity by deriving two Doppler centroid frequencies from the received signals. 6) An angular spectrum is defined for the first time to estimate the spin rate of rotors from the received signals. A feasible procedure is also proposed to meticulously separate the signals from fuselage, main rotors, and tail rotor to improve the estimation accuracy of the tail-rotor spin rate. 7) Scan image projected from a 3-D target model is proposed for the first time to facilitate the verification of acquired SAR images and identification of the target. 8) The weak signals from spinning rotors and strong signals from the fuselage are intertwined, and a feasible procedure to untwine these signals for better SAR imaging is proposed for the first time. The acquired images are verified and examined under various scenarios. Previous works in the literature rarely discussed the imaging of tail rotor under the influence of fuselage and main rotors. The technical difficulties in determining the tail-rotor spin rate have also been resolved.

VI. CONCLUSION
In this article, a complete and rigorous procedure was developed to estimate the translational and spinning motion parameters of a flying helicopter, as well as to acquire SAR images for identification. A modified RDA was proposed to acquire the fuselage image first. With duo airborne SAR platforms, the velocity vector of the helicopter can be accurately estimated from the received signals. Then, a two-stage template matching method was proposed to accurately estimate the spin rates of main rotors and tail rotor, as well as their initial phases.
The performance of the proposed procedure was verified under various combinations of velocity vector and spin rate of flying helicopter. The effects of noise, comparison with three SOTA algorithms, and the computational load were also elaborated to verify the efficacy of the proposed procedure.