Fast Solution of Scattering and Micro-Doppler Features from Moving Target Using a Tailored Shooting and Bouncing Ray Method

: In this paper, we present a tailored shooting and bouncing ray (SBR) method for the fast solution of electromagnetic (EM) scattering from a moving target. And, the micro-Doppler features of the moving target are investigated using a time-frequency analysis technique. In our method, a dynamic spatial division technique is employed to accelerate facet information processing and ray-tracing progress of the moving target. At ﬁrst, the two coordinate systems are established, which are the geodetic coordinate system (GCS) and the local coordinate system (LCS). In GCS, the target is moving with translation and rotation. The dynamic spatial division is established in LCS to store the facet information and remain relatively stationary to the target. In comparison with the traditional SBR method, this technique avoids repetitive spatial division at each moment in the GCS. Then, ray tracing is performed to ﬁnd the illuminated facets in the LCS. Finally, the scattering ﬁeld and the phase compensation are computed in the GCS. In numerical simulations, the veriﬁcation and computation efﬁciency comparison are provided using our method and other solutions (MLFMM, RL-GO, and traditional SBR). Moreover, the micro-Doppler features are extracted and analyzed using the time-frequency analysis technique, which includes the precession and spin of the missile, and the rotation of the aircraft. Meanwhile, the micro-Doppler spectra of the target is also compared with the theoretical Doppler of equivalent strong scattering points, which are obtained using the instantaneous high-resolution range proﬁle (HRRP).


Introduction
Research on the electromagnetic (EM) scattering of a target is essential for radar detection and target recognition.With the development of computational electromagnetics, the methods for solving EM scattering problems have matured considerably over the past decades.For example, the method of moments (MOM) [1], the finite difference time domain (FDTD) [2], and the multilevel fast multipole method (MLFMM) [3] can achieve highprecision calculation; however, these methods are incapable for electrically large structures due to computing resources.Therefore, some high-frequency approximation methods are becoming increasingly popular in the field of EM calculation due to faster calculation speed and less resource consumption.Among them, physical optics (PO), geometric optics (GO), and ray-launching geometrical optics (RL-GO) [4] are more common, as well as the shooting and bouncing ray (SBR) method [5], which combines the two methods of PO and GO.In the SBR method, GO regulates the rules of wave propagation, and PO integrations are implemented to calculate the scattering field.Therefore, the SBR method is often applied to solve the EM scattering of electrically large complex structural targets [6].In addition, the physical diffraction theory (PTD) combined with SBR can solve the problem of edge diffraction [7] to improve the calculation accuracy.A series of methods have been proposed to accelerate the ray-tracing process in SBR.The tree structure algorithm [8] and some of its improvements [9,10] were proposed to speed up the process of ray tracing.The time-domain SBR based on the beam tracing (BT) technique [11] was also proposed to improve the accuracy and computational efficiency.The Open Graphics Library (OpenGL) hybrid SBR method [12] was proposed to accelerate shadow removal.At the same time, some technical solutions [13] combined with SBR have been used to process complex scenes with remarkable effectiveness.
The methods mentioned above were focused on solving EM scattering of a stationary target.There are few works on researching the EM scattering and micro-Doppler features of a moving target based on high-frequency algorithms.The research of frequency characteristics based on EM scattering helps to provide theoretical support for radar detection and target identification [14,15], as well as SAR image processing [16,17].Since the micro-Doppler features generated using different micro-motions are different, they also can be utilized to classify and recognize targets [18][19][20][21].Chen has conducted a lot of excellent work on the micro-Doppler.He introduced the micro-Doppler phenomenon from microwave radar and established the micro-Doppler models caused by different motion forms [22].With the deepening of research, the time-frequency analysis technology has become the main method to study micro-Doppler features [23].Gao et al. [24] adopted this technology to analyze the micro-Doppler features of ballistic targets.The micro-Doppler features of rotor blades in different frequency bands were also investigated by [25,26].This joint time-domain and frequency-domain analysis technology [27] can reveal the distribution of the Doppler spectra in the time domain.
However, all the works [20][21][22][23][24] above were based on the point scattering model of a moving target.This method just considers the scattering echo of a few strong scattering points on the target, ignoring the influence of structural features on the scattering echo.And the traditional EM scattering algorithms, such as MLFMM and MOM [26], consume too much computational resources in dealing with EM scattering from a moving target.In this paper, the micro-Doppler features of a moving target are investigated using EM modeling, and the dynamic spatial division technique-based SBR method is adopted to accelerate the solution of EM scattering with a moving target.At first, the geodetic coordinate system (GCS) and the local coordinate system (LCS) of the target at each moment are established.In LCS, the boundary of the target at the initial moment is determined and the spatial bounding box is constructed.The spatial bounding box stores facet information and remains relatively stationary to the moving target.In comparison with the traditional SBR method, this process avoids the spatial division at each moment.This will save lots of computational resources.Then, ray tracing is performed to find the illuminated facets in LCS.The illuminated facets information is converted from the LCS to the GCS by the rotation matrix-based coordinate system transformation.Finally, the scattering field and the phase compensation are computed in the GCS.The micro-Doppler features are extracted using the short-time Fourier transform (STFT), which is one of the time-frequency analysis techniques.In the simulation section, the effectiveness of the proposed method is verified by comparing the backscattering results of our method and other solutions (MLFMM, RL-GO, and traditional SBR), and the computation efficiency comparison is also given.Meanwhile, the micro-Doppler features of the moving target are analyzed, which includes the precession and spin of the missile, and the rotation of the aircraft.Moreover, the micro-Doppler spectra of the target is also compared with the theoretical Doppler of equivalent strong scattering points, which are obtained by the instantaneous high-resolution range profile (HRRP).
The remainder of this paper is organized as follows: In Section 2, the fast solution of the EM scattering model for a moving target is established, and the micro-Doppler features of the moving target are extracted using STFT.In Section 3, the accuracy and efficiency are verified by comparison with the results of our method and other solutions (MLFMM, RL-GO, and traditional SBR).The simulation results are discussed, including the micro-Doppler spectra of the missile and aircraft.The conclusion is reported in Section 4.

Motion Modeling of the Target
Motion modeling is the foundation of solving the EM scattering of the moving target.Since the interaction time between the radar signal and targets is too short, target motion can be neglected.Therefore, the quasi-static method is utilized for motion modeling.Usually, the motion of the target can be divided into two parts: one is the translation relative to the radar irradiation direction, and the other is the rotation of the target relative to the radar coordinate system, as seen in Figure 1.
Remote Sens. 2023, 15, x FOR PEER REVIEW 3 of 16 the micro-Doppler spectra of the missile and aircraft.The conclusion is reported in Section IV.

Motion Modeling of the Target
Motion modeling is the foundation of solving the EM scattering of the moving target.Since the interaction time between the radar signal and targets is too short, target motion can be neglected.Therefore, the quasi-static method is utilized for motion modeling.Usually, the motion of the target can be divided into two parts: one is the translation relative to the radar irradiation direction, and the other is the rotation of the target relative to the radar coordinate system, as seen in Figure 1.In a three-dimensional space, the position vector of point  in   is    = [ _ ,  _ ,  _ ]  .At moment   , the position vector  _  in   is obtained by rotating    for an angle  around the rotation axis  = [  ,   ,   ]  .It can be written as where ℜ|  () is the rotation matrix from   to   at moment  , which can be expressed as where,  3×3 is the unit matrix;  ̃ is a skew symmetric matrix of ;  is the rotation angle of the target around the rotation axis  at moment   .Considering that the target also has a translation from   to   , then the motion of the target can be expressed in where  0 |   0 is the distance vector from   to   of the target at initial moment  0 ;   is the translational vector from the position of moment  0 to   in   .In this paper, the superscripts and subscripts of symbols such as (⋅)|   represent the transformation relationship from the coordinate system LCS to GCS. 0 and () represent the LCS of the target at the initial moment and the nth moment, respectively.In a three-dimensional space, the position vector of point p in XYZ LCS is r L p = r p_x , r p_y , r p_z ] T .At moment t n , the position vector r G p_n in xyz GCS is obtained by rotating r L p for an angle θ around the rotation axis ω = ω x , ω y , ω z ] T .It can be written as where R|

L(n) G
is the rotation matrix from XYZ LCS to xyz GCS at moment t n , which can be expressed as where, I 3×3 is the unit matrix; ∼ ω is a skew symmetric matrix of ω; θ is the rotation angle of the target around the rotation axis ω at moment t n .Considering that the target also has a translation from XYZ LCS to xyz GCS , then the motion of the target can be expressed in where R 0 | G 0 G is the distance vector from xyz GCS to UVW GCS of the target at initial moment t 0 ; vt n is the translational vector from the position of moment t 0 to t n in XYZ LCS .In this paper, the superscripts and subscripts of symbols such as (•)| L G represent the transformation relationship from the coordinate system LCS to GCS.L0 and L(n) represent the LCS of the target at the initial moment and the nth moment, respectively.

Fast Solution for EM Scattering of the Moving Target Using the Tailored SBR Method
In the SBR method, GO regulates the rules of wave propagation, and then PO solves the induced EM current of the illuminated triangular facets.For example, the ray is reflected three times between the different facets of the trihedral shown in Figure 2. The initial facet illuminated by the incident ray is T 0 and the facet illuminated by the ray after the first-order reflection is T 1 ; then, the facet illuminated by the second-order reflection is T 2 .The total EM scattering field is contributed to by these illuminated facets, which include first-order and higher-order scattering.So, the SBR method greatly improves the calculation accuracy compared to the PO method.

Fast Solution for EM Scattering of the Moving Target using the Tailored SBR Method
In the SBR method, GO regulates the rules of wave propagation, and then PO solves the induced EM current of the illuminated triangular facets.For example, the ray is reflected three times between the different facets of the trihedral shown in Figure 2. The initial facet illuminated by the incident ray is  0 and the facet illuminated by the ray after the first-order reflection is  1 ; then, the facet illuminated by the second-order reflection is  2 .The total EM scattering field is contributed to by these illuminated facets, which include first-order and higher-order scattering.So, the SBR method greatly improves the calculation accuracy compared to the PO method.In the traditional SBR method, the consumption of computational resources is mainly focused on two parts, creating spatial division to process facet information and ray tracing to find illuminated facets.In addition, the spatial division needs to be rebuilt for each moment in   , which will consume large computational resources.In this paper, the rotation matrix-based coordinate system transformation is employed to establish the dynamic spatial division, and bidirectional ray-tracing technology [28] is adopted to accelerate the faceted information process and the ray-tracing process.
As shown in Figure 2b, in   , the target is moving with translation and rotation.The dynamic spatial division is established in   to store facet information and remain relatively stationary to the target, which avoids the spatial division at each moment in   .The transformation relationship from the coordinate system   to   is given by Equation (3).
Figure 3 represents the flowchart for the dynamic spatial division.In this process, the spatial bounding box is moving with the target and remains relatively stationary.Thus, the facet number information stored in the leaf nodes does not change with the movement of the target.Only the coordinate information corresponding to the leaf nodes needs to be updated.Then, the new round of spatial division can be completed.Compared to traditional spatial division, there is no need to recognize the boundaries of the target and create a spatial bounding box at each moment.This will save lots of computational resources.

Moment t0
Moment tn In the traditional SBR method, the consumption of computational resources is mainly focused on two parts, creating spatial division to process facet information and ray tracing to find illuminated facets.In addition, the spatial division needs to be rebuilt for each moment in xyz GCS , which will consume large computational resources.In this paper, the rotation matrix-based coordinate system transformation is employed to establish the dynamic spatial division, and bidirectional ray-tracing technology [28] is adopted to accelerate the faceted information process and the ray-tracing process.
As shown in Figure 2b, in xyz GCS , the target is moving with translation and rotation.The dynamic spatial division is established in XYZ LCS to store facet information and remain relatively stationary to the target, which avoids the spatial division at each moment in xyz GCS .The transformation relationship from the coordinate system XYZ LCS to xyz GCS is given by Equation (3).
Figure 3 represents the flowchart for the dynamic spatial division.In this process, the spatial bounding box is moving with the target and remains relatively stationary.Thus, the facet number information stored in the leaf nodes does not change with the movement of the target.Only the coordinate information corresponding to the leaf nodes needs to be updated.Then, the new round of spatial division can be completed.Compared to traditional spatial division, there is no need to recognize the boundaries of the target and create a spatial bounding box at each moment.This will save lots of computational resources.
where êTE | are direction vectors of TE and TM, respectively, which can be obtained by where ki is the direction vectors of the incident wave; N The induced electric and magnetic current of illuminated facet Tri_t n can be solved using the PO method where is the reflection coefficient of the illuminated facet Tri_t n , which can be obtained by in Equation ( 9), Z TE = −µ r and Z TN = ε r , which are the relative permittivity and relative permeability of Tri_t n .
Remote Sens. 2023, 15, x FOR PEER REVIEW 6 of 16 in Equation ( 9),   = −  and   =   , which are the relative permittivity and relative permeability of _  .The far-field scattering of _  in   can be approximated using the PO integral The far-field scattering of Tri_t n in xyz GCS can be approximated using the PO integral where η = ε 0 /µ 0 ; k = ω √ ε 0 µ 0 ; î and ŝ are the unit vectors in the incident direction and the scattering direction, respectively; R is the distance from the source point to the observed point.Figure 5 represents the flowchart for solving the scattering field of Tri_t n .The χ is the Gordon surface integral [29], and it can be written as where S represents the area of triangular facet are the vector along the mth edge and the midpoint coordinate of the facet can be obtained by Equation (3).In Equation ( 12), ĥ can be obtained by ĥ (13) considering that the facet also has translational motion, and the translational vector is T = vt n , as shown in Figure 3. Hence, the phase compensation ϕ T for translation is performed then, the total phase function of the illuminated triangular facet moving from Tri_t 0 to Tri_t n can be expressed as Remote Sens. 2023, 15, x FOR PEER REVIEW 7 of 16 then, the total phase function of the illuminated triangular facet moving from _ 0 to _  can be expressed as Figure 5. Flowchart for solving the scattering field of .
In the SBR method, high-order ray tracing is also a critical step.The irradiation area of the reflected ray usually covers multiple triangular facets, as shown in Figure 6.In this paper, the bidirectional ray-tracing technology is adopted to judge whether the facets are illuminated.The facet 0 is irradiated by the incident wave.The reflected ray starts from the center of facet 0 and reaches facet 13.Facet 13 is identified using forward ray tracing as one of the illuminated facets.The backward ray tracing is performed in the opposite direction of the reflected rays, which starts from the neighboring centers of facet 13 (facets 7, 12, and 14).If the backward ray intersects facet 0, then the neighboring facet is also one In the SBR method, high-order ray tracing is also a critical step.The irradiation area of the reflected ray usually covers multiple triangular facets, as shown in Figure 6.In this paper, the bidirectional ray-tracing technology is adopted to judge whether the facets are illuminated.The facet 0 is irradiated by the incident wave.The reflected ray starts from the center of facet 0 and reaches facet 13.Facet 13 is identified using forward ray tracing as one of the illuminated facets.The backward ray tracing is performed in the opposite direction of the reflected rays, which starts from the neighboring centers of facet 13 (facets 7, 12, and 14).If the backward ray intersects facet 0, then the neighboring facet is also one of the facets irradiated by the reflected wave.For a facet of the illuminated area, its adjacent facets will be checked by backward ray tracing.This bidirectional ray-tracing technique combined with neighboring findings improves the efficiency of ray tracing.For the illuminated facets, the far-field scattering can be approximated by Equation (10).In the SBR method, high-order ray tracing is also a critical step.The irradiation area of the reflected ray usually covers multiple triangular facets, as shown in Figure 6.In this paper, the bidirectional ray-tracing technology is adopted to judge whether the facets are illuminated.The facet 0 is irradiated by the incident wave.The reflected ray starts from the center of facet 0 and reaches facet 13.Facet 13 is identified using forward ray tracing as one of the illuminated facets.The backward ray tracing is performed in the opposite direction of the reflected rays, which starts from the neighboring centers of facet 13 (facets 7, 12, and 14).If the backward ray intersects facet 0, then the neighboring facet is also one of the facets irradiated by the reflected wave.For a facet of the illuminated area, its adjacent facets will be checked by backward ray tracing.This bidirectional ray-tracing technique combined with neighboring findings improves the efficiency of ray tracing.For the illuminated facets, the far-field scattering can be approximated by Equation (10).

Time-Frequency Analysis for the Micro-Doppler
In Equations ( 10) and ( 15), the scattering field and phase function of the illuminated triangular facet moving with rotation and translation are derived.This can be simplified and written as

Time-Frequency Analysis for the Micro-Doppler
In Equations ( 10) and ( 15), the scattering field and phase function of the illuminated triangular facet moving with rotation and translation are derived.This can be simplified and written as where ŝ is the unit vector in the scattering direction; E s n represents the amplitudes of the scattering field.The change process of phase Φ n and micro-Doppler frequency shift induced by the motion of this facet to the incident direction can be obtained by where f is the incident frequency; v is the velocity vector of translation; R| G is the coordinate system transformation from the LCS to the GCS induced by rotation; and the amplitude of the scattered field corresponding to this Doppler frequency f m−D n is E s n .In order to solve the overall micro-Doppler shift of the target with motion, it is often necessary to resort to the STFT, which is one of the time-frequency analysis techniques.The steps of the STFT are performed as follows.

•
Step 1: The window function w(t) is moved to the starting position of the time sequence E s | n 0 .At this time, the center position of w(t) is at t = τ 0 , and the scattering field sequence is processed by the window function • Step 2: The fast Fourier transform (FFT) is performed on a small sequence covered by the window function to obtain the Doppler shift of this sequence and it can be expressed as where F[•] represents the FFT operation; • Step 3: After completing the FFT of the first sequence, the center of the window function is moved to t = τ 1 .The moving distance should be less than the width of the window function, to ensure that there is an overlap between the two windows; • Step 4: Repeat the above operation as shown in Figure 7, continuously slide the window function and process the sequence using FTT, and finally get the Doppler frequency spectra of all segments on τ 0 ∼ τ m .
In order to solve the overall micro-Doppler shift of the target with motion, it is often necessary to resort to the STFT, which is one of the time-frequency analysis techniques.The steps of the STFT are performed as follows. • Step 1: The window function () is moved to the starting position of the time sequence   | 0  .At this time, the center position of () is at  =  0 , and the scattering field sequence is processed by the window function • Step 2: The fast Fourier transform (FFT) is performed on a small sequence covered by the window function to obtain the Doppler shift of this sequence and it can be expressed as where [⋅]represents the FFT operation; • Step 3: After completing the FFT of the first sequence, the center of the window function is moved to  =  1 .The moving distance should be less than the width of the window function, to ensure that there is an overlap between the two windows; Step 4: Repeat the above operation as shown in Figure 7, continuously slide the window function and process the sequence using FTT, and finally get the Doppler frequency spectra of all segments on  0 ~ .Thus, the micro-Doppler features of a moving target based on EM scattering are extracted using the time-frequency analysis technique.This can be expressed as ( ) Thus, the micro-Doppler features of a moving target based on EM scattering are extracted using the time-frequency analysis technique.This can be expressed as where i = 0 ∼ m, S( f m−D , τ i ) is the micro-Doppler distribution at moment τ i .
For better understanding, Figure 8 represents the flowchart for the method proposed in this paper to accelerate the solution of EM scattering and the micro-Doppler time-frequency spectra with a moving target.
For better understanding, Figure 8 represents the flowchart for the method proposed in this paper to accelerate the solution of EM scattering and the micro-Doppler time-frequency spectra with a moving target.

Simulation Results and Discussion
In this section, first, the EM scattering of the moving target is simulated using our method and compared with the results of MLFMM, RL-GO, and the traditional SBR method [28].Meanwhile, the computation efficiency comparison is given.Then, the micro-Doppler features of the moving target are analyzed, which include the precession and spin of the missile and the rotation of the aircraft.Moreover, the micro-Doppler spectra of the target is also compared with the theoretical Doppler of the equivalent strong scattering points, which are obtained using the instantaneous HRRP.The geometric model adopted in this paper is shown in Figure 9.All simulations were performed on a computer with Intel(R) Xeon(R) Platinum 8275CL at 3.6 GHz and 384 GB RAM.

Simulation Results and Discussion
In this section, first, the EM scattering of the moving target is simulated using our method and compared with the results of MLFMM, RL-GO, and the traditional SBR method [28].Meanwhile, the computation efficiency comparison is given.Then, the micro-Doppler features of the moving target are analyzed, which include the precession and spin of the missile and the rotation of the aircraft.Moreover, the micro-Doppler spectra of the target is also compared with the theoretical Doppler of the equivalent strong scattering points, which are obtained using the instantaneous HRRP.The geometric model adopted in this paper is shown in Figure 9.All simulations were performed on a computer with Intel(R) Xeon(R) Platinum 8275CL at 3.6 GHz and 384 GB RAM.

EM Scattering of the Moving Target
In Figure 10, the backscattering of the dihedral angle and combined model are compared with our method and other solutions (MLFMM, RL-GO, and the traditional SBR). Figure 10a is the scattering of the dihedral angle composed of two 0.5 m × 0.5 m plates with HH polarization.The dihedral angle rotates around the  -axis with an angular velocity of  = /18 rad/s.Where the incident frequency is  = 10 GHz, the incident angle is   = 0°, and the azimuth angle is   = 0°.Figure 10b is the scattering of the combined model composed of two cubes with side lengths of 0.4 m and 0.2 m, respectively.The combined model rotating around the y -axis with an angular velocity is  = /2 rad/s.
Where the incident frequency is  = 8 GHz, the incident angle is   = 90°, the azimuth angle is   = 0°, and the polarization is VV.In this section, since there is need to evaluate the effectiveness of our method, some statistical concepts are utilized to measure the deviation degree of our method from other solutions.The mean absolute error  and its standard deviation  [7] are defined as where   are the differences between the predicted RCS values of our method and other solutions; and N is the total number of moments for target movement.Take the MLFMM and RL-GO results as the reference values, the  and  for our method are presented in Table 1.It is obvious that the  and  of our method are about 0.36~0.88 and 1.67~2.81.It indicates that these results using different methods match very well.

EM Scattering of the Moving Target
In Figure 10, the backscattering of the dihedral angle and combined model are compared with our method and other solutions (MLFMM, RL-GO, and the traditional SBR). Figure 10a is the scattering of the dihedral angle composed of two 0.5 m × 0.5 m plates with HH polarization.The dihedral angle rotates around the y-axis with an angular velocity of ω = π/18 rad/s.Where the incident frequency is f = 10 GHz, the incident angle is θ i = 0 • , and the azimuth angle is φ i = 0 • .Figure 10b is the scattering of the combined model composed of two cubes with side lengths of 0.4 m and 0.2 m, respectively.The combined model rotating around the y-axis with an angular velocity is ω = π/2 rad/s.Where the incident frequency is f = 8 GHz, the incident angle is θ i = 90 • , the azimuth angle is φ i = 0 • , and the polarization is VV.In this section, since there is need to evaluate the effectiveness of our method, some statistical concepts are utilized to measure the deviation degree of our method from other solutions.The mean absolute error µ and its standard deviation σ [7] are defined as where x i are the differences between the predicted RCS values of our method and other solutions; and N is the total number of moments for target movement.Take the MLFMM and RL-GO results as the reference values, the µ and σ for our method are presented in Table 1.It is obvious that the µ and σ of our method are about 0.36~0.88 and 1.67~2.81.It indicates that these results using different methods match very well.

EM Scattering of the Moving Target
In Figure 10, the backscattering of the dihedral angle and combined model are compared with our method and other solutions (MLFMM, RL-GO, and the traditional SBR). Figure 10a is the scattering of the dihedral angle composed of two 0.5 m × 0.5 m plates with HH polarization.The dihedral angle rotates around the  -axis with an angular velocity of  = /18 rad/s.Where the incident frequency is  = 10 GHz, the incident angle is   = 0°, and the azimuth angle is   = 0°.Figure 10b is the scattering of the combined model composed of two cubes with side lengths of 0.4 m and 0.2 m, respectively.The combined model rotating around the y -axis with an angular velocity is  = /2 rad/s.
Where the incident frequency is  = 8 GHz, the incident angle is   = 90°, the azimuth angle is   = 0°, and the polarization is VV.In this section, since there is need to evaluate the effectiveness of our method, some statistical concepts are utilized to measure the deviation degree of our method from other solutions.The mean absolute error  and its standard deviation  [7] are defined as where   are the differences between the predicted RCS values of our method and other solutions; and N is the total number of moments for target movement.Take the MLFMM and RL-GO results as the reference values, the  and  for our method are presented in Table 1.It is obvious that the  and  of our method are about 0.36~0.88 and 1.67~2.81.It indicates that these results using different methods match very well.SBR. Figure 11a  Where the incident frequency is  = 10 GHz, the incident angle is   = 90°, and the azimuth angle is   = 90°.Figure 11c is the scattering of the aircraft with HH polarization.The aircraft is performing a climbing motion at a speed of =120 m/s.Where the incident frequency is  = 10 GHz, the incident angle is   = 120°, and the azimuth angle is   = 90°.In Table 1, it is evident that our method shows good consistency with the MLFMM and RL-GO.In Table 2, the computation cost of Figures 10 and 11 is compared between our method and other solutions.It can be seen from Table 2 that the method adopted in this paper has significant improvement in computational speed compared with the traditional SBR and MLFMM or RL-GO.And with the number of facets increasing, the acceleration effect is more obvious.The speedup of our method relative to traditional SBR is also given in Table 2.In Figure 11a,b, the backscattering of the trihedral corner reflector and missile is compared using our method and other solutions (MLFMM, RL-GO, and the traditional SBR). Figure 11c is the backscattering of aircraft compared using our method and the traditional SBR. Figure 11a is the scattering of the trihedral corner reflector composed of three right-angled triangles (right-angle side 0.3 m) with HH polarization.The trihedral corner reflector rotates around the rotation axis n = [−1, 1, 0] with an angular velocity ω = π/18 rad/s.Where the incident frequency is f = 16 GHz, the incident angle is θ i = 0 • , and the azimuth angle is φ i = 45 • .Figure 11b is the scattering of the missile with VV polarization.The missile rotates around the z-axis with an angular velocity ω = π/180 rad/s.Where the incident frequency is f = 10 GHz, the incident angle is θ i = 90 • , and the azimuth angle is φ i = 90 • .Figure 11c is the scattering of the aircraft with HH polarization.The aircraft is performing a climbing motion at a speed of v= 120 m/s.Where the incident frequency is f = 10 GHz, the incident angle is θ i = 120 • , and the azimuth angle is φ i = 90 • .In Table 1, it is evident that our method shows good consistency with the MLFMM and RL-GO.In Table 2, the computation cost of Figures 10 and 11 is compared between our method and other solutions.It can be seen from Table 2 that the method adopted in this paper has significant improvement in computational speed compared with the traditional SBR and MLFMM or RL-GO.And with the number of facets increasing, the acceleration effect is more obvious.The speedup of our method relative to traditional SBR is also given in Table 2.

Micro-Doppler of the Wingless Missile
In Figure 12, the instantaneous HRRP and micro-Doppler of the wingless missile with precession are simulated using our method.The missile rotates counterclockwise with the angular velocity ω = 2π rad/s around the rotation axis n = [0.1961, 0.9806, 0].Where the incident frequency is f = 6 GHz, the incident angle is θ i = 45 • , the azimuth angle is φ i = 0 • , and the polarization is HH.In order to verify the micro-Doppler spectra in Figure 12d, the equivalent strong scattering point of the target can be extracted by calculating the instantaneous HRRP.It is the main component that causes the micro-Doppler at the time.The theoretical Doppler formula is given by Equation (19).

Micro-Doppler of the Wingless Missile
In Figure 12, the instantaneous HRRP and micro-Doppler of the wingless missile with precession are simulated using our method.The missile rotates counterclockwise with the angular velocity  = 2 rad/s around the rotation axis  ̂= [0.1961,0.9806,0].Where the incident frequency is  = 6 GHz, the incident angle is   = 45°, the azimuth angle is   = 0°, and the polarization is HH.In order to verify the micro-Doppler spectra in Figure 12d, the equivalent strong scattering point of the target can be extracted by calculating the instantaneous HRRP.It is the main component that causes the micro-Doppler at the time.The theoretical Doppler formula is given by Equation (19). Figure 12a-c is the instantaneous HRRPs of the missile with precession at  = 0.4 s, 0.6s, and 0.8s.In Figure 12a-c, the equivalent strong scattering points are concentrated Figure 12a-c is the instantaneous HRRPs of the missile with precession at t = 0.4 s, 0.6s, and 0.8s.In Figure 12a-c, the equivalent strong scattering points are concentrated with the head (Point 2, 4, 6) and tail (Point 1, 3, 5) of the missile, and the intensity of the Point 2 (Point 4, 6) is weaker than Point 1 (Point 3, 5). Figure 12d is the verification of the theoretical Doppler and micro-Doppler spectra.In Figure 12d, the red scattering points are the theoretical Doppler calculated by points at the head, and the black scattering points are the theoretical Doppler calculated by points at the tail.It can be seen from Figure 12d that the micro-Doppler spectra is consistent with the theoretical Doppler of the equivalent scattering point for most moments.

Micro-Doppler of the Winged Missile
In this example, the instantaneous HRRP and micro-Doppler of the wingless missile with spin motion are presented in Figure 11.The missile rotates counterclockwise with the angular velocity ω = 2π rad/s around the y-axis to simulate spin motion.Figure 11a is the instantaneous HRRP of the missile with spin motion at t = 0.25 s Figure 11b is the micro-Doppler spectra of the missile with spin motion.Where the incident frequency is f = 10 GHz, the incident angle is θ i = 90 • , the azimuth angle is φ i = 0 • , and the polarization is VV.
As shown in Figure 13a,b, the micro-Doppler spectra of the spin missile has two identical cycles and each cycle lasts 0.25 s.This is because the missile has four symmetrical tails.Thus, one cycle is generated every 0.25 s when the spin angular velocity is 2π rad/s.It can be also seen that the micro-Doppler frequency shift and the overall scattering intensity reaches the highest at t = 0.25 s.This is because at this moment the surface of the tail is vertical to the direction of the incident wave, as shown in Figure 13a.The theoretical Doppler of equivalent strong scattering points 1, 2, and 3 in Figure 13a are calculated, and compared with the micro-Doppler spectra shown in Figure 13b.12d is the verification of the theoretical Doppler and micro-Doppler spectra.In Figure 12d, the red scattering points are the theoretical Doppler calculated by points at the head, and the black scattering points are the theoretical Doppler calculated by points at the tail.It can be seen from Figure 12d that the micro-Doppler spectra is consistent with the theoretical Doppler of the equivalent scattering point for most moments.

Micro-Doppler of the Winged Missile
In this example, the instantaneous HRRP and micro-Doppler of the wingless missile with spin motion are presented in Figure 11.The missile rotates counterclockwise with the angular velocity  = 2 rad/s around the  -axis to simulate spin motion.Figure 11a is the instantaneous HRRP of the missile with spin motion at  = 0.25 s Figure 11b is the micro-Doppler spectra of the missile with spin motion.Where the incident frequency is  = 10 GHz, the incident angle is   = 90°, the azimuth angle is   = 0°, and the polarization is VV.
As shown in Figure 13a,b, the micro-Doppler spectra of the spin missile has two identical cycles and each cycle lasts 0.25.This is because the missile has four symmetrical tails.Thus, one cycle is generated every 0.25 s when the spin angular velocity is 2 rad/s.It can be also seen that the micro-Doppler frequency shift and the overall scattering intensity reaches the highest at  = 0.25 s.This is because at this moment the surface of the tail is vertical to the direction of the incident wave, as shown in Figure 13a.The theoretical Doppler of equivalent strong scattering points 1, 2, and 3 in Figure 13a are calculated, and compared with the micro-Doppler spectra shown in Figure 13b.In Figure 14a, the missile rotates around  ̂= [0.1961,0.9806,0] to simulate cone motion and the rotational angular velocity of cone motion is  = 2 rad/s.The precession in Figure 14b is a superposition of the spin motion in Figure 13 and the cone motion in Figure 14a.In Figure 14a, the missile rotates around n = [0.1961, 0.9806, 0] to simulate cone motion and the rotational angular velocity of cone motion is ω = 2π rad/s.The precession in Figure 14b is a superposition of the spin motion in Figure 13 and the cone motion in Figure 14a.
In Figure 14b, the micro-Doppler of the precession motion can be clearly distinguished into two parts: the missile body (feature 1) and the tail (feature 2).The time interval between every feature 2 is 0.125 s.This is because when the missile spins, the micro-motion of the tail will produce four features 2. At the same time, the cone motion of the missile also causes four features 2 produced by the tail.Thus, there are eight features 2 in one progression cycle.Since the main part of the missile is also involved in a coning motion, the scattering intensity of eight features 2 is weaker than feature 1.In addition, due to the influence of the motion of the missile body, the eight features 2 in the micro-Doppler spectra are not symmetrical.In Figure 14b, the micro-Doppler of the precession motion can be clearly distinguished into two parts: the missile body (feature 1) and the tail (feature 2).The time interval between every feature 2 is 0.125s.This is because when the missile spins, the micromotion of the tail will produce four features 2. At the same time, the cone motion of the missile also causes four features 2 produced by the tail.Thus, there are eight features 2 in one progression cycle.Since the main part of the missile is also involved in a coning motion, the scattering intensity of eight features 2 is weaker than feature 1.In addition, due to the influence of the motion of the missile body, the eight features 2 in the micro-Doppler spectra are not symmetrical.

Micro-Doppler of the Aircraft
Figure 15 is the micro-Doppler spectra of the missile with a flipping motion in different incident angles.The aircraft rotates with the angular velocity  = 0.2 rad/s around the y -axis to simulate a flipping motion.Where the incident frequency is  = 6 GHz, the incident angle is   = 90°, 60°, 45°, and 30°; the azimuth angle is   = 0°; and the polarization is HH.With the decrease in the incident angle, the variation in time-frequency ridges is reflected that the maximum of the Doppler shift is shifted to the left.The maximum of the Doppler shift is shifted to the left by 0.77~0.84s when the incident angle   decreases by 30°.In addition, as shown in Figure 15a, time-frequency ridges in the first half of the cycle contain more micro-Doppler components than the last half.This is because the Doppler shift is contributed to by the micro-motion of the two vertical tails and other physical structures on the front of the fuselage; in addition, the front of the fuselage is more complex than the back, so the amount of Doppler shift information embodied in time-frequency ridges is more.

Micro-Doppler of the Aircraft
Figure 15 is the micro-Doppler spectra of the missile with a flipping motion in different incident angles.The aircraft rotates with the angular velocity ω = 0.2π rad/s around the yaxis to simulate a flipping motion.Where the incident frequency is f = 6 GHz, the incident angle is θ i = 90 • , 60 • , 45 • , and 30 • ; the azimuth angle is φ i = 0 • ; and the polarization is HH.With the decrease in the incident angle, the variation in time-frequency ridges is reflected that the maximum of the Doppler shift is shifted to the left.The maximum of the Doppler shift is shifted to the left by 0.77 ∼ 0.84 s when the incident angle θ i decreases by 30 • .In addition, as shown in Figure 15a, time-frequency ridges in the first half of the cycle contain more micro-Doppler components than the last half.This is because the Doppler shift is contributed to by the micro-motion of the two vertical tails and other physical structures on the front of the fuselage; in addition, the front of the fuselage is more complex than the back, so the amount of Doppler shift information embodied in time-frequency ridges is more.In Figure 14b, the micro-Doppler of the precession motion can be clearly distinguished into two parts: the missile body (feature 1) and the tail (feature 2).The time interval between every feature 2 is 0.125s.This is because when the missile spins, the micromotion of the tail will produce four features 2. At the same time, the cone motion of the missile also causes four features 2 produced by the tail.Thus, there are eight features 2 in one progression cycle.Since the main part of the missile is also involved in a coning motion, the scattering intensity of eight features 2 is weaker than feature 1.In addition, due to the influence of the motion of the missile body, the eight features 2 in the micro-Doppler spectra are not symmetrical.

Micro-Doppler of the Aircraft
Figure 15 is the micro-Doppler spectra of the missile with a flipping motion in different incident angles.The aircraft rotates with the angular velocity  = 0.2 rad/s around the y -axis to simulate a flipping motion.Where the incident frequency is  = 6 GHz, the incident angle is   = 90°, 60°, 45°, and 30°; the azimuth angle is   = 0°; and the polarization is HH.With the decrease in the incident angle, the variation in time-frequency ridges is reflected that the maximum of the Doppler shift is shifted to the left.The maximum of the Doppler shift is shifted to the left by 0.77~0.84s when the incident angle   decreases by 30°.In addition, as shown in Figure 15a, time-frequency ridges in the first half of the cycle contain more micro-Doppler components than the last half.This is because the Doppler shift is contributed to by the micro-motion of the two vertical tails and other physical structures on the front of the fuselage; in addition, the front of the fuselage is more complex than the back, so the amount of Doppler shift information embodied in time-frequency ridges is more.

Conclusions
In this paper, the EM scattering of a moving target is solved using a tailored SBR method based on the dynamic spatial division technique.In this technique, the dynamic spatial division is established in the LCS to store the facet information and remain relatively stationary with the target.In comparison with the traditional SBR method, this technique avoids repetitive spatial division at each moment in the GCS, which greatly decreases the consumption of computing resources.In our numerical simulations, the verification and computation efficiency comparison are given by our method and other solutions (MLFMM, RL-GO, and traditional SBR method).Moreover, the micro-Doppler features are extracted and analyzed using the time-frequency analysis technique, which includes the precession and spin of the missile, and the rotation of the aircraft.Meanwhile, the micro-Doppler spectra of the target is also compared with the theoretical Doppler of the equivalent strong scattering points, which are obtained using the instantaneous HRRP.In our future work, the tailored SBR method in this paper will be improved to solve and analyze the scattering and micro-Doppler of multiple moving targets.

Conclusions
In this paper, the EM scattering of a moving target is solved using a tailored SBR method based on the dynamic spatial division technique.In this technique, the dynamic spatial division is established in the LCS to store the facet information and remain relatively stationary with the target.In comparison with the traditional SBR method, this technique avoids repetitive spatial division at each moment in the GCS, which greatly decreases the consumption of computing resources.In our numerical simulations, the verification and computation efficiency comparison are given by our method and other solutions (MLFMM, RL-GO, and traditional SBR method).Moreover, the micro-Doppler features are extracted and analyzed using the time-frequency analysis technique, which includes the precession and spin of the missile, and the rotation of the aircraft.Meanwhile, the micro-Doppler spectra of the target is also compared with the theoretical Doppler of the equivalent strong scattering points, which are obtained using the instantaneous HRRP.In our future work, the tailored SBR method in this paper will be improved to solve and analyze the scattering and micro-Doppler of multiple moving targets.

Figure 1 .
Figure 1.Geometry of moving target with translation and rotation.

Figure 1 .
Figure 1.Geometry of moving target with translation and rotation.

Figure 3 .
Figure 3. Flowchart of the dynamic spatial division.

Figure 4 Figure 3 .
Figure4illustrates the geometry of the illuminated facet with rotation and translation.The facet _  is illuminated by the incident electric field    in   .The TE and TM waves can be expressed as

Figure 4
Figure 4 illustrates the geometry of the illuminated facet with rotation and translation.The facet Tri_t n is illuminated by the incident electric field E G inc in xyz GCS .The TE and TM waves can be expressed as

Figure 4 .
Figure 4. Geometry of illuminated facet with rotation and translation.

Figure 4 .
Figure 4. Geometry of illuminated facet with rotation and translation.

Figure 5 .
Figure 5. Flowchart for solving the scattering field of Tri_t n .

Figure 5 .
Figure 5. Flowchart for solving the scattering field of .

Figure 7 .
Figure 7. Schematic of the sliding window function.

Figure 7 .
Figure 7. Schematic of the sliding window function.

Figure 8 .
Figure 8. Flowchart of the proposed method.The red dashed box indicates the calculation of the scattering field by the tailored SBR method.The green dashed box indicates the micro-Doppler obtained by the STFT.

Figure 8 .
Figure 8. Flowchart of the proposed method.The red dashed box indicates the calculation of the scattering field by the tailored SBR method.The green dashed box indicates the micro-Doppler obtained by the STFT.
is the scattering of the trihedral corner reflector composed of three rightangled triangles (right-angle side 0.3 m) with HH polarization.The trihedral corner reflector rotates around the rotation axis  = [−1,1,0] with an angular velocity  = /18 rad/s.Where the incident frequency is  = 16 GHz, the incident angle is   = 0°, and the azimuth angle is   = 45°.Figure 11b is the scattering of the missile with VV polarization.The missile rotates around the  -axis with an angular velocity  = /180 rad/s.

Table 1 .
Comparison of effectiveness between our method and other solutions.

Table 2 .
Computation cost comparison between our method and other solutions in Figures10 and 11.

Table 2 .
Computation cost comparison between our method and other solutions in Figures10 and 11.