### I. Introduction

_{1}) [15] and a neural network classifier (NNC

_{2}) [16] were compared. The classification results showed that the MD bandwidth was the most adequate for short observation time

*T*

_{ob}but that to exploit the time-periodic nature of MD,

*T*

_{ob}must be extended. The correct classification ratio Pc ≈ 100% was achieved by using only three features. NNC

_{1}was found to be more appropriate than NNC

_{2}for drone–bird classification.

### II. Mathematical Modeling and Proposed Method

### 1. Mathematical Modeling of the Rotating Drone Blade and Flapping Bird Wing

*a*

_{b}and length 2

*b*

_{b}, which are rotated around the x-axis, one by ±45° and another by −45° (Fig. 1). The coordinates of a blade tip rotating clockwise at frequency

*f*

_{b}around the z-axis are obtained as

##### (1)

*φ*and an elevation

*θ*(Fig. 2), the RCS of the flat plate rotated by 45° is analytically given as [13]:

*f*

_{w}, the analysis assumed that the lower muscle tip at (

*x*

_{1},

*y*

_{1},

*z*

_{1}) rotated by

*θ*

_{1}= 40°cos(2

*πf*

_{w}

*t*)

*+*15° and the upper muscle tip at (

*x*

_{2},

*y*

_{2},

*z*

_{2}) rotated by

*θ*

_{2}= 30° cos(2

*πf*

_{w}

*t*)

*+*40°. The forward–backward movement of the upper wing was modeled by rotating (

*x*

_{2},

*y*

_{2},

*z*

_{2}) by

*φ*

_{w}= 20°sin(2

*πf*

_{w}

*t*)

*+*40° in the

*x*-

*y*plane. Thus, the two-muscle position can be expressed as

##### (9)

*r*

_{1}and

*r*

_{2}are

*θ*

_{w}with respect to the major axis is analytically given as

*a*

_{w}and

*b*

_{w}are the lengths of the minor and major axes, respectively [13]. Assuming that the vector of the major axis of the lower and upper muscles is

*p̄*

_{i}(

*i*= 1 for the lower and 2 for the upper muscle),

*θ*

_{w}is calculated as

### 2. Radar Signal Modeling and Principle of MD

*T*

_{chirp}[13] (Fig. 5).

*T*

_{chirp}, is expressed as

*f*

_{0}is the minimum frequency and

*K*

_{r}is the chirp rate (Fig. 5). The bandwidth

*B*,

*K*

_{r}, and

*T*

_{chirp}are related as

*τ*= 2

*r*/

*c*and

*c*is the speed of the light. For the bistatic observation scenario, 2

*r*is replaced by

*r*

_{t}

*+ r*

_{r}, where

*r*

_{t}is the distance from the transmitting radar to the target and

*r*

_{r}is the distance to the target receiving radar.

*K*

_{r}

*τ*is the beat frequency:

*K*

_{r}

*τ*

^{2}/2 has been removed.

*A*

_{0}is observed for a CPI

*T*

_{CPI}, (11)–(15) are used for each slow time

*t*

_{s}=

*kT*

_{chirp}during

*T*

_{CPI}. Thus, the received signal can be expressed as

*R*(

*t*

_{s}) is the distance to the target at

*t*

_{s}. Dechirping changes

*s*

_{0}(

*t*,

*t*

_{s}) to

*s*

_{0d}(

*t*,

*t*

_{s}) as

*t*-domain to the

*f*-domain yields a dechirped signal in HRRP at

*t*

_{s}(i.e.,

*s*

_{r}(

*f*,

*t*

_{s}):

*p*

_{r}(·) is a sinc function formed by a rectangular window of length

*T*

_{chirp}. The down-range location of the target can be calculated by determining the peak of

*p*

_{r}(·).

*s*

_{MD}(

*t*

_{s}) =

*s*

_{r}(

*K*

_{r}

*τ*

_{p},

*t*

_{s}), where

*τ*

_{p}is the time delay of the signal reflected from the target.

*s*

_{MD}(

*t*

_{s}) is the time-varying signal of the target caused by the micro-motion and sampled at

*t*

_{s}. However,

*s*

_{MD}(

*t*

_{s}) is time-varying and thus should be analyzed in the TF domain rather than the frequency domain. In this study, we used the short-time Fourier transform (STFT), which is simple to implement and free from cross-term interference. STFT transforms

*s*

_{MD}(

*t*

_{s}) to MD [18, 19]:

*f*

_{s}is the frequency and

*W*is the window function. To reduce the sidelobe caused by the rectangular window, we used the Hamming window. In discrete form, (22) is expressed as

*k*is the time index,

*L*is the window length,

*p*is the sampled time index, and

*q*is the sampled frequency index.

### 3. Feature Vectors for Classification

*f*

_{rcs}of the MD signal variation in the time domain represents the rotation frequency of the drone blades; thus,

*f*

_{rcs}is used as a feature and is defined as

*rcs*(

*τ*) = |

*s*

_{MD}(

*τ*)|.

*f*

_{h}of the MD change in the time domain differs considerably between the drone and the bird; thus, this change is used as a feature. To obtain

*f*

_{h}, the frequency-domain image

*I*

_{ff}is obtained by the FT of the MD image in (23) in the time domain as

_{q}is the FT in

*q*for each

*p*. Then,

*f*

_{h}is obtained as

*f*

_{v}is used as a feature and can be defined as

*f*

_{max}(

*s*(

*k*),

*γ*) is the maximum

*k*corresponding to a value larger than

*γ*% of the maximum

*s*(

*k*).

*cor*± of ±MD frequencies is used as a feature. The correlation

*cor*

_{±}is defined as

*X*and

*Y*are the MD images clipped from

*I*

_{ff}for positive and negative MD frequencies, respectively;

*X̄*and

*Ȳ*are the averages of

*X*and

*Y*; and

*V*(·) is the variance.

*x-y*plane, so the monostatic/bistatic MD images of MD drones are similar, whereas bird wings flap in the

*z-y*plane. The correlation

*cor*

_{mb}obtained using the monostatic and bistatic images is used as a feature:

*I*

_{ff,m}and

*I*

_{ff,b}are monostatic and bistatic

*I*

_{ff}, respectively, and

*Ī*

_{ff,m}and

*Ī*

_{ff,b}are the averages of

*I*

_{ff,m}and

*I*

_{ff,b}.

*f*

_{pk,v}is used as a feature:

*f*

_{pk,h}is used as a feature and is defined as

*f*

_{pk,2D}in the normalized

*I*

_{ff}is

*ent*

_{v}of

*Z*

_{v}is used as a feature and can be expressed as

*Z*

_{h}is also used as a feature:

*ent*

_{2D}, as the 2D entropy of

*I*

_{N}, is used as a feature:

### 4. Construction of the Training Database, Classifier, and Overall Procedure

*θ*

_{t}in increments of Δ

*θ*

_{t}. In addition, for each flight direction,

*f*

_{b}of the drone blade was sampled between

*f*

_{b,1}and

*f*

_{b,2}in increments of Δ

*f*

_{b}, and

*f*

_{w}of the bird wing was sampled between

*f*

_{w,1}and

*f*

_{w,2}in increments of Δ

*f*

_{w}. In the bistatic observation scenario, a receiving radar was located at a distance

*r*

_{b}from the transmitting radar, and the training database was constructed using the same grid points and motion parameters.

_{1}and NNC

_{2}, and their classification accuracies were compared. NNC

_{1}uses a simple Euclidean norm as [15]

*χ̄*

_{u}is a test vector and

*f̄*

_{i}is a training vector that belongs to the

*i*th class. The class

*i*that yields the minimum Euclidean norm is the class that includes

*χ̄*

_{u}.

_{2}, we constructed a neural network composed of two hidden layers: the first layer had five nodes and the second had three [16] (Fig. 8). Training was conducted using a backpropagation algorithm; the output was set to “1” for drones and “0” for birds. Because the decrease in the training error after 5,000 iterations was slight, 10,000 iterations were conducted to prevent time consumption. The classification was performed by following the general target recognition procedure (Fig. 9). When the test data were obtained, the target location and motion parameter were not known; thus, the target was randomly positioned in the training space, and the micro-motion was modeled using a parameter that was randomly set within the range of that parameter in the training data. Then, an MD image was obtained following Fig. 6. Classification was conducted using the defined features (Table 1), and the importance of each feature and the classification accuracy of the two classifiers were evaluated.

### III. Simulation Result

### 1. Simulation Condition

*B*= 150 MHz; however, we reduced it to 0.150 MHz for fast computation, because the MD image is not relevant to range resolution but sampling the target using PRF, which exists in a range bin (Fig. 10).

*f*

_{b,1}= 40 Hz and

*f*

_{b,2}= 50 Hz. The blade dimensions were selected as

*a*

_{b}= 8 cm and

*b*

_{b}= 1.5 cm (Fig. 2). For the bird, we assumed a wing span = 60 cm (i.e.,

*a*

_{w}=

*b*

_{w}= 15 cm) at a flapping frequency

*f*

_{w,1}= 4 Hz or

*f*

_{w,2}= 7 Hz. As mentioned above, the MD of a bird is very different from that of a drone (Fig. 11). Due to the asymmetric motion of a bird, its MD is not symmetric and blade flashes are not observed. In addition, the period of MD is much smaller than that of the blade.

*x-*axis between 300 and 1,500 m was divided into ten equal intervals; the range on the

*z*-axis between 10 and 1,000 m was divided into five equal intervals; and the azimuth angle between −45° and

*+*45° was divided into five equal intervals (Fig. 12(a)). At each grid point, the flight direction was set between −

*θ*

_{t}= −30° and

*+θ*

_{t}= 30° in increments of Δ

*θ*

_{t}= 12°. For each flight direction on each grid point, the MD image was obtained by varying

*f*

_{b}and

*f*

_{w}in increments of 1 Hz. The test data were randomly obtained. Identical simulations were conducted for bistatic observations with a receiving radar located at 1,000 m on the

*y-*axis from the transmitting radar at the origin.

*x*≤ 1,500 m, 10 m ≤

*z*≤ 100 m, 0° ≤ azimuth angle ≤ 45°), the target was flown in a random direction at −30 ≤

*θ*

_{t}≤ 30° with a random micro-motion frequency (Fig. 12(b)). For

*T*

_{CPI}= 0.0385 seconds, 2,500 test images were used per target; for

*T*

_{CPI}= 1.001 seconds, this number was decreased to 300 due to the increased computation time. As in the training phase, the bistatic test data were obtained using the same parameters. The classification accuracy was expressed as the correct classification percentage:

*M*

_{c}is the number of correct classifications and

*M*

_{te}is the number of test samples.

### 2. Classification Results

#### 1) Classification results for the stationary target

_{1}was used as the classifier. The

*P*

_{c}values for both

*T*

_{CPI}were proportional to the SNR, except for

*f*

_{21}for

*T*

_{CPI}= 0.0385 seconds.

*P*

_{c}was approximately 100% at SNR = 20 dB (Fig. 13). Comparing the features for each

*T*

_{CPI}, the highest

*P*

_{c}s at low SNRs were obtained using

*f*

_{2}and

*f*

_{5}(i.e., MD frequency

*f*

_{h}for monostatic and bistatic scenarios), which is a rather unexpected result. These results are attributed to the fact that for the short CPI, the frequency of MD change of the bird is approximately 0 Hz due to the slow variation in bird flapping, whereas the blade rotation is fast enough to yield a certain value of

*f*

_{h}. As a result,

*P*

_{c}was high.

*T*

_{CPI}values, the features that exploited the MD periodicity and the projected sum onto the horizontal axis (i.e.,

*f*

_{1},

*f*

_{7},

*f*

_{10},

*f*

_{13},

*f*

_{16}) were significantly improved at

*T*

_{CPI}= 1.0 second because the increased observation time allowed the periodicity to be represented by FT along the time domain (= horizontal axis).

*f*

_{2}and

*f*

_{5}are also the features representing MD periodicity; however, the

*P*

_{c}s of these vectors for

*T*

_{CPI}= 1.0 second were lower than those for

*T*

_{CPI}= 0.0385 seconds because of the large zero–nonzero difference.

*f*

_{21}increased for

*T*

_{CPI}because the increased observation time represented the MD image in detail. The entropy-related features (

*f*

_{7}–

*f*

_{12}) yielded a high

*P*

_{c}value but were sensitive to noise.

#### 2) Classification results for the moving target

*v*and acceleration

*a*of a moving target by randomly selecting 0 ≤

*v*≤ 10 m/s and 0 ≤

*a*≤ 10 m/s

^{2}.

*v*of the rigid body shifted the MD and

*a*tilted the MD upward or downward, the test MD image was significantly different from the training image (Fig. 14). Other simulation conditions were the same as those mentioned in Section III-2-1.

*v*and

*a*should be accurately estimated and compensated for (Fig. 15). Compared to the result of the stationary target (Fig. 13),

*P*

_{c}of the features decreased significantly. Due to the change in MD frequency caused by

*v*and

*a*, the RCS frequency varied; thus,

*f*

_{1}and

*f*

_{4}yielded poor

*P*

_{c}. In addition,

*P*

_{c}of the features (

*f*

_{7},

*f*

_{10},

*f*

_{13},

*f*

_{16}) using the horizontally projected sum decreased considerably due to the tilted MD image in the TF domain; the amount of decrease was larger at

*T*

_{CPI}= 1.001 seconds than that at

*T*

_{CPI}= 0.0385 seconds because the tilt of the MD frequency in the TF domain was larger at

*T*

_{CPI}= 1.001 seconds than that at

*T*

_{CPI}= 0.0385 seconds. Thus, the test

*I*

_{ff}became totally different from the training

*I*

_{ff}s after FT on the horizontal axis.

*T*

_{CPI}affected several

*P*

_{c}values. At a short

*T*

_{CPI}, the

*P*

_{c}values were high because these features are nonzero for a bird and close to zero for a drone. However, at a long

*T*

_{CPI}, the frequency of the MD change affected by

*v*and

*a*was represented by FT, so

*P*

_{c}decreased to approximately 50%. The MD bandwidth (

*f*

_{3}and

*f*

_{6}) was less affected because the MD bandwidth of the drone was much larger than that of the drone. The entropy obtained using the projected sum onto the vertical axis (

*f*

_{7}and

*f*

_{10}, “fast”) was less affected by

*v*and

*a*than the MD bandwidth because the entropy is not determined by the MD bandwidth but by the overall distribution. The

*P*

_{c}values of

*f*

_{19}and

*f*

_{20}were also affected by

*v*and

*a*, and those of

*f*

_{21}increased for

*T*

_{CPI}due to the similarity between the monostatic and bistatic MD images for the increased observation time.

*v*and

*a*were assumed to have been perfectly compensated for, the number of targets was randomly selected between 1 and 5 inclusive, and the classifications were conducted using the same parameters as those used in the first classification.

*v*and

*a*were completely removed,

*P*

_{c}(Fig. 16) was slightly lower than the values shown in Fig. 13. The

*P*

_{c}values for both

*T*

_{CPI}values were proportional to SNR except for

*f*

_{21}for

*T*

_{CPI}= 0.0385 seconds.

*f*

_{2}and

*f*

_{5}, which represent MD frequency

*f*

_{h}, yielded the highest

*P*

_{c}at low SNRs due to zero and nonzero MD frequencies.

*f*

_{1},

*f*

_{7},

*f*

_{10},

*f*

_{13}, and

*f*

_{16}, which represent the projected sum onto the horizontal axis, improved significantly for

*T*

_{CPI}= 1.0 second due to the increased observation time.

*f*

_{2}and

*f*

_{5}for

*T*

_{CPI}= 1.0 second yielded lower

*P*

_{c}than those for

*T*

_{CPI}= 0.0385 seconds because of the large zero–nonzero difference. The entropy-related features (

*f*

_{7}–

*f*

_{12}) yielded high

*P*

_{c}but were sensitive to noise, and

*f*

_{21}increased for

*T*

_{CPI}due to the increased observation time.

#### 3) Classification results for the fusion of features

*P*

_{c}s of NNC

_{1}and NNC

_{2}were compared for various SNRs using the training and test data same as those presented in Section III-2-1. The number of combinations was large, many of which had

*P*

_{c}

*≈*100%, so combinations of two and three among the 21 features were analyzed, and the optimal combinations that yielded

*P*

_{c}

*≈*100% were found.

*P*

_{c}at SNR = 0 dB and 20 dB, the best five combinations obtained

*P*

_{c}≥ 95% for SNR = 0 dB and

*P*

_{c}≥ 98% for SNR = 20 dB (Fig. 17). As in the classifications conducted above,

*P*

_{c}was proportional to SNR. For

*T*

_{CPI}= 0.0385 seconds, NNC

_{1}using

*f*

_{2}

*+ f*

_{5},

*f*

_{2}

*+ f*

_{15},

*f*

_{2}

*+ f*

_{18},

*f*

_{5}

*+ f*

_{15}, and

*f*

_{5}

*+ f*

_{18}yielded high

*P*

_{c}

*≈*100% for SNR = 20 dB, and for NNC

_{2},

*f*

_{2}

*+ f*

_{5},

*f*

_{2}

*+ f*

_{7},

*f*

_{5}

*+ f*

_{10}, and

*f*

_{5}

*+ f*

_{13}were more effective than other features. Comparing NNC

_{1}and NNC

_{2}, the

*P*

_{c}of NNC

_{1}was higher than that of NNC

_{2}. This is because the designed NNC

_{2}(Fig. 8) was rather simple and the training data were over-fitted to NNC

_{2}; thus,

*P*

_{c}was lower because of the lack of the generalization capability.

*P*

_{c}s for

*T*

_{CPI}= 1.001 seconds were much higher for both classifiers than for

*T*

_{CPI}= 0.0385 seconds. For NNC

_{1},

*P*

_{c}s of

*f*

_{2}

*+ f*

_{3},

*f*

_{2}

*+ f*

_{6}, and

*f*

_{3}

*+ f*

_{7}were approximately 100% for all SNRs, and

*f*

_{3}

*+ f*

_{13}was sensitive to noise. For NNC

_{2},

*P*

_{c}s of

*f*

_{1}

*+ f*

_{4},

*f*

_{1}

*+ f*

_{20},

*f*

_{4}

*+ f*

_{15},

*f*

_{4}

*+ f*

_{19}, and

*f*

_{4}

*+ f*

_{20}were slightly lower than those of NNC

_{2}at SNR ≥ −10 dB but close to 100%. The reason for the improvement was that for

*T*

_{CPI}= 1.001 seconds, the features

*f*

_{2}and

*f*

_{7}for NNC

_{1}and the features

*f*

_{1},

*f*

_{4},

*f*

_{19}, and

*f*

_{20}better represented the periodicity of MD than that at

*T*

_{CPI}= 0.0385 seconds.

*P*

_{c}value compared to those using combinations of two features, so the requirement to sift the classification results was changed. For

*T*

_{CPI}= 0.0385 seconds, the best five combinations that satisfied

*P*

_{c}≥ 90% for SNR = −10 dB and

*P*

_{c}≥98% for SNR = 20 dB were displayed (Fig. 18). For

*T*

_{CPI}= 1.001 seconds, the combinations for NNC

_{1}that yielded

*P*

_{c}= 100% for all SNRs were used, and for NNC

_{2}, the combinations that yielded

*P*

_{c}≥ 95% for all SNRs were used.

*T*

_{CPI}= 0.0385 seconds

*P*

_{c}improved only slightly at low SNRs because of the poor representation of the target by the features that exploit periodicity. Using NNC

_{1},

*f*

_{2}

*+ f*

_{5}

*+ f*

_{15}and

*f*

_{2}

*+ f*

_{5}

*+ f*

_{15}satisfied the requirement; using NNC

_{2},

*f*

_{2}

*+ f*

_{5}

*+ f*

_{10},

*f*

_{2}

*+ f*

_{5}

*+ f*

_{13},

*f*

_{2}

*+ f*

_{5}

*+ f*

_{16},

*f*

_{2}

*+ f*

_{5}

*+ f*

_{20}, and

*f*

_{2}

*+ f*

_{7}

*+ f*

_{13}satisfied the requirement.

*P*

_{c}s of NNC

_{1}were slightly higher than those of NNC

_{2}, although the number of feature combinations for NNC

_{1}was smaller than that for NNC

_{2}. For both classifiers, features

*f*

_{2}and

*f*

_{5}worked well because of the zero–nonzero MD relationship (Section III-2-1).

*T*

_{CPI}= 1.001 seconds, 23 feature combinations for NNC

_{1}yielded

*P*

_{c}s = 100% for all SNRs:

*f*

_{2}

*+ f*

_{3}

*+ f*

_{7},

*f*

_{2}

*+ f*

_{3}

*+ f*

_{9},

*f*

_{2}

*+ f*

_{3}

*+ f*

_{10},

*f*

_{2}

*+ f*

_{3}

*+ f*

_{12},

*f*

_{2}

*+ f*

_{3}

*+ f*

_{13},

*f*

_{2}

*+ f*

_{3}

*+ f*

_{14},

*f*

_{2}

*+ f*

_{3}

*+ f*

_{15},

*f*

_{2}

*+ f*

_{3}

*+ f*

_{16},

*f*

_{2}

*+ f*

_{3}

*+ f*

_{18},

*f*

_{2}

*+ f*

_{3}

*+ f*

_{19},

*f*

_{2}

*+ f*

_{3}

*+ f*

_{20},

*f*

_{2}

*+ f*

_{3}

*+ f*

_{21},

*f*

_{2}

*+ f*

_{7}

*+ f*

_{9},

*f*

_{2}

*+ f*

_{7}

*+ f*

_{10},

*f*

_{2}

*+ f*

_{7}

*+ f*

_{12},

*f*

_{2}

*+ f*

_{7}

*+ f*

_{13},

*f*

_{2}

*+ f*

_{7}

*+ f*

_{16},

*f*

_{2}

*+ f*

_{7}

*+ f*

_{21},

*f*

_{2}

*+ f*

_{9}

*+ f*

_{13},

*f*

_{2}

*+ f*

_{10}

*+ f*

_{13},

*f*

_{2}

*+ f*

_{13}

*+ f*

_{16},

*f*

_{2}

*+ f*

_{13}

*+ f*

_{21},

*f*

_{5}

*+ f*

_{6}

*+ f*

_{7},

*f*

_{5}

*+ f*

_{6}

*+ f*

_{10},

*f*

_{5}

*+ f*

_{6}

*+ f*

_{12},

*f*

_{5}

*+ f*

_{6}

*+ f*

_{13},

*f*

_{5}

*+ f*

_{6}

*+ f*

_{16},

*f*

_{6}

*+ f*

_{10}

*+ f*

_{12},

*f*

_{6}

*+ f*

_{10}

*+ f*

_{16}, and

*f*

_{6}

*+ f*

_{10}

*+ f*

_{21}. The increased observation time made the features related to the horizontal component work very well, and the increased dimension provided additional information to separate the two classes further. For NNC

_{2},

*f*

_{2}

*+ f*

_{3}

*+ f*

_{16}satisfied the requirement; at SNR = 5 dB,

*P*

_{c}was only 95.2% due to the lack of generalization capability, so further study is required to achieve high

*P*

_{c}.

### IV. Conclusion

*v*and

*a*,

*P*

_{c}decreased considerably, so the effects of

*v*and

*a*should be removed before the features are extracted. The proposed features were also robust to the existence of multiple targets, yielding a small decrement in

*P*

_{c}. High

*P*

_{c}values were obtained by combining the features; as the observation time increased, the discriminant ability of the features related to the horizontal axis increased, so the improvement in

*P*

_{c}increased. NNC

_{1}yielded a higher

*P*

_{c}than NNC

_{2}, which is attributable to the simple structure of the neural network. Further studies should be conducted to increase the number of neurons or use convolutional neural network structures.